1. 地理坐标系统定义

1.1. 人类对地球形状的描述

人类发现地球是个“赤道稍胖”的椭球后，就打算用一些数学或者物理的手段描述地球的形状。

早期，是用一个叫“大地水准面”的概念描述地球的，这个概念的说法是，地球海水静止后，海水面的形状就是地球的形状（陆地部分则想象海水穿过）。

后来，又提出了“似大地水准面”这一概念，它用的就不是海水面了，而是每个地方的重力线的顶点构成的面。

最后，为了便于数学计算，采用“椭球面”这一数学概念来描述地球形状。

在大地测量学中，“大地水准面”、“似大地水准面”所对应的“正高”、“正常高”是必须熟背于心的，但是在GIS中，本篇只讨论最后一个椭球面。

1.2. 旋转椭球面方程

根据解析立体几何，一个旋转椭球面的方程为：

它是个什么玩意儿呢？它是：

一个椭圆，这个椭圆以短轴为z轴，椭圆心为原点，然后绕z轴旋转而成的曲面。

(网络图片, http://xuxzmail.blog.163.com/blog/static/251319162009618102642971/ )

用平行于xOy平面的面切这个椭球面，相交的形状是一个圆。

1.3. 球面坐标系与经纬度

根据解析立体几何，常用三种三维空间坐标系，笛卡尔空间直角坐标系、球面坐标系、柱面坐标系。

本节回答为什么三维的经纬度只有两个分量的问题。

球面坐标系的定义是怎么样的呢？

球面坐标是三维坐标，自然有三个分量：r、θ、φ

r表示该点到原点的距离；θ表示该点与原点连线和z轴的夹角；φ表示该点与原点连线在xOy平面的投影和x轴的夹角。

那么，经纬度呢？

我们假想x轴是赤道面上这么一根半径所在的直线：这根半径线段与0度经线相交，也即：

同理，y轴、z轴也有类似的定义。

但是，点P（经度，纬度，第三个分量）究竟是什么呢？

其实，经度就是φ，纬度就是θ。

“经度(φ)就是椭球上的点与原点连线这一线段，在赤道平面（xOy平面）上投影与x轴的夹角”——只不过加了东经和西经，并不是0到360°。

“纬度就是椭球上的点与原点连线这一线段，与z轴的夹角的余角。”——赤道上的点与原点连线和z轴的夹角是90度，但是纬度是0度，所以是余角的关系。

所以，第三个分量就十分明确了：r，表示点到原点（椭球心）的距离。但是，为什么平时只用经纬度呢？

那是因为这个r非常大，通常我们谈高度只谈海拔高度，并不谈到地心的距离，所以这个r是被忽略的，这就解释了明明是三维坐标，却只有经纬度两个分量。

如果文字啃得太生硬，可以看下图：

1.4. 椭球与地理坐标系统

根据1.2，得知椭球面方程有两个参数a和b。

根据1.1，得知地球的形状是椭球体，表面是椭球面。

所以，描述地球通常只需要这两个参数即可，我们下一个定义：

定义a为赤道半径，即椭球的长半轴长；

定义b为极半径，即椭球的短半轴长。

赤道半径为地心（椭球心）到赤道任意一点的距离，极半径为地心（椭球心）到任意一个极点的距离。

有这两个参数后，还可以延伸出扁率和偏心率这两个概念。扁率有1个，偏心率则有两个。公式定义如下：

e和e'分别是第一偏心率和第二偏心率。

有了椭球，我们就有了地球的形状。实际上，地理坐标系统（GCS）的定义绝大部分就是由椭球体这两个参数定义的，那么地理坐标系统又是如何定义的呢？

给个公式吧：

GCS = f（椭球体）

f是椭球体的球心对于地球实际中心的偏移。为什么要做偏移？见下节讲解。

1.5. 参心地理坐标系统与地心地理坐标系统

根据1.4，我们知道地理坐标系统是定义在一个数学椭球面上的，具体方程已经给出。

但是还有一个小问题：偏移。

虽然椭球面方程“决定”了地球的形状，但是原点的位置却没有指定。按理说，是统一使用地心才对的，还是处于“懒”，为了方便计算，会直接使用椭球的球心当原点。

事实上，如果地心≠椭球心，椭球面就会比较靠近某个地区，这当然是认为的，这种“靠近”就便于某个国家或地区的计算，因为一旦靠近，很多地方的位置偏差就很小。

我们说，

地心地理坐标系统：椭球的球心=地球的质心

参心地理坐标系统：椭球的球心≠地球的质心

当今为了全球计量需要，有两个我们熟知的地心地理坐标系：WGS84和CGCS2000。

也就是说，北京54和西安80实际上是两个参心坐标系，它们的椭球体分别是克拉索夫斯基1940椭球体和IUGG1975椭球体。

1.6. WKT举例

还是老话，WKT的文章太多了，不再赘述，只摘取一些比较简单的属性讲解。

①WGS84

GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]

GEOGCS定义了一个地理坐标系统，内部第一个属性是字符串"WGS 84"是这个地理坐标系的名字。

然后，这个地理坐标系统有基准面"DATUM"，基准面下的"SPHEROID"是椭球体的意思，椭球体下的第二个、第三个属性是长半轴长和扁率的倒数。

最后AUTHORITY属性是这个地理坐标系的WKID信息，是4326.

②CGCS2000

GEOGCS["China Geodetic Coordinate System 2000",
    DATUM["China_2000",
        SPHEROID["CGCS2000",6378137,298.257222101,
            AUTHORITY["EPSG","1024"]],
        AUTHORITY["EPSG","1043"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4490"]]

和WGS84类似，不讲了。

1.7. 常见地理坐标系具体信息

这里不得不说的是，国家2000和WGS84几乎可以兼容，但是得先确定拿到的是真的国家2000的经纬度哦。

轶闻：其实还有一个新北京54坐标系的，WKID是4555，有兴趣的朋友可以查查这个坐标系的历史。

2. 投影坐标系统定义

2.1. 详细定义公式

PCS|x = f1（GCS|经纬度）

PCS|y = f2（GCS|经纬度）

简单解释一下：投影坐标系统的x坐标和y坐标分别由两个计算法则f1和f2计算，需要的参数有经度、纬度、椭球的参数。

2.2. 正算公式与反算公式

根据2.1，查阅资料，以4326做3857投影为例，以及CGCS2000做高斯克吕格投影为例。

不附代码。

① 网络墨卡托投影坐标系统

此处设网络墨卡托的地理坐标系统基于正球体，半径为R，点P的经纬度均为弧度十进制数：

x=R×弧度十进制经度

y=R×ln(tan(π/4 + 弧度十进制纬度/2))

此时，反算公式比较容易推导，不讲了。

② 高斯克吕格基于国家2000投影坐标系统

• 预备参数：椭球长半轴a；椭球扁率f；椭球短半轴b；椭球的第一第二偏心率e1、e2。

• 必备参数：经度J，纬度W

=====分割线=====

第一步，计算辅助量R、t、η、p、X、dL

• （子午圈（就是所在投影带的中央经线圈）半径）

• t=tanB
• 
• p=180*3600/π
• （子午线弧长）
• dL=B-中央经线度数

第二步，计算辅助常量a0、a2、a4、a6、a8和m0、m2、m4、m6、m8：

（这里e就是e1）

第三步，计算xy坐标：

反算公式即从x、y坐标算经纬度坐标。

此处不做展开，有兴趣的朋友可以查阅文末的参考文档。

2.3. 投影带问题

①换带操作

在arcgis中操作，其实只需要重投影即可。

一种方法是使用“投影”工具，将投影坐标系统的数据重新投影到它原本的地理坐标系统上，然后再用一次“投影”工具将地理坐标系统的数据再次投影到目标坐标系统上，完成换带。

另一种方法是直接用“投影”工具，将投影坐标系统的数据投影到目标PCS上即可。

具体操作见第4节。

②高斯克吕格投影坐标的判断

附一个坐标判断例子：

（41569821,4590855），已知在中国境内，已知地理坐标是国家2000.

横坐标是八位数，那么前两位一定是带号，41度带，那么就不可能是六度带，结果是三度带的高斯克吕格投影坐标系统，WKID是4529.

2.4. WKT举例

①网络墨卡托

PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["X",EAST],
    AXIS["Y",NORTH],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]

• 最外层是PROJCS，即投影坐标系统。
• 第一个属性"WGS 84 / Pseudo-Mercator"是这个坐标系的名称。
• 第二个属性GEOCS是这个投影坐标系统的地理坐标系统，详见上文。
• 第三个属性PROJCTION是投影方法"Mercator_1SP"。
• 第四~七个属性是其他属性，顺序下来是中央经线经度、比例因子、假东、假北。
• 第八个属性是单，第九个、第十个属性分别指示X和Y的方向是东和北。
• 第11个属性是此投影坐标系统在PROJ4中的定义。
• 第12个属性是此投影坐标系统在EPSG中的WKID。

②国家2000的高斯投影

以WKID=4547为例：

PROJCS["CGCS2000 / 3-degree Gauss-Kruger CM 114E",
    GEOGCS["China Geodetic Coordinate System 2000",
        DATUM["China_2000",
            SPHEROID["CGCS2000",6378137,298.257222101,
                AUTHORITY["EPSG","1024"]],
            AUTHORITY["EPSG","1043"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4490"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",114],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AUTHORITY["EPSG","4547"]]

• 最外层是PROJCS，即投影坐标系统。
• 第一个属性"CGCS2000 / 3-degree Gauss-Kruger CM 114E"是这个坐标系的名称。
• 第二个属性GEOCS是这个投影坐标系统的地理坐标系统，详见上文。
• 第三个属性PROJCTION是投影方法"Transverse_Mercator"，横轴墨卡托的意思。
• 第四~八个属性是其他属性，顺序下来是起始经线经度、中央经线经度、比例因子、假东、假北。
• 第九个属性是单位。
• 第十个属性是此投影坐标系统在EPSG中的WKID。

假东是什么意思？因为如果用赤道和中央经线的交点作为原点，投影得到的原始坐标会有负值。

我们记原始坐标为P，则给y坐标（经度方向）加500km后的P'就不会是负值了。

在P'的y坐标值（经线方向）加上带号，例如上图中的红色数字20，就成了带带投影带的坐标。

x方向的坐标一般不变，除非在地方坐标系中有需要，则设置假北（False North）。

2.5. 投影坐标系统的xy和ArcGIS的xy

在测量学的规定中，投影坐标系统上，x方向是指南北方向，y方向则是东西方向；

而在ArcGIS中，x方向则是东西方向，y方向是南北方向，正好颠倒。

所以，获取一份投影坐标系统的数据时，如果是正统的测量数据，那么y值应该在导入ArcGIS时被用于x，x值则用于y。

ps：我一直觉得，x和y只是一个记号，但是人就是那么喜欢用，换ab也可以，用uv也可以，切记：只是个符号，不要把xy的方向绝对化。

3. 高程系统

3.1. 1985国家高程基准

由1.3小节，我们知道球面坐标的第三个参数，点到椭球心（原点）的距离一般来说没什么用，我们听到更多的是“海拔高度”。

什么是海拔高度？珠穆朗玛峰海拔高度8844.43米，这个就是海拔高度。

那这个海拔高度的起点，也就是0米，是以那个地方的地面作为依据的呢？

答案就是，我国的“1985国家高程基准”，它的基准点位于青岛市某个地方，基准点高度为72.260m。

这个72.26m是什么意思呢？就是指，这个地方作为我国所有高程测量的起始值，别处测量的高度再加上72.26m即海拔高度。

3.2. GPS的高度——大地高度vs海拔高度

我们在文章开头1.1小节处，提及了正高和正常高两个概念，我们不引入太多测量学里的定义，但是，我国的高程系统一律是使用“正常高”的。

我们定义一个高度：大地高度H

大地高度是什么意思呢？大地高度就是点到椭球面的距离（沿着法线）。与1.3节定义的到椭球心的距离r相比，少了好长一截（暂停5分钟，读者可以想象一下）。

由卫星测得的高度就是大地高度。

那么大地高度H和我们说的海拔有什么关系呢？我们说我国高程测量是用“正常高”这个方法的，即重力等位面。

而我国的海拔高度又是基于正常高的，记作H'，那么H和H'的关系是：

H=H'+a

这里的a代表的意义是，“正常高”为零时的那个面距离椭球面的高度。回忆一下1.1节的内容，这个面是什么？

“正常高”的面是重力等位面，也即似大地水准面。

我们画个图表示表示：

当然，大地水准面也有类似的图：

此时大地高度H=H'+N，N即大地水准面到椭球面的距离，H'即正高（实际点到大地水准面的距离）。

plus：美国GPS的经纬度定位精度是不错的，但是高程的测量就比较差。

4. 坐标系统转换

4.1. n参数（n=3，4，7）与地理转换

①n参数

一个坐标系统挪到另一个坐标系统，有哪些情况呢？

最简单的是平移原点，只需要给出三个方向的平移量dx、dy、dz，此时，称之为三参数转换；

复杂的还可以加上4个量：三个方向的旋转角度α、β、γ+统一的缩放比例k，称之为七参数转换；

另外，如果是平面上二维坐标系的转换，可以使用两个平移量dx、dy，一个旋转角度α，一个缩放比例k来完成。

举个例子，在珠海既有基于北京54的投影坐标系统又有珠海的自己的地方投影坐标，在这两种坐标之间转换就用到四参数。四参数的获取需要有两个公共已知点。

如果区域范围不大，最远点间的距离不大于 30Km（ 经验值 ） ，这可以用三参数或者四参数。

坐标系统转换的实质就是地理坐标系统的转换，也即椭球体的转换。

当然，在书本上，会有投影坐标系统直接转换而不经过地理坐标系统的算法（《地理信息系统概论》黄杏元第三版），但是那个比较难。

②地理转换

在ArcGIS中，允许用户自定义七参数或三参数来进行不同椭球体（不同地理坐标系）的转换，当然，这些所谓的七参数和三参数的获取，至少在国内的转换中，是保密的，需要到有关部门购买相同位置的三个点的两个不同坐标系下的坐标，然后自己计算得到七参数。

有关这些参数的计算，参考更丰富的测量专业的书籍或者博客。

假设已经获取了七参数/三参数，那么可以在ArcMap中，使用“创建自定义地理(坐标)转换”工具为这些参数定义一个“地理转换”：

方法参数有很多，选一个需要的即可，不懂是啥的可以百度一下（我也没用过，大家可以边搜边试）。

4.2. ArcGIS中重投影操作

使用“地理转换”工具和“投影”/“投影栅格”工具。以下以矢量数据为例，使用“投影工具”。

①PCS1转PCS2（不同GCS）（使用投影工具）

跨不同地理坐标系统的转换，需要使用4.1提及的自定义地理(坐标)转换工具创建地理转换。

②PCS1转PCS2（相同GCS）（使用投影工具）

③PCS1回算PCS1.GCS（使用投影工具）

④GCS1转GCS2

两个不同地理坐标系的数据进行坐标系转换，需要使用4.1提及的自定义地理(坐标)转换工具创建地理转换：

此处为WGS84到国家2000，椭球不同，必须使用地理转换。

我们发现，需要地理转换的操作，通常就意味着跨地理坐标系统转换；

反过来说，跨地理坐标系统的转换就需要一个地理转换定义，也即n参数。

4.3. 前端转换计算之turf.js

turf.js只支持3857和4326的互转。

①使用turf.toWgs84()转换网络墨卡托的xy坐标到经纬度

②使用turf.toMercator()转换经纬度到xy网络墨卡托坐标

4.4. 前端转换计算之openlayers(6.x)

主要功能都在ol/proj模块下，另外在自定义坐标系和转换时会用到第三方库proj4.js，但本文非开发类的博客，不细展开。

①ol/proj.fromLonLat(coordinate, opt_projection)方法

fromLonLat方法将经纬度coordinate转换到目标坐标系opt_projection下，opt_projection默认值是"EPSG:3857"，是“ProjectionLike”类型的参数。

对应方法是ol/proj.toLonLat()。

②ol/proj.get(string)

获取坐标系信息，string是"EPSG:3857"的字符串，必须大写EPSG。这个字符串在openlayer6中叫做“ProjectionLike”类型。

返回一个ol/proj/Projection类型的对象

③ol/proj.addCoordinateTransforms(source, destination, forward, inverse)

添加两个坐标系之间的转换方法，source是待转换坐标系，destination是目标坐标系，二者均以"EPSG:XXXX"的字符串传入。

forward是

④ol/proj.proj4.register(proj4)

让openlayer知道你注册了一个自定义坐标系统。详情请参考proj4.js有关资料。

⑤ol/proj.getTransform(source, destination)

给定待转换坐标系source和目标坐标系destination，返回二者之间的转换方法。

⑥ol/proj.transform(coordinate, source, destination)

将坐标点从source坐标系到destination坐标系转换，source和destination均为"EPSG:xxxx"的字符串（即“ProjectionLike”类型），EPSG四个字母大写。

4.5. 前端转换计算之cesium

cesium只支持4326和3857的互相转换。常用的类有如下几个：

①Cesium.MapProjection类

属性：

ellipsoid。Ellipsoid类型，即椭球。

方法：

project()和unproject()。一个用于将地理坐标转换为投影坐标，一个用于将投影坐标转回地理坐标。详见API。

②Cesium.GeographicProjection(ellipsoid)类

表示地理坐标系统的一个类，使用Ellipsoid类型的参数进行实例化。方法与MapProjection类相同。

默认构造参数是Ellipsoid.WGS84

③Cesium.WebMercatorProjection(ellipsoid)类

表示网络墨卡托投影坐标系统的一个类，使用Ellipsoid类型的参数进行实例化。

默认构造参数是Ellipsoid.WGS84（是不是很奇怪，和上面那个一样）

也拥有project()和unproject()两个方法。详见API。

④Cesium.Cartographic(longitude, latitude, height)类

这个类的意思就是一个地理坐标系统下的点，包括经度longitude，纬度latitude，和大地高度height

静态方法：

• Cesium.Cartographic.fromCartesian(Cartesian3对象, ellipsoid, result)：将投影坐标实例Cartesian3转换到地理坐标系统ellipsoid上，通常ellipsoid参数是Ellipsoid.WGS84。
• Cesium.Cartographic.fromDegrees(经度,纬度,大地高度,result)：创建一个地理坐标点
• Cesium.Cartographic.fromRadians()：同上只不过用弧度制
• Cesium.Cartographic.toCartesian()：将地理坐标转换为投影坐标

⑤Cesium.Cartesian3(x, y, z)类

笛卡尔坐标点，即投影坐标点。

该类也提供了类似Cartographic类的转换方法，详情请自行查阅API文档。

4.6. *硬改数据坐标系的定义

在gis软件中为数据重新定义一个坐标系，这有可能出现极大问题。通常不推荐做这种非精确的转换。

曾经在实践中遇到过类似的问题，就是很多情况下，有的人并不在意坐标系有多么精确，甚至有时候，能把数据强硬编辑挪到喜欢的位置上就罢了。

事实上，在精度不高的情况下（例如一个城市，或者一个城市群这么大级别的区域），直接改动数据的坐标系统的定义，而不是经过精确的地理转换、坐标转换计算，有时候在这么大的尺度下可能看不出来什么。

有个特例，WGS84和国家2000坐标系的改动——因为这两个坐标系的的确确很接近。什么？你跟我说硬改还是很大偏差？

那你考虑一下你是否拿到了真的国家2000坐标，而不是什么所谓的GCJ02和BD09。

碎碎念

又熬夜了，能在2019年结束前重写完坐标系这三篇博客，也算是对自己的一个承诺的实现了。

我知道在大地测量学专业上有更加精妙的计算，有更为严苛的定义和转换，但是，作为一个GIS从业者，能用上测量学和地图学的坐标系统成果，已经游刃有余了。

我希望我的读者也能明白这点，未来加油。