GIS小知识-由对GCJ-02的疑问引出的坐标系

坐标系

在介绍GCJ-02之前，我们先了解一下常见关于坐标系的术语

直角坐标系

​ 直角坐标系又叫笛卡尔坐标系，它通过一对数字坐标在平面中唯一地指定每个点，该坐标系是以相同的长度单位测量的两个固定的垂直有向线的点的有符号距离。

平面坐标系

​ 表示点的平面位置。中国一般采用以高斯-克吕格投影分带的中央子午线为纵轴和赤道的投影为横轴的高斯-克吕格平面直角坐标系,简称高斯平面坐标系。坐标纵轴为x,自原点向北为正;坐标横轴为y,自原点向东为正。点的平面坐标为(x,y)。选任意子午线为坐标纵轴和高斯投影 (opens new window)面的坐标系或选高斯-克吕格投影 (opens new window)分带的中央子午线为纵轴和任意高程面的坐标系 (opens new window),则属于地方(矿区)平面坐标系。如果任意选定坐标原点和x轴方向,则称独立平面坐标系。

大地坐标系

地球并不是一个标准的球，而是一个南北两极稍扁赤道略胖的胖子。

大地水准面：用平静的海面描述地球就是假设地球表面都是水，且风平浪静。

旋转椭球面：椭球面在每个坐标平面上的投影都是椭圆,你可以用它的方程去验证。而旋转椭球面是可以用一个椭圆绕对称轴旋转得到,所以它在某个坐标平面上的投影是椭圆。

地心坐标系：以地球质心为旋转椭球面的中心的坐标系，又名协议地球坐标系，且唯一。

参心坐标系：即椭球中心不在地球质心的坐标系。（地图是为人民服务的，每个国家都想让地图尽可能准确地描述本国的地形地貌，因此就有国家把质心“移走”，让局部的表面“贴到”该国的国土，使之高程误差尽量减小到最小）

大地坐标系：（Geodetic coordinate System）是大地测量中以参考椭球面为基准面的坐标，地面点P的位置用大地经度L、大地纬度B和大地高H表示。地心坐标系、参心坐标系都是大地坐标系的一种。

GIS坐标系中的椭球，如果加上高程系，在其内涵上就是GCS（地理坐标系）。其度量单位就是度分秒。

中国常用大地坐标系

2000国家大地坐标系

2000国家大地坐标系，是我国当前最新的国家大地坐标系，英文名称为China Geodetic Coordinate System 2000，英文缩写为CGCS2000。2000国家大地坐标系的原点为包括海洋和大气的整个地球的质量中心；2000国家大地坐标系的Z轴由原点指向历元2000.0的地球参考极的方向，该历元的指向由国际时间局给定的历元为1984.0的初始指向推算，定向的时间演化保证相对于地壳不产生残余的全球旋转，X轴由原点指向格林尼治参考子午线与地球赤道面（历元2000.0）的交点，Y轴与Z轴、X轴构成右手正交坐标系。

西安80坐标系

西安80坐标系是指1980年西安坐标系，又简称西安大地原点。该坐标系的大地原点设在我国中部的陕西省泾阳县永乐镇，位于西安市西北方向约60公里，故称1980年西安坐标系，又简称西安大地原点，基准面采用青岛大港验潮站1952－1979年确定的黄海平均海水面（即1985国家高程基准）。西安80坐标系，属参心坐标系，长轴6378140m，短轴6356755，扁率1/298.25722101。

北京54坐标系

北京54坐标系(BJZ54)是指北京54坐标系为参心大地坐标系，大地上的一点可用经度L54、纬度M54和大地高H54定位，它是以克拉索夫斯基椭球为基础，经局部平差后产生的坐标系。1954年北京坐标系可以认为是前苏联1942年坐标系的延伸。它的原点不在北京而是在前苏联的普尔科沃。

WGS84坐标系

WGS84:World Geodetic System 1984，是为GPS全球定位系统使用而建立的坐标系统。通过遍布世界的卫星观测站观测到的坐标建立，其初次WGS84的精度为1-2m，在1994年1月2号，通过10个观测站在GPS测量方法上改正，得到了WGS84（G730），G表示由GPS测量得到，730表示为GPS时间第730个周。1996年，National Imagery and Mapping Agency (NIMA) 为美国国防部 (U.S.Departemt of Defense, DoD)做了一个新的坐标系统。这样实现了新的WGS版本：WGS（G873）。其因为加入了USNO站和北京站的改正，其东部方向加入了31-39cm 的改正。所有的其他坐标都有在1分米之内的修正。第三次精化：WGS84（G1150），于2002年1月20日启用。

投影坐标系

PCS(Projection Coordinate System)

在大地坐标系中是不能精准计算面积的，因为一度经度在不同的纬度表示的弧长是不一样的，在赤道附近弧长最长，在两极附近弧长最短。

为了建立一个坐标系可以使地图分析和空间分析可以定量计算，人们想到了投影坐标系。

地球是一个球面，而地图必须是一个平面，因此将地球表面展开成地图平面必然会产生裂隙或褶皱，必须采用一定的数学方法将曲面展成平面，而且使其变形较小，这就需要某种数学方法实现。

投影变换（地图投影）

​ 投影变换（projection transformation）是将一种地图投影点的坐标变换为另一种地图投影点的坐标的过程。研究投影点坐标变换的理论和方法。在地图投影的过程中，将球面上的元素（距离、角度、图形）投影到平面上，就会和原来的距离、角度、图形呈现出差异，这一差异就是投影变形。

投影变形的形式：角度变形、长度变形和面积变形。

地图投影的方式：

• 等角投影——投影前后的角度相等，但长度和面积有变形；
• 等距投影——投影前后的长度相等，但角度和面积有变形；
• 等积投影——投影前后的面积相等，但角度和长度有变形。

所以可以认为：PCS=GCS+投影方式

投影方式就是上面说的数学方法，不同的投影有不同的用途，也就有了不同的名字。

投影坐标系的基础是地理坐标系，没有地理坐标系，也就没有所谓的投影坐标系，投影坐标系是地理坐标系上的地物投射到具体投影面上的一种结果。

墨卡托投影

它的投影面是竖着的椭圆柱面，并且投影面与地轴方向一致，所以也叫正轴等角切/割圆柱投影。

意思就是既可以切圆柱，就是球体和椭圆柱面相切；也可以割圆柱，就是球体和椭圆柱面相割。

百度地图和Google Maps使用的投影方法都是墨卡托投影。

高斯-克吕格投影 高斯-克吕格投影是由德国数学家、物理学家、天文学家高斯于19 世纪20 年代拟定，后经德国大地测量学家克吕格于1912 年对投影公式加以补充，故称为高斯-克吕格投影，它的投影面是椭圆柱面，假设椭圆柱躺着，和地轴垂直，而且投影面与之相切，就是横轴墨卡托了，也可以称作等角横轴切椭圆柱投影。

竖着的有三条线，中间的那条就是投影中心线，根据取法不同，可以分为3度带和6度带。

3度带和6度带的起算经线是不一样的。

6°分带法:从格林威治零度经线起,每6°分为一个投影带,全球共分为60个投影带。

3°分带法:从东经1°30′起,每3°为一带,将全球划分为120个投影带。

高斯-克吕格投影的特点主要有三个：

• 投影后的地图，角度不变，面积会变。离中央经线越远的地区，面积变化越大。此投影合适用于导航。
• 投影椭圆柱面是横着的；
• 投影椭圆柱面与椭球体相切。

我国1：2.5万~1：50万地形图使用6度分带法；1：5000~1：10000地形图使用3度分带法。

通用横轴墨卡托投影

大家应该还听过它的江湖别称：UTM投影，它和高斯克吕格投影特别像，只不过它是割圆柱，就是球体与椭圆柱面相割。因此也被称作横轴等角割圆柱投影。

我国遥感影像通常会用UTM投影。

坐标系转换

​ 坐标转换是空间实体的位置描述，是从一种坐标系统变换到另一种坐标系统的过程。通过建立两个坐标系统之间一一对应关系来实现。是各种比例尺地图测量和编绘中建立地图数学基础必不可少的步骤。两个及以上的坐标转换时由极坐标相对参照确定维数空间。

七参数

​ 两个不同的三维空间直角坐标系（投影坐标系）之间转换时，通常使用七参数模型（数学方程组），在该模型中有七个未知参数。两个椭球间的坐标转换，一般而言比较严密的是用七参数布尔莎模型，即 X 平移， Y 平移， Z 平移， X 旋转（WX）， Y 旋转（WY）， Z 旋转（WZ），尺度变化（DM ）。要求得七参数就需要在一个地区需要 3 个以上的已知点。如果区域范围不大，最远点间的距离不大于 30Km（ 经验值 ） ，这可以用三参数，即 X 平移， Y 平移， Z 平移，而将 X 旋转， Y 旋转， Z 旋转，尺度变化面DM视为 0 。

火星坐标GCJ-02

了解了前置的坐标系知识后，这里我们终于可以开始讨论GCJ-02，又称火星坐标

火星坐标系，也叫国测局坐标系（GCJ02），正式名称为「地形图非线性保密处理技术」。是由中国国家测绘局2002年制订的地理信息系统的坐标系统。

我国规定国内出版的各种地图系统(包括电子形式)，必须至少采用“GCJ02”对地理位置进行首次加密。

所以在中国所有的地图必须使用“GCJ02”对地理位置进行首次加密。比如谷歌中国、高德、腾讯都在用这个坐标系。

它主要目的不是对GPS定位进行偏移，而是对高精度地图进行偏移，降低其精度。定位信息不是重中之重，高精度地图则是国家机密。而对定位进行偏移只不过是为了与脱密后的地图保持同步性，从而不影响导航等公开领域的应用。

目标是通过保密处理，对高精度数据进行偏移，使其成为低精度数据，从而降低国外导弹精准打击的准度，以保护我国、我军的重要基础设施。

GCJ-02坐标系 / 火星坐标系的存在没有现时意义，因为火星坐标系对内（测绘领域内部），偏移距离过大，导致某些数据产生的成果都是错的，也就阻碍了脱密数据的应用；对外，由于现阶段公式已被破解，实际上已丧失了对高精度地图的保密效果。

扩展文档 (opens new window)

java版坐标转换方法

WGS84、GCJ02、BD09坐标之间的转换方法

/**
 * 
 * 坐标系转换工具类
 */
public class GPSUtil {
	public static double pi = 3.1415926535897932384626;
	public static double x_pi = 3.14159265358979324 * 3000.0 / 180.0;
	public static double a = 6378245.0;
	public static double ee = 0.00669342162296594323;
 
	public static double transformLat(double x, double y) {
		double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y
				+ 0.2 * Math.sqrt(Math.abs(x));
		ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
		ret += (20.0 * Math.sin(y * pi) + 40.0 * Math.sin(y / 3.0 * pi)) * 2.0 / 3.0;
		ret += (160.0 * Math.sin(y / 12.0 * pi) + 320 * Math.sin(y * pi / 30.0)) * 2.0 / 3.0;
		return ret;
	}
 
	public static double transformLon(double x, double y) {
		double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1
				* Math.sqrt(Math.abs(x));
		ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
		ret += (20.0 * Math.sin(x * pi) + 40.0 * Math.sin(x / 3.0 * pi)) * 2.0 / 3.0;
		ret += (150.0 * Math.sin(x / 12.0 * pi) + 300.0 * Math.sin(x / 30.0
				* pi)) * 2.0 / 3.0;
		return ret;
	}
 
	public static double[] transform(double lat, double lon) {
		if (outOfChina(lat, lon)) {
			return new double[] { lat, lon };
		}
		double dLat = transformLat(lon - 105.0, lat - 35.0);
		double dLon = transformLon(lon - 105.0, lat - 35.0);
		double radLat = lat / 180.0 * pi;
		double magic = Math.sin(radLat);
		magic = 1 - ee * magic * magic;
		double sqrtMagic = Math.sqrt(magic);
		dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
		dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
		double mgLat = lat + dLat;
		double mgLon = lon + dLon;
		return new double[] { mgLat, mgLon };
	}
 
	public static boolean outOfChina(double lat, double lon) {
		if (lon < 72.004 || lon > 137.8347)
			return true;
		if (lat < 0.8293 || lat > 55.8271)
			return true;
		return false;
	}
 
	/**
	 * 84 to 火星坐标系 (GCJ-02) World Geodetic System ==> Mars Geodetic System
	 * 
	 * @param lat
	 * @param lon
	 * @return
	 */
	public static double[] gps84_To_Gcj02(double lat, double lon) {
		if (outOfChina(lat, lon)) {
			return new double[] { lat, lon };
		}
		double dLat = transformLat(lon - 105.0, lat - 35.0);
		double dLon = transformLon(lon - 105.0, lat - 35.0);
		double radLat = lat / 180.0 * pi;
		double magic = Math.sin(radLat);
		magic = 1 - ee * magic * magic;
		double sqrtMagic = Math.sqrt(magic);
		dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
		dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
		double mgLat = lat + dLat;
		double mgLon = lon + dLon;
		return new double[] { mgLat, mgLon };
	}
 
	/**
	 * * 火星坐标系 (GCJ-02) to 84 * * @param lon * @param lat * @return
	 * */
	public static double[] gcj02_To_Gps84(double lat, double lon) {
		double[] gps = transform(lat, lon);
		double lontitude = lon * 2 - gps[1];
		double latitude = lat * 2 - gps[0];
		return new double[] { latitude, lontitude };
	}
 
	/**
	 * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 将 GCJ-02 坐标转换成 BD-09 坐标
	 * 
	 * @param lat
	 * @param lon
	 */
	public static double[] gcj02_To_Bd09(double lat, double lon) {
		double x = lon, y = lat;
		double z = Math.sqrt(x * x + y * y) + 0.00002 * Math.sin(y * x_pi);
		double theta = Math.atan2(y, x) + 0.000003 * Math.cos(x * x_pi);
		double tempLon = z * Math.cos(theta) + 0.0065;
		double tempLat = z * Math.sin(theta) + 0.006;
		double[] gps = { tempLat, tempLon };
		return gps;
	}
 
	/**
	 * * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法 * * 将 BD-09 坐标转换成GCJ-02 坐标 * * @param
	 * bd_lat * @param bd_lon * @return
	 */
	public static double[] bd09_To_Gcj02(double lat, double lon) {
		double x = lon - 0.0065, y = lat - 0.006;
		double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * x_pi);
		double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * x_pi);
		double tempLon = z * Math.cos(theta);
		double tempLat = z * Math.sin(theta);
		double[] gps = { tempLat, tempLon };
		return gps;
	}
 
	/**
	 * 将gps84转为bd09
	 * 
	 * @param lat
	 * @param lon
	 * @return
	 */
	public static double[] gps84_To_bd09(double lat, double lon) {
		double[] gcj02 = gps84_To_Gcj02(lat, lon);
		double[] bd09 = gcj02_To_Bd09(gcj02[0], gcj02[1]);
		return bd09;
	}
 
	public static double[] bd09_To_gps84(double lat, double lon) {
		double[] gcj02 = bd09_To_Gcj02(lat, lon);
		double[] gps84 = gcj02_To_Gps84(gcj02[0], gcj02[1]);
		// 保留小数点后六位
		gps84[0] = retain6(gps84[0]);
		gps84[1] = retain6(gps84[1]);
		return gps84;
	}
 
	/**
	 * 保留小数点后六位
	 * 
	 * @param num
	 * @return
	 */
	private static double retain6(double num) {
		String result = String.format("%.6f", num);
		return Double.valueOf(result);
	}
 
	public static void main(String[] args) {
		String str = "113.2362898,21.9053547;113.2353993,21.9076491;113.2352312,21.9079071;113.2345464,21.9086146;113.2342782,21.9089879;113.2341989,21.9093424;113.2342173,21.9101219;113.2345064,21.9118989";
		String[] arr = str.split(";");
		String res = "";
		for(String s : arr) {
			String[] lntlat = s.split(",");
			Double lnt = Double.parseDouble(lntlat[0]);
			Double lat = Double.parseDouble(lntlat[1]);
			
			double[] bd09 = gps84_To_bd09(lnt, lat);
			res += bd09[0] + "," + bd09[1] + ";";
		}
		res = res.substring(0, res.lastIndexOf(";"));
		System.out.println(res);
		
	}
}