經緯度轉換TWD97

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

麥卡脫投影

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

地圖的變形

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

projection

如上圖所示,我們可以發現,我們用紙圓筒將地球包起來,可是接觸地球表面的點只有一個,而其它地方離地球表面越來越遠,這造成從地心投出的光線越偏越多,也就是說,離這條線越遠的地方,地圖的變形越大,所以其實每張地圖並不是所有地方的大小尺吋都是正確的,只有那條線是正確的,往兩邊移變形會越來越大,因此,如果我們以日本為此線的中央的話,來到台灣地圖可能已經嚴重變形,我們都知道土地吋土吋金,想想看如果我們的地圖存在著很嚴重的變形問題,可能原本小小的圖地在地圖上變形成很大的一塊,這樣是沒辦法接受的,為此,每個國家通常都會有他們自己的投影參數,以適合他們自己的國家,不會讓地圖變形太大,當灣當然也有自己的地圖,因此也需要自己的投影參數基準,而台灣自己的投影座標基準就叫做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)(n2 - n3) + (55/64)(n4 - n5) ...]

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

E’ = (315 tan4/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網友,我一個外行人寫這非我專業領域的東西光是背景知識就弄得一頭霧水時,如果沒有他的熱心幫忙,我可能一輩子也寫不出來

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

This entry was posted in 中文文章, 分享 and tagged , , , , , , , . Bookmark the permalink.

14 Responses to 經緯度轉換TWD97

  1. Summer says:

    您好:小弟近日也在試著寫一寫有關LBS的東西,剛好也是需要用到座標轉換,很高興看到您已經將公式實作出來,不知道是否也可以分享TMParameter這一支CLASS的內容一併讓小弟參考呢?我是想要做TWD67 TM2 -> TWD97 經緯度的轉換。

  2. Spring says:

    您好,個人對TMParameter這個class深感興趣,
    與您請教TMParameter程式中的變數,謝謝。

    1. 請問 double e = tm.getE();
    請問e的值為何?

    2. 其的變數的值是否如下所示?

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

    以上二個問題,煩請您解惑,感恩。

    • Joe says:

      您好

      我也有跟Spring一樣的問題,有幾個變數不知道確切的值,不知道您方便回覆給我嗎?

      以下是不確定的值 :
      dx =
      dy =
      lon0 = 121
      k0 = 0.9999
      a = 6378137.0公尺
      b = 6356752.314245公尺
      e =

      感謝您~

  3. dego says:

    您好,小弟最近也是在測,可是改您的經緯轉97
    我是用VB,有些運算元好像順訊跟優先計算次序不同,所以會有誤差的數值出來,可以請問一下MAIL原始的數學運算方程是參考嗎,謝謝!!!

  4. Skai says:

    您好,這邊的 k0=0.9999
    與網站裡的0.9996不同
    請問是筆誤還是另有含意呢?

  5. victor says:

    @Skai

    k0是一個參數,可依照你想要的比例來條整經緯度的比例,不同的大地基準有不同的參數,TWD97給的k0就是0.99999

    我的認知是這樣,不過我不是這方面專家,關於k0的正確定義我也不清楚

  6. Skai says:

    後來查到是投影方式的差異造成的變化,感謝!

  7. Skai says:

    貼出去後才發現你已經回我了,謝謝你喔。

  8. compileerrors says:

    不好意思我想請問一下……..
    您JAVA程式所表示的TMParameter
    tm.getDx()…..的dx
    應該跟Python裡定義的一樣?
    我不了解的是JAVA裡有tm.getDy()
    Python確沒有……
    是否可以解惑一下?
    謝謝……

  9. sky says:

    下載了您的電力座標定位app
    但目前在戶外行程中
    更需要wgs84轉換成twd67的app
    能越多轉換格式當然更好

    不知道版大是否會寫相關app
    而造福廣大的戶外玩家呢
    或將電力座標的app修改成其他欄位也可輸入
    而算出其他座標呢

    非常感謝

  10. mercator says:

    投影程式可以參考 Proj.4
    足夠你用了

  11. mercator says:

    版主文中有提到ko這個參數,恐怕很多人搞不清楚是甚麼東西。
    ko的全名是中央經線尺度比。
    在TM投影中常用的帶寬是2度,3度,6度。
    因為不同的帶寬就會產生不同的變形量(因為把圓弧伸展到平面上必然的會產生變形)。
    中央經線就橫切面來看,就是那一條圓弧和平面接觸的點。
    沿著中央經線,尺度比是1.0,越向兩側變形量就越大。為了讓整幅圖可以達到一個比較平均的變形量,就加上一個尺度比。讓圓弧和投影面從切線變成割線。
    尺度比一般規定
    2度 0.9999
    3度 1.0000
    6度 0.9996
    尺度比通常出現在美系地圖。
    中國大陸叫高斯克呂格投影,是不用尺度比的。

  12. Pingback: 座標轉換 ← Tom's Dev Blog

  13. Cedric says:

    您好, 由於援引了您的程式創作了 GPS 羅盤及座標轉換的 Android app, 特此知會: https://play.google.com/store/apps/details?id=com.ubiris.twdcompass

    雖然此 app 純以山難預防及救助角度開發, 完全免費也不含任何廣告, 但若您認為援引您的創作有不妥之處, 或不願在 app 說明頁被提及來源出處, 還請不吝告知.

    也非常感謝您願意分享此智慧財產!

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>