最近內政部釋出 20 m DTM，志展(aka. Rex)收集了資料並做了許多處理應用(參見應用內政部20公尺網格數值地形模型資料)，綬草北三也提到幾個參數轉換的更正，一併整理如下：
MillerLiu's comments:
「關於坐標轉換參數請務必更正, 引述參數文章關於轉換部分恰好是我多年前蒐集及撰寫的， OSG1是TWD67轉TWD97、OSG2則是Hu-Tzu-Shan轉TWD97、OSG3則是Bessel轉TWD97， 三者是台灣自日治時期有測量以來使用的三種datum, 請務必說明補正, 不然參數誤用就完了, OSG1是成大教授論文所提供參數, 精度極高。OSG2則是中正理工學院論文提供, OSG3則是山友利用日治文獻座標計算而得。 因此僅有OSG1是目前多數應用需要的正確轉換參數。
上述網頁SINICA參數則是Garmin手持機中舊的三參數，已經在數年前由我去信修改成-685, -470, -237(台灣)，
及-752, -349, -179(澎湖)，此參數都經過共同點驗證過誤差，目前為所有Garmin 手持機的官方內建參數，
精度提升了約10米。OSG4則很有問題，與澎湖的三參數居然所差無幾。」
[original post](http://mutolisp.logdown.com/posts/207563-taiwan-geodetic-coordinate-system-conversion] 這是一個非常，非常，討厭而煩人的問題。但滿手歷史 GIS 資料的人，又不得不好好去面對這個歷史遺跡 XD
##原理 轉換的方式有很多種，如果是平面座標系統，就直接用平移的就好，但地球是個橢球體，地表又不平(XD)，所以球面座標系統轉換好麻煩，最常見的有三參數、四參數或是六參數、七參數的轉換，各有其優缺點我們就不一一介紹，這邊只稍微提一下 Bursa-Wolf 七參數(Bursa-Wolf transformation model)的轉換。簡單的說就是把原本的座標矩陣乘上轉換矩陣再加上各別維度(X, Y, Z)的偏移量，其數學計算式如下：
\left[ \begin{array}{lll}
X_t \\
Y_t \\
Z_t
\end{array} \right] =
\left[ \begin{array}{lll}
\Delta X \\
\Delta Y \\
\Delta Z
\end{array} \right] + \lambda \;
\left[ \begin{array}{rrr}
1 & \gamma & -\beta \\
-\gamma & 1 & \alpha \\
\beta & -\alpha & 1
\end{array} \right] \:
\left[ \begin{array}{lll}
X \\
Y \\
Z
\end{array} \right]
若以矩陣形式可以寫成
\mathbf{X}_{t} = \Delta\mathbf{X} + \lambda\mathbf{X_o}
其中 $X_t$ 為轉換過的三維座標系統(假設為 TWD1997 TM2)，$X_o$ 為原來的座標系統(假設為 TWD1967 TM2)，$\lambda$ 為尺度常數(scale factor)， $\alpha, \beta, \gamma$ 分別代表 x,y,z 旋轉的角度。求解的話可以用最小平方法(Least square approximation)來逼近參數值，但是需要有許多線性代數計算的技巧，所以就省略不提（要有 design matrix，transpose 來 transpose 去，算 residual errors 之類的）。圖示說明如下：
##實作 下表收集了一些常見的參數，而這些參數我們轉換的時候都會使用到：
| 來源 | dX (m) | dY (m) | dZ (m) | rX (radian) | rY (radian) | rZ (radian) | dS | 適用座標轉換 | 來源 |
|---|---|---|---|---|---|---|---|---|---|
| MOI | -750.739 | -359.515 | -180.510 | 0.00003863 | 0.00001721 | 0.00000197 | 0.99998180 | ||
| OSG1 | -730.160 | -346.212 | -472.186 | -0.00003863 | -0.0000172 | -0.00000197 | 0.99998180 | TWD67轉TWD97 | |
| OSG2 | 541.3141 | 48.2482 | 57.1666 | -0.0000140116 | 0.0001078381 | -0.0001818323 | -0.000002847 | 虎子山(Hutzushan)轉 TWD97 | |
| OSG3 | -1104.5 | 229.5 | 71.2 | -0.0000000485 | 0.0000000969 | 0.0000001939 | -0.000001 | Bessel 轉 TWD97 | |
| OSG4 | -752 | -358 | -179 | -0.0000011698 | 0.0000018398 | 0.0000009822 | 0.00002329 | *可能有問題 | |
| SINICA(臺灣本島) | -685 | -470 | -237 | ||||||
| SINICA(澎湖) | -752 | -349 | -179 |
MOI = 內政部86年二等衛星控制點測量平差工作報告; OSG[1-3] 為 osgeo wiki [6] 上的所列的參數，由 MillerLiu (綬草北三) 整理 OSG4 是以前 QGIS 1.8 之前附註在座標系統設定中說明的範例 SINICA 中研院的座標轉換程式參數（尚待查證）
不過如果你是要 TWD1997 TM2 和 TWD1967 TM2 互相轉換，若使用 QGIS 2.0 以上的版本內建的 EPSG Proj.4 參數，誤差則非常大，所以我們來看一下大概是怎麼回事： ###測試環境： Proj.4 版本 4.8.0
# 我們用基隆的一個衛星控制點來轉換，使用 cs2cs 這個程式來處理，
# 語法是cs2cs +init=原本的座標系統(EPSG:代碼) +to +init=要轉換過去的座標系統(EPSG:代碼）
$ echo "319685.630 2778228.552" | cs2cs +init=EPSG:3828 +to +init=EPSG:3826
319685.38 2778218.94 0.00
但和實際上的 TWD97 TM2 座標(320516.475 2778024.41)差距很大，因為 TWD1997 TM2 和 TWD1967 TM2 實際上大概直線距離會差 800 多公尺。所以我們來看一下其詳細的參數為何[5–6]：
TWD 1967 TM2 Zone 121 (EPSG:3828) 的 Proj.4 參數：
+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA
TWD 1997 TM2 Zone 121 (EPSG:3826) 的 Proj.4 參數：
+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
翻譯成一般人能比較理解的語言就是，平面投影座標系統是橫麥卡托投影(Transverse Mercator = tmerc a.k.a. TM)，原點的緯度是 0 度，中央經線是東經 121 度，尺度常數是 0.9999。因為台灣的橫麥卡托投影是二度分帶，由東經 121 度切成兩個區域，原點左方的 x 座標變成負值，我們希望調整成正值，因此東移 250000 m，+ellps 代表的是橢球體類型。比較了 TWD1997 和 TWD1967 的 Proj.4 參數，原來 EPSG:3828 (即 TWD1967 TM2 Zone 121) 並沒有設定座標系統參數的轉換！所以我們必須手動設置，就採用上表 OSG4 的參數（-752,-358,-179,-0.0000011698,0.0000018398,0.0000009822,0.00002329）
$ echo "319685.630 2778228.552" | cs2cs +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +towgs84=-752,-358,-179,-0.0000011698,0.0000018398,0.0000009822,0.00002329 +units=m +no_defs +to +init=EPSG:3826
320514.30 2778022.28 28.88
因此我們在 QGIS 裡頭的設定 -> 自定座標系統(Custom coordinate reference system definition) 中新增一個 TWD67 TM2 座標系統，並設定好參數，這樣轉換就能夠
+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +towgs84=-752,-358,-179,-0.0000011698,0.0000018398,0.0000009822,0.00002329 +units=m +no_defs
後來我用 R 寫了轉換，並且比較，內政部的就全台灣的控制點來說相對上是誤差比較小的，大概都誤差 2–4 m 內，但區域性可能還要再用 TLS/LS 來逼近參數
參考資料與延伸閱讀
- Bursa, M. (1962) The theory for the determination of the non-parallelism of the minor axis of the reference ellipsoid and the inertial polar axis of the Earth, and the planes of the initial astronomic and geodetic meridians from observations of artificial Earth satellites. Studia Geophysica et Geodetica, 6:209–214.
- Wolf, H. (1963) Geometric connection and re-orientation of three-dimensional triangulation nets. Bulletin Géodésique, 68:165–169.
- 中央研究院 GIS 應用支援工具集 #shptrans
- Spatial reference (EPSG) TWD1967 TM2 參數
- Spatial reference (EPSG) TWD1997 TM2 參數
- OSGEO wiki Taiwan datums
- Proj.4 general parameters
- http://spatialreference.org/ref/sr-org/7992/
- 應用內政部20公尺網格數值地形模型資料
- Facebook 討論
@mutolisp
gist有個小地方需要修正
根據proj.4的手冊，cs2cs 使用的init file名稱應為小寫(
epsg)，在一般環境下用大寫似乎會產生錯誤