PROJ.4學習——坐標系轉換
前言
PROJ可以做任從最簡單的投影到許多參考數據非常復雜的轉換。PROJ最初是作為地圖投影工具開發的，但隨着時間的推移，它已經發展成為一個強大的通用坐標轉換引擎，可以同時進行大規模地圖投影和高精密度的坐標轉換。
在PROJ中，有兩個用於大地測量轉換的框架，proj框架和cs2cs框架。第一個是PROJ中用於進行大地測量轉換的原始且有限的框架，第二個是一個新添加的框架，旨在成為一個更完整的轉換框架。
在描述這兩個框架的細節之前，讓我們首先注意到，大多數大地測量轉換的情況都可以表示為一系列基本操作，一個操作的輸出是下一個操作的輸入。例如，當從UTM區域32，基准ED50，到UTM區域32，基准ETRS89時，在最簡單的情況下，必須經歷5個步驟：
-
- 1. 將UTM坐標反投影到地理坐標
- 2. 將地理坐標轉換為3D笛卡爾地心坐標
- 3. 應用Helmert轉換，從ED50轉換到ETRS89
- 4. 從笛卡爾坐標轉換回地理坐標
- 5. 最后將地理坐標投影到UTM 32區平面坐標
# ED50 / UTM zone 32N <23032> +proj=utm +zone=32 +ellps=intl +towgs84=-87,-98,-121,0,0,0,0 +units=m +no_defs <>
# ETRS89 / UTM zone 32N
<25832> +proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs <>
通過管道轉換
PROJ框架在做投影，坐標系轉換時，其操作風格類似於Linux Shell命令。管道框架是通過一個特殊的投影實現的，該投影以用戶提供的一系列基礎操作作為參數，並將這些操作串聯在一起，以實現所需要的完成轉換。
此外，一些基本的大地測量操作，包括Helmert轉換、一般的高階多項式位移和Molodensky模型轉換，都可以作為管道的一部分。 （詳見下面：基於方程式方法的坐標系轉換）
Molodensky變換直接從一個基准面的大地坐標轉換到另一個基准面的大地坐標，而Helmert變換(通常更准確)則從3D笛卡爾坐標轉換到3D笛卡爾坐標。
因此，當使用Helmert變換時，通常需要做一個從大地坐標到笛卡爾坐標的初始轉換，以及反過來的最終轉換，以得到期望的結果。幸運的是，這個三步復合變換具有吸引人的特性，即每一步只依賴於前一步的輸出。
因此，我們可以構建一個geodetic-to-geodetic Helmert轉換，通過捆綁在一起的輸出和輸入，3個步驟為：geodetic-to-cartesian→Helmert→cartesian-to-geodetic。
管道驅動程序通過這種鏈式轉換實現，該實現非常緊湊，每一步只包含一個偽投影，稱之為：pipeline。它以基本的投影字符串作為參數。
所有的管道（偽投影轉換）均由基礎的轉換構成，所有這些轉換均提供了框架，用於為廣泛的大地測量任務構建高精度的解決方案。
基於前言中的 案例，我們看看使用PROJ是如何實現 geodetic → Cartesian → Helmert → geodetic（大地坐標 → 秒迪卡 → 赫爾默特變換 → 大地坐標，即2-4步）的。
proj=pipeline
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
x=-81.0703 y=-89.3603 z=-115.7526
rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645
step proj=cart inv ellps=GRS80
完整的pipeline語法如下：
proj=pipeline
step init=./s45b.pol:s45b_tc32
step proj=utm inv ellps=intl zone=32
step proj=cart ellps=intl
step proj=helmert convention=coordinate_frame
x=-81.0703 y=-89.3603 z=-115.7526
rx=-0.48488 ry=-0.02436 rz=-0.41321 s=-0.540645
step proj=cart inv ellps=GRS80
step proj=utm ellps=GRS80 zone=33
這樣就實現了將一個分米級別轉換為厘米級別的操作。
pipeline這塊涉及到PROJ的4D轉換（cct）。在CentOs執行上面的命令嘗試。發現inv不支持等。也是很頭大，后續有機會再嘗試嘗試。
cs2cs表達式
| 參數 | 描述 |
| +datum | 基准面名稱，cmd中輸入：proj -ld |
| +geoidgrids | 用於垂直數據轉換的GTX網格文件的文件名 |
| +nadgrids | 用於數據轉換的NTv2網格文件的文件名 |
| +towgs84 | 3參數或7參數基面轉換 |
| +to_meter | 將水平單位轉換為米計算輸出轉換參數，如：1英尺= |
| +vto_meter | 將垂直單位轉換為米計算輸出轉換參數 |
cs2cs框架提供了管道框架中可用的大地測量轉換的子集。
坐標轉換在cs2cs框架中需要經過兩個轉換步驟。第一步，以WGS84為軸心為基准面，將需要轉換的數據轉換為WGS84；第二步，將轉換后的WGS84數據再轉換為特定需要的坐標系（通過利用Helmert赫爾默特變換或基准面改變，或兩者結合）。
基准面的改變可以通過“ proj-string ”中的 +towgs84，+nadgrids，+geoidgrids 參數來定義。
對於這三個參數都可以逆轉換，如果在“ proj-string ”中輸入指定的參數。
+towgs84 參數是三參數、七參數轉換（Helmert 赫爾默特變換）中，以WGS84作為中間轉換時，WGS84的參數。
如果沒有指定WGS84的具體實現參數，則在轉換的過程中會產生相當多的不確定性。
+nadgrids參數可以應用由校正網格插值得到的非線性平面校正。最初，這是作為轉換北美基准NAD27和NAD83之間的坐標的一種方法來實現的，但是可以對存在校正網格的任何基准進行校正。
+geoidgrids參數用於垂直組件中的網格矯正。
兩種網格校正方法都允許在同一轉換中包含多個網格
與transformation pipeline相比，cs2cs實際上是執行了兩個“ proj-string ”。源坐標系轉換為WGS84，WGS84再轉換為目標坐標系。兩個坐標系轉換，中間用 +to 連接。
# Greek GGRS87 基准面 轉 wgs84
cs2cs +proj=latlong +ellps=GRS80 +towgs84=-199.87,74.79,246.62 +to +proj=latlong +datum=WGS84 20 35 20d0'5.467"E 35d0'9.575"N 8.570
#EPSG文件提供 WGS72 轉 WGS84 的7參數
cs2cs +proj=latlong +ellps=WGS72 +towgs84=0,0,4.5,0,0,0.554,0.219 \
+to +proj=latlong +datum=WGS84
4 55
4d0'0.554"E 55d0'0.09"N 3.223
基於網格的基准面調整 Grid Based Datum Adjustments
在許多地方(特別是北美和澳大利亞)，國家大地測量組織提供網格移位文件，以便在不同的基准點之間進行轉換，如NAD27到NAD83。這些網格移位文件包括在每個網格位置應用的移位。實際上，網格移位通常是基於包含四個網格點之間的插值來計算的。
PROJ支持使用網格文件在不同的坐標系之間進行轉換。網格移位表的格式有ctable (PROJ nad2bin程序生成的二進制格式)、NTv1(舊的加拿大格式)和NTv2(.gsb -新的加拿大和澳大利亞格式)。
cs2cs和pipeline均可以用於網格校准轉換，只是cs2cs先轉換為WGS84，而pipeline只要存在網格，則可以任意轉換。
使用cs2cs轉換需要指定 +nadgrids 關鍵字。如：
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=ntv1_can.dat +to +proj=latlong +ellps=GRS80 +datum=NAD83 -111 50 111d0'2.952"W 50d0'0.111"N 0.000
ntv1_can.dat位於：E:\SvnWorkspace\LY_WEB_GIS\branches\Documents\ms4w-mapserver-for-wimdows\release-1911-x64-gdal-2-3-3-mapserver-7-2-1\bin\proj\SHARE 文件夾下面（根據你自己實際安裝位置），他是一個網格偏移文件，用於獲取所選點的網格位移值。
如果列出多個grid shift文件，在這種情況下，將依次嘗試每個文件，直到找到一個包含要轉換的點。
cs2cs +proj=latlong +ellps=clrk66 \
+nadgrids=conus,alaska,hawaii,stgeorge,stlrnc,stpaul +to +proj=latlong +ellps=GRS80 +datum=NAD83
-111 44
111d0'2.788"W 43d59'59.725"N 0.000
跳過丟失網格Skipping Missing Grids
@前綴加在網格文件前面，表示其可選。如果帶@符號的數據文件不存在，則搜索后面的文件。
通常任何未找到的網格都會導致錯誤。比如下面的案例會使用@ntv2_0.gsb文件，如果找到ntv2_0.gsb文件則使用，如果找不到則使用ntv1_can.dat。
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=@ntv2_0.gsb,ntv1_can.dat +to +proj=latlong +ellps=GRS80 +datum=NAD83 -111 50 111d0'3.006"W 50d0'0.103"N 0.000
空網格 The null Grid
PROJ在4.4.6版本后提供一個空網格文件，這個文件用於為世界提供一個零點轉換。如果您希望對所有其他網格的有效區域之外的點應用零位移，那么可以在+nadgrid文件列表的末尾列出它。通常，如果沒有找到包含要轉換的點的網格，將會發生錯誤。
cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null +to +proj=latlong +ellps=GRS80 +datum=NAD83 -111 45 111d0'3.006"W 50d0'0.103"N 0.000 cs2cs +proj=latlong +ellps=clrk66 +nadgrids=conus,null +to +proj=latlong +ellps=GRS80 +datum=NAD83 -111 44 -111 55 111d0'2.788"W 43d59'59.725"N 0.000 111dW 55dN 0.000
更多詳細信息，強查看：https://proj4.org/resource_files.html#transformation-grids
基於方程式方法的坐標系轉換
主要是介紹一下坐標轉換原理。
1. 三參數變換
基准面變換方法是地心（或三參數）變換。地心變換在 XYZ 或 3D 直角坐標系中對兩個基准面間的差異情況進行建模。定義一個基准面使其中心為 0,0,0。相距一定距離定義另一個基准面（dx,dy,dz 或 ΔX,ΔY,ΔZ，單位為米）。
通常，變換參數被定義為“從”區域基准面“到”1984 世界坐標系 (WGS) 或另一個地心基准面。
三個參數是線性平移量並且始終以米為單位。
2. 七參數轉換——Helmert 赫爾默特變換
從參考系1轉換到參考系2可以通過Δx，Δy，Δz三個偏移量，Rx，Rz，Rz三個旋轉角度和縮放比例參數 μ （圖二公式中是s）來實現。
赫爾默特變換(以弗里德里希·羅伯特·赫爾默特的名字命名，1843-1917)是三維空間中的一種變換方法。在大地測量中，它經常用於從一個基准到另一個基准的無失真轉換。Helmert變換也被稱為一個七參數變換，是一個相似變換。
這個就是空間坐標系中七參數的轉換原理。
3. Molodensky莫洛金斯基轉換
莫洛金斯基方法直接在兩種地理坐標系之間轉換，實際上無需轉換到 XYZ 系統。莫洛金斯基方法需要三個平移量 (dx,dy,dz) 以及兩個旋轉橢球體的長半軸 (Δa) 和扁率 (Δf) 的差。投影引擎根據相關基准面自動計算旋轉橢球體差。
- h = 橢球體高（米）
- Φ = 緯度
- λ = 經度
- a = 橢球體長半軸（米）
- b = 橢球體短半軸（米）
- f = 橢球體扁率
- e = 橢球體偏心率
M 和 N 分別是給定緯度下的子午線和卯酉圈曲率半徑。M 和 N 的方程如下：
求解 Δλ 和 ΔΦ。
簡化莫洛金斯基方法是莫洛金斯基方法的精簡版。請參見下面的方程：
基於網格的轉換
基於格網的變換方法包括下列方法：
NADCON 和 HARN 方法
美國使用基於格網的方法在多個地理坐標系之間進行轉換。基於格網的方法可用於對不同坐標系之間的差異進行建模，這些方法可能是最為精確的方法。感興趣區域被划分為多個單元。美國國家大地測量局 (NGS) 發布了用於在北美洲基准面 (NAD) 1927 (North American Datum (NAD) 1927) 和其他較早的地理坐標系與 NAD 1983 之間進行轉換的格網。這些變換均屬於 NADCON 方法。主 NADCON 格網（即 CONUS 格網）用於轉換美國本土毗鄰的 48 個州。其他 NADCON 格網用於針對以下地區將較早的地理坐標系轉換為 NAD 1983：
- 阿拉斯加
- 夏威夷群島
- 波多黎各和維爾京群島
- 阿拉斯加的聖喬治島、聖勞倫斯島和聖保羅島
美國本土毗鄰各州的精度約為 0.15 米，阿拉斯加及其所屬島嶼的精度約為 0.50 米，夏威夷的精度約為 0.20 米，波多黎各和維爾京群島的精度約為 0.05 米。上述精度可根據計算格網時區域中大地數據的准確度而發生變化（NADCON，1999）。
NAD 1927 中尚不存在夏威夷群島。夏威夷群島是使用被統稱為舊版夏威夷基准面的多個基准面繪制而成的。
NGS 和美國各州可以使用新的測量和衛星測量技術來更新大地控制點網絡。各州完成更新后，NGS 將發布用於在 NAD 1983 與更為精確的控制點坐標系之間進行轉換的格網。最初，此項工作被稱為“高精度大地網”(HPGN)。現在稱為“高精度參考網絡”(HARN)。截至 2004 年 1 月，美國的四塊領地和 46 個州已發布了 HARN 格網。HARN 變換的精度約為 0.05 米（NADCON，2000）。
以十進制秒為單位的差值存儲在兩個文件中：一個文件用於存儲經度差值，另一個用於存儲緯度差值。雙線性插值用於計算兩種地理坐標系在某一點處的精確差異。格網是二進制文件，但是，NGS 的 NADGRD 方案可將格網轉換為美國信息交換標准代碼 (ASCII) 格式。頁面底部顯示的是 CSHPGN.LOA 文件的標題和第一“行”。這是南加利福尼亞州的經度格網。第一行數字所表示的內容依次為：列數、行數、z 值數（始終為 1）、最小經度、單元大小、最小緯度、單元大小和未使用值。
此實例中隨后列出的 37 個值是以 0.25°（或 15 分）經度為間隔、在 32° N 位置處從 -122° 平移到 -113° 所得的經度位移值。
NADCON EXTRACTED REGION NADGRD 37 21 1 -122.00000 .25 32.00000 .25 .00000 .007383 .004806 .002222 -.000347 -.002868 -.005296 -.007570 -.009609 -.011305 -.012517 -.013093 -.012901 -.011867 -.009986 -.007359 -.004301 -.001389 .001164 .003282 .004814 .005503 .005361 .004420 .002580 .000053 -.002869 -.006091 -.009842 -.014240 -.019217 -.025104 -.035027 -.050254 -.072636 -.087238 -.099279 -.110968
國家坐標系變換第 2 版
與美國一樣，加拿大也使用基於格網的方法在 NAD 1927 與 NAD 1983 之間進行轉換。國家坐標系變換第 2 版 (NTv2) 方法與 NADCON 非常相似。兩種地理坐標系之間的差異包含在一組二進制文件中。雙線性插值用於計算點的精確值。
與每次只能使用一個格網的 NADCON 不同，NTv2 可檢查多個格網以獲得最精確的位移信息。存在一組針對加拿大的低密度基礎格網。某些區域（如城市）具有高密度局部子格網，它們與基礎格網（或父格網）的一部分相疊加。如果某點位於其中一個高密度格網中，NTv2 將使用高密度格網；否則，該點將位於低密度格網中。
NTv2 格網位移文件中的次格網插圖
如果某點位於上圖左下方的星形標記之間，則將使用高密度次格網計算位移。位於其他坐標位置的點將使用低密度基礎格網來計算其位移。軟件將自動計算所要使用的基礎格網或次格網。
加拿大的父格網的間距范圍從 5 分到 20 分不等。高密度格網的單元大小通常為 30 秒或 0.08333333°。
與 NADCON 格網不同，NTv2 格網將列出每個點的精度。精度值的范圍從幾厘米到大約一米不等。高密度格網的精度通常在厘米級以下。
澳大利亞和新西蘭也采用 NTv2 格式在不同的地理坐標系之間進行轉換。澳大利亞已發布了多種基於州的格網，用以在 1966 年澳大利亞大地基准面 (AGD 1966) 或 AGD 1984 與 1994 年澳大利亞地心基准面 (GDA 1994) 之間進行轉換。這些州格網已合並為國家格網。新西蘭已發布了用於在 1949 年新西蘭大地基准面 (NZGD 1949) 與 NZGD 2000 之間進行轉換的國家格網。
國家坐標系變換第 1 版
與 NADCON 一樣，國家坐標系變換第 1 版 (NTv1) 使用單個格網來對加拿大的 NAD 1927 與 NAD 1983 之間的差異進行建模。該版本就是 ArcInfo Workstation 中的 CNT。就精度而言，74% 的點的實際差異在 0.01 米以內，而對於所有情況的 93% 來說，實際差異則在 0.5 米以內。