太硬了，大家悠着點...

本篇講坐标系統的詳細定義，有關坐标系的變換公式，以及簡單說說高程坐标系統。

本文約6000字，閱讀時間建議45分鐘。硬内容比較多，如有疏漏錯誤請指出，建議有興趣的朋友進一步閱讀。

作者：部落格園/B站/知乎/csdn/小專欄 @秋意正寒

版權：轉載請告知，并在轉載文上附上轉載聲明與原文連結（https://www.cnblogs.com/onsummer/p/12082454.html）。

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.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

聊聊GIS中的坐标系|再版詳細定義、計算及高程系統
（子午圈（就是所在投影帶的中央經線圈）半徑）

t=tanB
聊聊GIS中的坐标系|再版詳細定義、計算及高程系統
p=180*3600/π
聊聊GIS中的坐标系|再版詳細定義、計算及高程系統
（子午線弧長）
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"]]

第一個屬性"CGCS2000 / 3-degree Gauss-Kruger CM 114E"是這個坐标系的名稱。
第三個屬性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.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.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從業者，能用上測量學和地圖學的坐标系統成果，已經遊刃有餘了。

我希望我的讀者也能明白這點，未來加油。

參考文檔

[1] 高斯正反算公式：https://wenku.baidu.com/view/5776611cd4d8d15abf234e14.html

[2] 資訊工程大學ppt：https://wenku.baidu.com/view/88fb6e0d84868762cbaed50d.html

[3] 扒一扒坐标轉換之七參數：http://www.sohu.com/a/318537831_689260

[4] 寫給測繪小白，講解四參數與七參數坐标轉換含義及差別：https://rtkhome.com/?p=1210

[5] 布爾莎-沃爾夫轉換模型的幾何證明：https://wenku.baidu.com/view/11bbf607ba1aa8114431d97f.html

聊聊GIS中的坐标系|再版 詳細定義、計算及高程系統

目錄