使用 proj 在 itrf 與 twd97 間轉換

日常我都是使用 proj 作座標轉換，從經緯度座標轉到平面座標或地心卡氏座標。最近因為使用了 network-rtk 相對定位，需要處理 gps 的 itrf 座標系統與台灣的 twd97 系統座標不同的問題，所以搞了幾天，終算弄懂怎麼用 proj 作不同 itrf 系統和與 twd97 系統的轉換。

itrf 框架概念

itrf 框架是一個由多個 gnss 固定站定義出來的座標系統， itrf 每幾年會依據新的觀觀資料，發布新版本，例如 itrf94 itrf2000 itrf2014 等。同一個點在 itrf 不同版本間座標會不同，同時每一個 itrf 版本也都是動態基準，同一個點在 itrf94 1996 年時和 itrf94 2002 年時座標也會不同。

一般這類座標的變化是每年數公厘等級，由於精度需求不高，其時間會以西元年加上小數點表示，例如 itrf94 在 2010 年，和 2020-04-01，差不多是 2020.25。

itrf 發布新版時，會提供新舊版間的轉換參數，每個轉換是不同系統在同一時間點的轉換。例如 itrf94 在 2010 年測得的座標，以 itrf 的參數轉到 itrf2000 後，也是 itrf2000 在 2010 年的座標。

itrf 的時變參數

itrf 提供的轉換是 helmert 轉換，其中每個參數是時間的線性函數，也就是不同系統間的差異，是會隨時間改變的。 itrf 是以每個參數每年會改變多少來提供轉換參數的，也就是一維線性擬合。

例如上述的 itrf 參數網站中，提供的是在 2010 時的 itrf2014 座標，要轉到 2010 時的 itrf94，其 Tx 參數會是 7.4。若要把 2020 時的 itrf2014 轉到 2020 的 itrf94，因為差了 10 年，Tx 變成 7.4 加上 Tx 的年變化率 0.1 乘上 10 年，也就 Tx 變成 8.4。

不同時間點的轉換

上述都是同一時間點在不同系統間的轉換，如果要做不同時間點的轉換，據 euref 提供的線上轉換功能，需要有該站在該系統中的速度。也就是說，同一點在座標系統中是會移動的，以台灣來說是板塊運動活躍地帶，約為一年 3 公分。

要注意的是，在不同系統上同一點的移動速度不同，因為不同系統間也是有相對速度的，不然轉換參數怎麼會依時間變化。

台灣的的座標系統

目前台灣最新的是 twd97 2010，也就是同樣使用 twd97 的座標系統 itrf94，但把時間點拉到 2010 年。

另外國土測繪中心提供的即時動態定位服務 e-gnss，可以在幾分鐘內達到數公分的精度。因為其使用的是 network-rtk，雖然可以獲得絕對座標，但本質上是相對於 egnss 固定站網的 相對定位 。所以 egnss 會有自己的一套座標系統，與 twd97 2010 不同，也與 gps 現今的座標系統 itrf2014 不同。故國土測繪中心有提供將 egnss 座標轉為 twd97 2010 的服務。

gps 的座標系統

現今大部份的 ppp 服務，結果是在 itrf2014 的瞬時時間點下；也就是 2020 年的觀測，解算使用 2020 年發布的精密星曆，所以解出的座標也會是 itrf2014 在 2020 時的座標。但早期 gps 與 itrf 的座標系統是不一致的，處理起來比較麻煩，我也沒有碰過。

proj 座標系統轉換

proj 在第 5 版提供了 cct 指令，並在隨附的檔案中提供了 itrf 各系統間的轉換參數，可以處理不同 itrf 系統間的轉換。其輸入輸出需要為地心卡氏座標。

可以將時間以年份放在座標的第四個欄位，或是用 -t 參數指令。例如從 itrf2014 2020 轉到 itrf94 2020，參數如下：

~ $ cct +init=ITRF2014:ITRF94 <<CCT
> -3025437.00     5035172.72 2476758.40 2020
> CCT
-3025437.0131   5035172.7359  2476758.3166     2020.0000

若要反向轉換，從 itrf94 轉到 itrf2014，加上 -I 參數即可， cct -I +init=ITRF2014:ITRF94 。

其中 itrf2014 到各系統的轉換參數，是使用 proj 的參數檔案儲存，在 linux 下是在 /usr/share/proj/ITRF2014 或 /usr/local/share/proj/ITRF2014 。

新竹市交通大學周遭的轉換

我在交大同時使用了即時 ppp 與 egnss 定位，由於二者存在偏差，即時 ppp 是在 itrf2014 瞬時座標所以打算進行轉換。

由於不知道 egnss 座標的座標系統，所以直接使用國土測繪中心提供的轉換服務，將 egnss 轉到 twd97 2010。

交通大學在 itrf 中的速度

再來需要讓 twd97 2010 和 itrf2014 瞬時能夠比較，我選擇將 itrf2014 轉到 twd97 2010。因為我手上有交大站的歷年座標，所以我可以算出交大站的速度，並假設區域內的速度接近，直接把交大站的速度當作是 ppp 定位結果位置的速度。因為在 2014 年之前的 gnss 座標所在座標系統我不確定懶得處理，所以直接用 2014-2020 的資料算出交大站在 itrf2014 中的速度。

並且為了減少誤差，我是計算在平面座標的速度，而不是地心卡氏座標的速度。因為平面座標中，高程方向幾乎沒有變，代表新竹地區的變化幾乎都是在平面上的位移，只有東西和南北方向有速度。但投影到地心卡氏座標的話，會造成 xyz 三方向都有位移。

具體轉換過程

ppp 結果為經緯度
投影到 twd97 平面座標
用東西方向和南北方向速度，把平面座標倒回 2010 年。
轉回經緯度
轉到地心卡氏座標
從 itrf2014 轉到 itrf94
轉回經緯度
投影到 twd97 平面座標
此時結果即是 twd97 2010 座標

寫成 proj 的 pipeline 指令的話，我是把速度寫成 helmert 中的 dx dy，（也就是 bash 中的 east north 變數，單位是每年變化的公尺數。）搭配上時間讓他自動轉回 2010。但因為 helmert 預設是用來轉地心卡氏座標的，直接把平面座標用 helmert 轉 proj 會認為你亂來，所以加一層 unitconvert 讓 proj 認為沒問題。

cct -d 5 -t 2020.25 +proj=pipeline \
    +step +init=epsg:3826 \
    +step +proj=unitconvert +xy_in=m +xy_out=m \
    +step +proj=helmert +dx=$east +dy=$north +t_epoch=2010 \
    +step +proj=unitconvert +xy_in=m +xy_out=m \
    +step +init=epsg:3826 +inv \
    +step +proj=cart \
    +step +init=ITRF2014:ITRF94 \
    +step +proj=cart +inv \
    +step +init=epsg:3826

或是可以分成多個指令，但缺點是會把一個指令可以做到的事拆到多個指令，每個指令都要重新格式化、解析一次數字，且每個指令都要加上 -f %.5f 指定精度， cct 則是用 -d 5 指定要到小數下五位。（這裡我懶的加。）

proj +init=epsg:3826 \
| cct -t 2020.25 +proj=helmert +dx=$east +dy=$north +t_epoch=2010 \
| cs2cs +init=epsg:3826 +to +proj=cart \
| cct -t 2020.25 +init=ITRF2014:ITRF94 \
| cs2cs +proj=cart +to +init=epsg:3826