前言

最近因為工作需求碰到地理空間搜尋與儲存的技術，其中包含 PostGIS 這個 PostgreSQL 實作的地理物件函式庫，裏頭涵蓋各地區座標系統的 SRID(Spatial Reference Identifier)，令人欣慰的是台灣地區常用的 TWD97 與 TWD67 兩種座標系統也涵蓋如此，但在實作過程發現 TWD67 的座標點位有高達數百公尺的誤差，追溯源頭後意外發現，原來是 PostGIS 內的參數設定有誤，解決之後就順手紀錄一下踩坑的過程。

座標參考系格式 - WKT vs proj.4 vs SRID

開始以前要先提一下地理知識，因為地球是圓的不是平的，所以要將一個橢圓形的物件投影到紙面上勢必會有投影扭曲的問題，而為了降低投影扭曲的影響，學術上有很多轉換方式，包含等面積投影、等角投影(麥卡托投影)…，但不管選用哪一種投影都只能降低扭曲的幅度，並不存在絕對好的投影方式。

為了描述座標參考系(CRS, Coordinate Reference System)，訂了幾種常用的描述格式，包含:

  • WKT(Well-Known Text):

    • 是一種描述點、線、多邊形、幾何形狀的文本格式，例如: POINT(121.5 22.5)，相較於其他格式更 human-readable。

    • 舉例來說，Google map 使用的坐標系 WGS84 的 WKT 格式: GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

    • 在 PostGIS 中被應用在 spatial_ref_sys 表內的 srtext 欄位中

  • proj.4:

    • 是另一種描述投影的文本格式，例如: +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs，用幾個參數如橢圓角度、基準、長度單位…來定義投影的方式。
    • Proj.4 是開源社群常用的地圖投影資料庫，據我所知 PostGIS 和 GDAL/OGR 這些軟體都直接或間接的使用到 Proj.4 資料庫
    • 在 PostGIS 中被應用在 spatial_ref_sys 表內的 proj4text 欄位中

  • SRID:

    • SRID 是 Spatial Reference Identifier 的縮寫，代表某地區的座標參考系，是一串數字組成。例如: WGS84 的 SRID 是 4326
    • 而每個 SRID 都會有一個授權機構(authority)，以 PostGIS 為例，內建的 SRID 都是在於 EPSG 這個機構，類似的機構還有 ESRI
    • 可以在 spatialreference 網站上查到 EPSG 上有哪些 SRID
    • 在 PostGIS 中被應用在 spatial_ref_sys 表內的 srid 欄位中

想了解更多比較可以看這篇 文章

重現 TWD67 座標偏差問題

了解這些描述格式後，為了重現遇到的 TWD67 問題，我找了台北 101 景觀台的座標點來 demo。

打開 Google Map，把台北 101 的經緯度座標右鍵複製下來，經緯度是 25.03375, 121.56487，這邊的格式是 緯度,經度。

101-WGS84

接著，打開 成大水工試驗所座標轉換程式 網站，將剛剛複製的經緯度貼上來，並選擇 TWD97(WGS84)經緯度，按下轉換後可以得到 TWD67 二度分帶坐標值: 306171.309, 2769839.581

有了台北 101 景觀大樓在 TWD67 的座標後，我們來用 PostGIS 將他轉回 WGS84 看看 (TWD67 的 SRID 是 3828, WGS84 的 SRID 是 4326)

SELECT ST_Transform(ST_GeomFromText('POINT(306171.309 2769839.581)', 3828), 4326)

結果如下，左邊紅色圈圈是用 PostGIS 內建的 TWD67 轉換的座標，右邊圈圈是台北 101 景觀大樓真實位置。可以看到兩者差了不知道多遠…

TWD67-wrong

解決方法: 重新調整 TWD67 投影參數

為了解決這個問題，我直覺的先去看一下 PostGIS 內建的 TWD67 座標參數，用 SQL 查一下，TWD67 的 SRID 是 3828

SELECT * FROM spatial_ref_sys WHERE srid=3828;

查詢結果中 proj4text 欄位的值是 +proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs

這個 proj.4 文字描述的是 PostGIS 內 TWD67 坐標系如何做參數的轉換，上網查了一下發現其他座標轉換軟體如 QGIS 也有類似的問題，詳細可見 這篇。

研究了一陣子後，發現錯誤的原因是 PostGIS 的 proj4 投影參數少了towgs84，應該要修正成以下內容:

原本的 proj4 參數:

+proj=tmerc +lat_0=0 +lon_0=121 +k=0.9999 +x_0=250000 +y_0=0 +ellps=aust_SA +units=m +no_defs

應修改為:

+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

主要差異是:

+towgs84=-752,-358,-179,-0.0000011698,0.0000018398,0.0000009822,0.00002329

這其實就是七參數轉換法，包含 x, y, z 三向的平移與旋轉，再加上尺度縮放，所以共有七個參數，來決定如何將座標點位的參考系轉換成 WGS84。

OK，找到解法後，接下來就是想辦法修改 PostGIS 內的 proj4 內容即可，但我還是想保留原本 TWD67(3828)的座標參考系，所以選擇額外新增一個客製化的座標參考系，就選用最近沒被用到的 SRID 3830。

其中 srtext 的內容可以用另一個座標轉換軟體 GDAL 來生成

gdalsrsinfo.exe mycrs.prj > mycrs.txt

執行完後會產出 WKT 並寫進 mycrs.txt 內，複製出來貼到 PostGIS 內的 rtext 即可

'PROJCS["TWD67 / TM2 zone 121",GEOGCS["TWD67",DATUM["Taiwan_Datum_1967 based on Australian Natl & S. Amer. 1969 ellipsoid",ELLIPSOID["Australian Natl & S. Amer. 1969",6378160,298.25,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["TWD67",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",121,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9999,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",250000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

整體 SQL 如下:

/* Add custom TWD67' with 7 params transformation */
INSERT INTO spatial_ref_sys
(srid, auth_name, auth_srid, proj4text, srtext)
VALUES (
    3830,
    'custom-twd-67',
    3830,
    '+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'
    'PROJCS["TWD67 / TM2 zone 121",GEOGCS["TWD67",DATUM["Taiwan_Datum_1967 based on Australian Natl & S. Amer. 1969 ellipsoid",ELLIPSOID["Australian Natl & S. Amer. 1969",6378160,298.25,LENGTHUNIT["metre",1,ID["EPSG",9001]]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["TWD67",METHOD["Transverse Mercator",ID["EPSG",9807]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",121,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",0.9999,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",250000,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'
);

接著，實測客製化 TWD67(3830)的座標轉換

SELECT ST_Transform(ST_GeomFromText('POINT(306171.309 2769839.581)', 3830), 4326)

結果如下，可以看到座標轉換的結果比預設的 TWD67 好太多了!

twd67-correct

以上就是這次踩坑 PostGIS 的紀錄，沒想到依照國際地理標準的這些軟體工具，背後的投影參數也是有可能需要被調整的，幸好這些標準都依循既定的格式來描述，所以修改起來負擔不大，擴充上仍保有彈性。