第二十四個夏天後: 常用的座標轉換筆記(TWD67, TWD97, WGS84)

台灣現在常看到的座標系統有三種，分別是 TWD67、TWD97 和 WGS84，其中 WGS84 就是網路應用上常用的 GPS (lat, lon) 格式(歐洲近幾年有在推另一個座標系統)。至於為啥會接觸到呢？實在是這幾年 ITS 的服務越來越興盛，而台灣區的資料慢慢地釋出(如 Taipei Open Data、交通部運輸研究所等)，然而就開始碰到資料格式問題，那就是台灣不少資料採用的地理座標系統不是 GPS 格式。如果你不是 ITS 領域，那大概就跟我一樣從摸 Google Maps 而接觸 GPS 座標系統，例如只要在 Google Maps 上打上 25.033661, 121.564815 後，顯示的就是 Taipei 101，但從台灣區取得的地理座標系統，丟到 Google Maps 卻行不通 :P

這幾年都有接觸 ITS 應用，但一直很懶得整理導致每次都花差不多的時間去處理，所以打算來寫篇筆記

目前碰過的座標格式：

TWD67 (TM: 2-degree Transverse Mercator, 二度分帶)

採用 1967年國際地球原子參數(Geodetic Reference System 1967，GRS67)，早期還未有衛星時，透過天文觀測、三角定位量測，埔里虎子山為測量原點，而 TWD 全名為 Taiwan Datums

TWD97

採用 1980年國際地球原子參數(Geodetic Reference System 1980, GRS80)，以 GPS 衛星定位重新測量，國內是在 1997 年啟用，所以稱作 TWD97

WGS84

TWD67 與 TWD97 轉換：

從 TWD97 轉 TWD67

X67=X97-807.8-A*X97-B*Y97
Y67=Y97+248.6-A*Y97-B*X97
A=0.00001549, B=0.000006521

從 TWD67 轉 TWD97

X97=X67+807.8+A*X67+B*Y67
Y97=Y67-248.6+A*Y67+B*X67
A=0.00001549, B=0.000006521

TWD67 轉 WGS84 (GPS)：

$ echo TWD67_X TWD67_Y | proj -I +proj=tmerc +ellps=aust_SA +lon_0=121 +x_0=250000 +k=0.9999

TWD97 轉 WGS84 (GPS)：

$ echo TWD97_X TWD97_Y | proj -I +proj=tmerc +ellps=GRS80 +lon_0=121 +x_0=250000 +k=0.9999

目前測試的心得：

先把 TWD67 轉 TWD97 後，在用 TWD97 轉 WGS84 的成果跟 Google Maps 的標記就接近了不少，但仍差大概約25m的距離，也還有一條街的距離。

Python Code：

def WGS84FromTWD67TM2(x,y):
out = {'status':False}
lat = None
lon = None

# TWD67 to TWD97
A = 0.00001549
B = 0.000006521
x = float(x)
y = float(y)
x = x + 807.8 + A * x + B * y
y = y - 248.6 + A * y + B * x

# TWD97 to WGS84
result = os.popen('echo '+str(x)+' '+str(y)+' | proj -I +proj=tmerc +ellps=GRS80 +lon_0=121 +x_0=250000 +k=0.9999 -f "%.8f"').read().strip() # lat, lng 格式, 不必再轉換
process = re.compile( '([0-9]+\.[0-9]+)', re.DOTALL )
for item in process.findall(result):
if lon == None:
lon = float(item)
elif lat == None:
lat = float(item)
else:
break

# result = os.popen('echo '+str(x)+' '+str(y)+' | proj -I +proj=tmerc +ellps=GRS80 +lon_0=121 +x_0=250000 +k=0.9999').read().strip() # 分度秒格式

# 分度秒格式轉換
#process = re.compile( "([0-9]+)d([0-9]+)'([0-9\.]+)\"E\t([0-9]+)d([0-9]+)'([0-9\.]+)", re.DOTALL )
#for item in process.findall(result):
# lon = float(item[0]) + ( float(item[1]) + float(item[2])/60 )/60
# lat = float(item[3]) + ( float(item[4]) + float(item[5])/60 )/60
# break
if lat == None or lon == None:
return out
out['lat'] = lat
out['lng'] = lon
out['status'] = True
return out

更多 proj 資訊：

$ proj -le
MERIT a=6378137.0 rf=298.257 MERIT 1983
SGS85 a=6378136.0 rf=298.257 Soviet Geodetic System 85
GRS80 a=6378137.0 rf=298.257222101 GRS 1980(IUGG, 1980)
IAU76 a=6378140.0 rf=298.257 IAU 1976
airy a=6377563.396 b=6356256.910 Airy 1830
APL4.9 a=6378137.0. rf=298.25 Appl. Physics. 1965
NWL9D a=6378145.0. rf=298.25 Naval Weapons Lab., 1965
mod_airy a=6377340.189 b=6356034.446 Modified Airy
andrae a=6377104.43 rf=300.0 Andrae 1876 (Den., Iclnd.)
aust_SA a=6378160.0 rf=298.25 Australian Natl & S. Amer. 1969
GRS67 a=6378160.0 rf=298.2471674270 GRS 67(IUGG 1967)
bessel a=6377397.155 rf=299.1528128 Bessel 1841
bess_nam a=6377483.865 rf=299.1528128 Bessel 1841 (Namibia)
clrk66 a=6378206.4 b=6356583.8 Clarke 1866
clrk80 a=6378249.145 rf=293.4663 Clarke 1880 mod.
CPM a=6375738.7 rf=334.29 Comm. des Poids et Mesures 1799
delmbr a=6376428. rf=311.5 Delambre 1810 (Belgium)
engelis a=6378136.05 rf=298.2566 Engelis 1985
evrst30 a=6377276.345 rf=300.8017 Everest 1830
evrst48 a=6377304.063 rf=300.8017 Everest 1948
evrst56 a=6377301.243 rf=300.8017 Everest 1956
evrst69 a=6377295.664 rf=300.8017 Everest 1969
evrstSS a=6377298.556 rf=300.8017 Everest (Sabah & Sarawak)
fschr60 a=6378166. rf=298.3 Fischer (Mercury Datum) 1960
fschr60m a=6378155. rf=298.3 Modified Fischer 1960
fschr68 a=6378150. rf=298.3 Fischer 1968
helmert a=6378200. rf=298.3 Helmert 1906
hough a=6378270.0 rf=297. Hough
intl a=6378388.0 rf=297. International 1909 (Hayford)
krass a=6378245.0 rf=298.3 Krassovsky, 1942
kaula a=6378163. rf=298.24 Kaula 1961
lerch a=6378139. rf=298.257 Lerch 1979
mprts a=6397300. rf=191. Maupertius 1738
new_intl a=6378157.5 b=6356772.2 New International 1967
plessis a=6376523. b=6355863. Plessis 1817 (France)
SEasia a=6378155.0 b=6356773.3205 Southeast Asia
walbeck a=6376896.0 b=6355834.8467 Walbeck
WGS60 a=6378165.0 rf=298.3 WGS 60
WGS66 a=6378145.0 rf=298.25 WGS 66
WGS72 a=6378135.0 rf=298.26 WGS 72
WGS84 a=6378137.0 rf=298.257223563 WGS 84
sphere a=6370997.0 b=6370997.0 Normal Sphere (r=6370997)