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

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

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

作者：博客園/B站/知乎/csdn/小專欄 @秋意正寒

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

目錄

1. 地理坐標系統定義

2. 投影坐標系統定義

3. 高程系統

4. 坐標系統轉換

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從業者，能用上測量學和地圖學的坐標系統成果，已經游刃有余了。

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

 

參考文檔

[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