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

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

本文約<>字，閱讀時間建議45分鐘。

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

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

目錄

1. 地理座標系統定義

2. 投影座標系統定義

3. 高程座標系統

4. 座標系統轉換

1. 地理座標系統定義

1.1. 人類對地球形狀的描述

人類發現地球是個「赤道稍胖」的橢球后，就打算用一些數學或者物理的手段描述地球的形狀。

早期，是用一個叫「大地水準面」的概念描述地球的，這個概念的說法是，地球海水靜止後，海水面的形狀就是地球的形狀（陸地部分則想象海水穿過）。

後來，又提出了「似大地水準面」這一律念，它用的就不是海水面了，而是每一個地方的重力線的頂點構成的面。

最後，爲了便於數學計算，採用「橢球面」這一數學概念來描述地球形狀。

在大地測量學中，「大地水準面」、「似大地水準面」所對應的「正高」、「正常高」是必須熟背於心的，可是在GIS中，本篇只討論最後一個橢球面（第3節討論高程時會簡單說）。

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；橢球的第一第二偏愛率e一、e2。

必備參數

經度J，緯度W

=====

第一步，計算輔助量R、t、η、p、X、dL

子午圈（就是所在投影帶的中央經線圈）半徑R

 

 

t=tanB

 

 

p=180*3600/π

（子午線弧長）

 

 

dL=B-中央經線度數

 第二步，計算輔助常量a0、a二、a四、a六、a8和m0、m二、m四、m六、m8：

 

 

 

 

 

 第三步，計算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. 高程座標系統【未寫完】

我國的高程座標系統，瞭解瞭解就能夠了

通常用的是黃海85

要寫一下GPS高度、大地高度、正高、正常高的含義。

4. 座標系統轉換

咱們在座標系統轉換的時候，用的圖形學變換方法是仿射變換。在三維空間中，進行座標系統的轉換就至關於進行了一次向量空間變換，須要一個轉換矩陣。

座標系統轉換的實質就是地理座標系統的轉換。

固然，在書本上，會有投影座標系統直接轉換而不通過地理座標系統的算法（《地理信息系統概論》黃杏元第三版），可是那個比較難。

4.1. 轉換矩陣與n參數【未寫完】

引入布爾莎模型等轉換模型

4.2. ArcGIS中重投影操做【未寫完】

使用「地理轉換」工具和「投影」/「投影柵格」工具。

①PCS1轉PCS2（不一樣GCS）

②PCS1轉PCS2（相同GCS）

③PCS1回算PCS1.GCS

④GCS1轉GCS2

咱們發現，須要地理轉換的操做，一般就意味着跨地理座標系統轉換，反過來講，跨地理座標系統的轉換就須要一個地理轉換定義，也即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下，默認值是"EPSG:3857"。

對應方法是ol/proj.toLonLat()。

②ol/proj.get(string)

獲取座標系信息，string是"EPSG:3857"的字符串，必須大寫EPSG。

返回一個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"的字符串，EPSG大寫。

4.5. 前端轉換計算之cesium【未寫完】

cesium只支持4326和3857的互相轉換。經常使用的類有以下幾個：

①Cesium.MapProjection類

②Cesium.GeographicProjection(ellipsoid)類

③Cesium.WebMercatorProjection(ellipsoid)類

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

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

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