GIS坐标系、投影和转换

整理资料信息如下：

坐标转换工具：几个常用坐标系，在投影坐标和高斯坐标之间相互转换。

空间直角坐标系转换：使用七参数转换不同的坐标系。

GIS坐标系和ArcGIS坐标系介绍：GIS地理坐标系（大地坐标系）、投影坐标系（主要是高斯克吕格投影、测量坐标系）、地理坐标和投影坐标之间转换。同一个坐标系的经纬度坐标和投影坐标之间可以直接转换，不需要其它参数。

高斯克吕格投影：高斯克吕格投影坐标系中X轴和Y轴的定义，投影坐标值的构成方式、3度分带和6度分带的划分方式。

Yaw（偏航角）、Pitch（倾斜角）、Roll（旋转角）：三维坐标系中这个三个角度的定义。

 

Java使用七参数做坐标变换：

package gis.util;

public class XYTransformor {
	
	private double dx = 0 ;
	
	private double dy = 0 ;
	
	private double dz = 0 ;
	
	private double wx = 0 ;
	
	private double wy = 0 ;
	
	private double wz = 0 ;
	
	private double m = 0 ;
	
	private double targetX = 0 ;
	
	private double targetY = 0 ;
	
	private double targetZ = 0 ;
	
	/**
	 * 设置坐标转换的七参数
	 * @param dx ：x方向平移量
	 * @param dy ：y方向平移量
	 * @param dz ：z方向平移量
	 * @param wx ：x方向旋转量
	 * @param wy ：y方向旋转量
	 * @param wz ：z方向旋转量
	 * @param m ：比例因子
	 */
	public XYTransformor(double dx,double dy,double dz,double wx,double wy,double wz,double m){
		this.dx = dx ;
		this.dy = dy ;
		this.dz = dz ;
		this.wx = wx ;
		this.wy = wy ;
		this.wz = wz ;
		this.m = m ;
	}
        //........set and get method	
	public boolean transform(double x,double y ,double z){
		targetX = dx + (1+m)*x + wz*y - wy*z ;
		
		targetY = dy + (1+m)*y - wz*x + wx*z ;
		
		targetZ = dz + (1+m)*z + wy*x - wx*y ;
		
		return true ;
	}

}

 坐标投影、距离计算、三维坐标角度计算、经纬度与弧度相互转换帮助类：

package gis.util;

public class GISUtil {
	/**
	 * 默认地球长半轴长度
	 */
	public static double EARTH_RADIUS = 6378140.0;
	/**
	 * 弧度转度参数
	 */
	public static double RADIAN_CONSTANT = 180.0/Math.PI ;
	/**
	 * 度转弧度参数
	 */
	public static double ANGLE_CONSTANT = Math.PI/180.0 ;
	/**
	 * 计算地球上两点之间的实际距离，单位米
	 * @param longitude：第一个点经度，单位度
	 * @param latitude：第一个点纬度，单位度
	 * @param longitude2：第二个点经度，单位度
	 * @param latitude2：第二个点纬度，单位度
	 */
	public static double getGeographicDistance(double longitude,double latitude,double longitude2,double latitude2,double semimajorAxis){
		double tempLongitude = angleToRadian(longitude);
		double tempLatitude = angleToRadian(latitude);
		double tempLongitude2 = angleToRadian(longitude2);
		double tempLatitude2 = angleToRadian(latitude2);
		double angleValue = Math.sin(tempLatitude)*Math.sin(tempLatitude2)+Math.cos(tempLatitude)*Math.cos(tempLatitude2)*Math.cos(Math.abs(tempLongitude-tempLongitude2)) ;
		return semimajorAxis * Math.acos(angleValue) ;
	}
	/**
	 * angle to radian
	 * @param angle
	 * @return the length ,unit is meter
	 */
	public static double angleToRadian(double angle){
		return angle * ANGLE_CONSTANT ;
	}
	/**
	 * 获取偏航角
	 * @param x：单位米
	 * @param y：单位米
	 * @param x2：单位米
	 * @param y2：单位米
	 * @return
	 */
	public static double getYawAngle(double x, double y, double x2, double y2){
		return Math.atan2(y2-y,x2-x) * RADIAN_CONSTANT ;
	}
	/**
	 * 获取倾斜角
	 * @param x：第一个点经度，单位度
	 * @param y：第一个点纬度，单位度
	 * @param z：第一个点高程，单位米
	 * @param x2：第二个点经度，单位度
	 * @param y2：第二个点纬度，单位度
	 * @param z2：第二个点高程，单位米 
	 * @return
	 */
	public static double getPitchAngle(double x, double y, double z, double x2, double y2, double z2,double semimajorAxis){
		double length = getGeographicDistance(x,y,x2,y2,semimajorAxis);
	    return Math.atan2(-length,z2-z)* RADIAN_CONSTANT ;
//		return Math.atan2(z2-z,length)* RADIAN_CONSTANT ;
	}
	/**
	 * @param longitude：经度 <br>
	 * @param latitude：纬度 <br>
	 * @param zoneWide：经度分带类型 <br>
	 * @param semimajorAxis：地球长半轴长度，单位米
	 * @param inverseFlattening：扁心率的倒数
	 */
	public static double[] gaussKrugerProject(double longitude, double latitude,int zoneWide,double semimajorAxis,double inverseFlattening){
		//默认为6的分带
		if(zoneWide != 3&&zoneWide != 6){
			zoneWide = 6 ;
		}
//		//计算所在位置分带的中央经线
		double primeMeridian = 0 ;
		int projNo= (int)((longitude-1.5) / zoneWide) ; 
		if(zoneWide==3){
			primeMeridian = projNo * zoneWide + zoneWide/2.0 + 1.5; 
		}else{
			primeMeridian = projNo * zoneWide + zoneWide/2.0 ; 
		}
		//度转弧度
		primeMeridian = angleToRadian(primeMeridian);
		longitude = angleToRadian(longitude);
		latitude = angleToRadian(latitude);
		//第一扁心率
		double e2=(2-1/inverseFlattening)/inverseFlattening ;
		//第二扁心率
		double ee=e2*(1.0-e2);
		
		double N=semimajorAxis/Math.sqrt(1.0-e2*Math.sin(latitude)*Math.sin(latitude));
		double T = Math.pow(Math.tan(latitude),2);
		double C = ee*Math.pow(Math.cos(latitude),2);
		double A=(longitude-primeMeridian)*Math.cos(latitude);
		double M=semimajorAxis*((1-e2/4-3*e2*e2/64-5*e2*e2*e2/256)*latitude-(3*e2/8+3*e2*e2/32+45*e2*e2
				*e2/1024)*Math.sin(2*latitude)+(15*e2*e2/256+45*e2*e2*e2/1024)*Math.sin(4*latitude)-(35*e2*e2*e2/3072)*Math.sin(6*latitude));
		double x = N*(A+(1-T+C)*A*A*A/6+(5-18*T+T*T+72*C-58*ee)*A*A*A*A*A/120);
		double y = M+N*Math.tan(latitude)*(A*A/2+(5-T+9*C+4*C*C)*A*A*A*A/24 + (61-58*T+T*T+600*C-330*ee)*A*A*A*A*A*A/720);
		
		//正北方向为x轴，正东方向为y轴，为了避免x轴出现负值，将坐标原点向西偏移500千米，并添加前缀分带编号
		x += 1000000L*(projNo+1)+500000L;
		return new double[]{y,x} ;
	}
}