經緯度轉換TWD97

最近因為專題上的需要，需要把經緯度轉成TWD97的座標，而在更早之前我寫的電力座標轉換程式，則是需要把TWD97的座標轉換成經緯度，今天我們來談談關於這些地理座標的轉換

麥卡脫投影

首先，我們都知道，地球是圓型的，更精確的來說，是橢圓型的，更正確來說，是不完美的橢圓型，而我們常用的地圖是平面的，為了能把地圖顯示在平面上，我們需要把地圖投影到平面上，而投影的方式有很多種，我們常見的地圖通常是使用麥卡脫投影，想像一下，用一張紙，捲起來成圓筒將地球包覆，在地球的中心擺一個光源，而從這光源發射的光線，透過地球表面映在紙面上，我們把這些映出來的影像畫下來，接著把紙攤平開來，如此一來我們就有了平面的地圖，這個過程就稱為麥卡脫投影，其實就只是把曲面的幾合圖型對應到平面的過程而已

地圖的變形

而經緯度是將地球以角度的方式分割的一種座標系，那什麼又是TWD97? 如我們剛才所說的，經緯度是以角度的方式來計算的座標，如果我們有一張地圖，光有經緯度我們沒辦法指出這個經緯度是在地圖上的哪裡，這時候我們就要經過麥卡脫投影，將經緯度透過某些特定的參數投影到地圖上，換算成平面的座標，單位為公尺，如此一來我們就可以較準確地指出，這個經緯度對應到地圖上的座標是在哪裡，而投影是會變形的，因為我們把曲面的東西投影在平面上，多少都會有變形

如上圖所示，我們可以發現，我們用紙圓筒將地球包起來，可是接觸地球表面的點只有一個，而其它地方離地球表面越來越遠，這造成從地心投出的光線越偏越多，也就是說，離這條線越遠的地方，地圖的變形越大，所以其實每張地圖並不是所有地方的大小尺吋都是正確的，只有那條線是正確的，往兩邊移變形會越來越大，因此，如果我們以日本為此線的中央的話，來到台灣地圖可能已經嚴重變形，我們都知道土地吋土吋金，想想看如果我們的地圖存在著很嚴重的變形問題，可能原本小小的圖地在地圖上變形成很大的一塊，這樣是沒辦法接受的，為此，每個國家通常都會有他們自己的投影參數，以適合他們自己的國家，不會讓地圖變形太大，當灣當然也有自己的地圖，因此也需要自己的投影參數基準，而台灣自己的投影座標基準就叫做TWD97(Taiwan Datum 1997)，其實在這之前有更早的TWD67，也只是較早期在用的投影基準，因為兩者投影參數的不同，所以即使同樣的座標，在實際上的位置其實是不一樣的，我們在此只介紹TWD97

TWD97

那TWD97又是怎麼樣的一個大地基準呢? 他是一個以東經121度為中央經線的大地基準，也就是，想像一下拿一張紙，將地球包起來，而接觸到地球的那條線經過東經121度，TWD97使用的就是這樣的投影方式，而我們先前有提到地球是不完美的橢圓，為了要計算投影的座標，TWD97一樣有一些參數是地球的橢圓半徑用來計算投影，因為不同的參數算出來的數字會不同，因此每個大地座標基準的參數一樣要規定

公式

現在我們知道了TWD97是從地球投影到地圖的座標系的基準，可是光知道這些其實沒什麼用處，重點是我們要能算出來，這時我們就需要投影的公式，這投影公式可以參考這個網頁

Converting UTM to Latitude and Longitude (Or Vice Versa)

從這網頁裡我們可以發現這個公式實在是又臭又長，我在打成程式時就出了不少錯，看到眼睛都快花掉，首先先看他的參數，lat是我們輸入的緯度，long是我們輸入的經度，接著就是關於投影的一些參數，a是地球赤道半徑(Equatorial Radius)，而b是兩極半徑(Polar Radius)，而long0是指中央經線，也就是我所說的卷紙接觸地球的線，而k0是延著long0的縮放比例，長度單位都是公尺，而角度都是用弧度，我之前一直算錯就有因為是用成角度來計算，全部都要是弧度才對，我們在此用的參數是

a = 6378137.0公尺
b = 6356752.314245公尺
long0 = 121度
k0 = 0.9999

而這網頁上打的公式有一些錯誤，像是

B’ = (3 tan/2)[1 – n + (7/8)(n² – n³) + (55/64)(n⁴ – n⁵) …]

一開始我弄得一頭霧水，tan/2到底是什麼樣的標記方式，一開始我猜是tan(1/2)，但是錯了，後來猜是tan代入後面那串除以二，也是錯，又或著代入後面那串後再除以二，還是錯，找了其它資料才發現是他打錯了，那t是多餘的，所以3tan/2其實是3an/2，光為了這點浪費了我不少時間，另一個錯誤是

E’ = (315 tan⁴/512)[1 – n  …]

512應該是51才對

計算程式

有了參數和公式，我們就可以寫程式來幫我們算，因為這用手算實在是太累人了，光看公式就眼花了，更何況是手算，計算是電腦的強項，當然是閃開，讓專業的來，以下是Lat/Lon轉TWD97的Python原始碼

from math import tan, sin, cos, radians

class LatLonToTWD97(object):
    """This object provide method for converting lat/lon coordinate to TWD97
    coordinate

    the formula reference to
    http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm (there is lots of typo)
    http://www.offshorediver.com/software/utm/Converting UTM to Latitude and Longitude.doc

    Parameters reference to
    http://rskl.geog.ntu.edu.tw/team/gis/doc/ArcGIS/WGS84%20and%20TM2.htm
    http://blog.minstrel.idv.tw/2004/06/taiwan-datum-parameter.html
    """

    def __init__(self,
        a = 6378137.0,
        b = 6356752.314245,
        long0 = radians(121),
        k0 = 0.9999,
        dx = 250000,
    ):
        # Equatorial radius
        self.a = a
        # Polar radius
        self.b = b
        # central meridian of zone
        self.long0 = long0
        # scale along long0
        self.k0 = k0
        # delta x in meter
        self.dx = dx

    def convert(self, lat, lon):
        """Convert lat lon to twd97

        """
        a = self.a
        b = self.b
        long0 = self.long0
        k0 = self.k0
        dx = self.dx

        e = (1-b**2/a**2)**0.5
        e2 = e**2/(1-e**2)
        n = (a-b)/(a+b)
        nu = a/(1-(e**2)*(sin(lat)**2))**0.5
        p = lon-long0

        A = a*(1 - n + (5/4.0)*(n**2 - n**3) + (81/64.0)*(n**4  - n**5))
        B = (3*a*n/2.0)*(1 - n + (7/8.0)*(n**2 - n**3) + (55/64.0)*(n**4 - n**5))
        C = (15*a*(n**2)/16.0)*(1 - n + (3/4.0)*(n**2 - n**3))
        D = (35*a*(n**3)/48.0)*(1 - n + (11/16.0)*(n**2 - n**3))
        E = (315*a*(n**4)/51.0)*(1 - n)

        S = A*lat - B*sin(2*lat) + C*sin(4*lat) - D*sin(6*lat) + E*sin(8*lat)

        K1 = S*k0
        K2 = k0*nu*sin(2*lat)/4.0
        K3 = (k0*nu*sin(lat)*(cos(lat)**3)/24.0) * 
            (5 - tan(lat)**2 + 9*e2*(cos(lat)**2) + 4*(e2**2)*(cos(lat)**4))

        y = K1 + K2*(p**2) + K3*(p**4)

        K4 = k0*nu*cos(lat)
        K5 = (k0*nu*(cos(lat)**3)/6.0) * 
            (1 - tan(lat)**2 + e2*(cos(lat)**2))

        x = K4*p + K5*(p**3) + self.dx
        return x, y
if __name__ == '__main__':
    from math import degrees

    c = LatLonToTWD97()
    lat = radians(float(raw_input('lat:')))
    lon = radians(float(raw_input('lon:')))
    print 'input lat/lon', degrees(lat), degrees(lon)
    x, y = c.convert(lat, lon)
    print (x, y)

經過我找網路上的三角點資料來計算後，跟那些資料比起來的誤差都在公分以下，不過這個程式只做為參考用，不負任何責任

TWD97轉lat/lon

由於我先前寫的電力座標定位程式，台電的電力座標用的是地圖座標，因此我需要將TWD97轉換成Lat/Lon才能在google地圖上定位，為此我也寫了TWD97轉成lat/lon的JAVA程式，

package com.ez2learn.android.powergrid.geo;

public class TMToLatLon {
	/**
		@brief TM coordinate to lat lon formula

		reference to http://www.uwgb.edu/dutchs/UsefulData/UTMFormulas.htm

		@param TMParameter object
		@param x
		@param y
	**/
	static public double[] convert(TMParameter tm, double x, double y) {
		double dx = tm.getDx();
		double dy = tm.getDy();
		double lon0 = tm.getLon0();
		double k0 = tm.getK0();
		double a = tm.getA();
		double b = tm.getB();
		double e = tm.getE();

		x -= dx;
		y -= dy;

		// Calculate the Meridional Arc
		double M = y/k0;

		// Calculate Footprint Latitude
		double mu = M/(a*(1.0 - Math.pow(e, 2)/4.0 - 3*Math.pow(e, 4)/64.0 - 5*Math.pow(e, 6)/256.0));
		double e1 = (1.0 - Math.pow((1.0 - Math.pow(e, 2)), 0.5)) / (1.0 + Math.pow((1.0 - Math.pow(e, 2)), 0.5));

		double J1 = (3*e1/2 - 27*Math.pow(e1, 3)/32.0);
		double J2 = (21*Math.pow(e1, 2)/16 - 55*Math.pow(e1, 4)/32.0);
		double J3 = (151*Math.pow(e1, 3)/96.0);
		double J4 = (1097*Math.pow(e1, 4)/512.0);

		double fp = mu + J1*Math.sin(2*mu) + J2*Math.sin(4*mu) + J3*Math.sin(6*mu) + J4*Math.sin(8*mu);

		// Calculate Latitude and Longitude

		double e2 = Math.pow((e*a/b), 2);
		double C1 = Math.pow(e2*Math.cos(fp), 2);
		double T1 = Math.pow(Math.tan(fp), 2);
		double R1 = a*(1-Math.pow(e, 2))/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), (3.0/2.0));
		double N1 = a/Math.pow((1-Math.pow(e, 2)*Math.pow(Math.sin(fp), 2)), 0.5);

		double D = x/(N1*k0);

		// lat
		double Q1 = N1*Math.tan(fp)/R1;
		double Q2 = (Math.pow(D, 2)/2.0);
		double Q3 = (5 + 3*T1 + 10*C1 - 4*Math.pow(C1, 2) - 9*e2)*Math.pow(D, 4)/24.0;
		double Q4 = (61 + 90*T1 + 298*C1 + 45*Math.pow(T1, 2) - 3*Math.pow(C1, 2) - 252*e2)*Math.pow(D, 6)/720.0;
		double lat = fp - Q1*(Q2 - Q3 + Q4);

		// long
		double Q5 = D;
		double Q6 = (1 + 2*T1 + C1)*Math.pow(D, 3)/6;
		double Q7 = (5 - 2*C1 + 28*T1 - 3*Math.pow(C1, 2) + 8*e2 + 24*Math.pow(T1, 2))*Math.pow(D, 5)/120.0;
		double lon = lon0 + (Q5 - Q6 + Q7)/Math.cos(fp);

		return new double[] {Math.toDegrees(lat), Math.toDegrees(lon)};
	}
}

先前也有寫一個python版本的做測試，不過弄丟了，雖然後來我再寫了一份，不過沒有檢查過，所以在此就不提供，此程式一樣僅供參考，不負任何責任

特別感謝

特別感謝ptt上的nehex網友，我一個外行人寫這非我專業領域的東西光是背景知識就弄得一頭霧水時，如果沒有他的熱心幫忙，我可能一輩子也寫不出來

最後，因為這不是我的專業，如果上面的說法有任何錯誤，麻煩請指正一下，謝謝