網頁

2010年8月11日 星期三

TWD97轉經緯度WGS84 經緯度WGS84轉TWD97

在台灣做地理資訊系統相關工作常常會遇到坐標轉換的問題,而有一件很有趣的事情,在網路上面有很多坐標轉換公式、很多的轉換方法、很多已經寫好的轉換程式,但是就是很難查到如果要寫成程式到底該怎麼做?如果現在必須要在TWD97的坐標系統下進行經緯度定位,我到底該如何在自己的程式中做這件事情?

總之,我們可以查到一堆公式、一堆論文、一堆已經寫好的程式、什麼三參數六參數,這個老師那個老師,然後切換一下VS的視窗,發現查了一個早上也沒看到程式多一行。

今天我也遇到一樣的問題,必須要轉換經緯度成TWD97。

有幸,在網路上查到[Victor的程式設計遇上小提琴]這篇文章,而他所參考的原始公式為一個Steve Dutch學者的文章,在Victor的部落格中已經提供了TWD97與WGS84兩者互轉的Python與JAVA程式,但我所執行的專案大部分為VB或C#語法,所以在這邊貼出改寫他程式碼後的C#語法,希望可以幫到找不到轉換程式寫法的人。

為了方便看到這篇的人複製貼上,以下的程式碼可以直接在VS中於App_Code裡面新增一個CoordinateTransform.cs的class,然後直接全部貼上。

共有兩個Public的函式
1. lonlat_To_twd97輸入XY度分秒或是XY的弧度值返回TWD97 "X,Y"字串。
2. TWD97_To_lonlat輸入TWD97 XY值與type(1表示傳回度分秒字串、2表示傳回弧度字串)。
當然還有兩個重點就是計算轉換的程式碼。
1. Cal_lonlat_To_twd97
2. Cal_TWD97_To_lonlat

所以要使用的時候很簡單,在你的程式裡面先把這個Class new起來,然後傳入你要轉換的值,他會回你一個字串,再自行用字串拆解拆開。
例:

CoordinateTransform CoordinateTransform = new CoordinateTransform();
this.Label1.Text = CoordinateTransform.TWD97_To_lonlat(248170.826, 2652129.977, 2);

Class程式碼:

public class CoordinateTransform
{
double a = 6378137.0;
double b = 6356752.314245;
double lon0 = 121 * Math.PI / 180;
double k0 = 0.9999;
int dx = 250000;

public CoordinateTransform()
{
//
// TODO: 在此加入建構函式的程式碼
//
}

//給WGS84經緯度度分秒轉成TWD97坐標
public string lonlat_To_twd97(int lonD,int lonM,int lonS,int latD,int latM,int latS)
{
double RadianLon = (double)(lonD) + (double)lonM / 60 + (double)lonS / 3600;
double RadianLat = (double)(latD) + (double)latM / 60 + (double)latS / 3600;
return Cal_lonlat_To_twd97(RadianLon, RadianLat);
}
//給WGS84經緯度弧度轉成TWD97坐標
public string lonlat_To_twd97(double RadianLon, double RadianLat)
{
return Cal_lonlat_To_twd97(RadianLon, RadianLat);
}

//給TWD97坐標 轉成 WGS84 度分秒字串 (type1傳度分秒 2傳弧度)
public string TWD97_To_lonlat(double XValue, double YValue, int Type)
{

string lonlat = "";

if (Type == 1)
{
string[] Answer = Cal_TWD97_To_lonlat(XValue, YValue).Split(',');
int LonDValue = (int)double.Parse(Answer[0]);
int LonMValue = (int)((double.Parse(Answer[0]) - LonDValue) * 60);
int LonSValue = (int)((((double.Parse(Answer[0]) - LonDValue) * 60) - LonMValue) * 60);

int LatDValue = (int)double.Parse(Answer[1]);
int LatMValue = (int)((double.Parse(Answer[1]) - LatDValue) * 60);
int LatSValue = (int)((((double.Parse(Answer[1]) - LatDValue) * 60) - LatMValue) * 60);

lonlat = LonDValue + "度" + LonMValue + "分" + LonSValue + "秒," + LatDValue + "度" + LatMValue + "分" + LatSValue + "秒,";
}
else if (Type == 2)
{
lonlat = Cal_TWD97_To_lonlat(XValue, YValue);
}

return lonlat;
}



private string Cal_lonlat_To_twd97(double lon ,double lat)
{
string TWD97 = "";

lon = (lon/180) * Math.PI;
lat = (lat/180) * Math.PI;

//---------------------------------------------------------
double e = Math.Pow((1 - Math.Pow(b,2) / Math.Pow(a,2)), 0.5);
double e2 = Math.Pow(e,2)/(1-Math.Pow(e,2));
double n = ( a - b ) / ( a + b );
double nu = a / Math.Pow((1-(Math.Pow(e,2)) * (Math.Pow(Math.Sin(lat), 2) ) ) , 0.5);
double p = lon - lon0;
double A = a * (1 - n + (5/4) * (Math.Pow(n,2) - Math.Pow(n,3)) + (81/64) * (Math.Pow(n, 4) - Math.Pow(n, 5)));
double B = (3 * a * n/2.0) * (1 - n + (7/8.0)*(Math.Pow(n,2) - Math.Pow(n,3)) + (55/64.0)*(Math.Pow(n,4) - Math.Pow(n,5)));
double C = (15 * a * (Math.Pow(n,2))/16.0)*(1 - n + (3/4.0)*(Math.Pow(n,2) - Math.Pow(n,3)));
double D = (35 * a * (Math.Pow(n,3))/48.0)*(1 - n + (11/16.0)*(Math.Pow(n,2) - Math.Pow(n,3)));
double E = (315 * a * (Math.Pow(n,4))/51.0)*(1 - n);

double S = A * lat - B * Math.Sin(2 * lat) +C * Math.Sin(4 * lat) - D * Math.Sin(6 * lat) + E * Math.Sin(8 * lat);

//計算Y值
double K1 = S*k0;
double K2 = k0*nu*Math.Sin(2*lat)/4.0;
double K3 = (k0*nu*Math.Sin(lat)*(Math.Pow(Math.Cos(lat),3))/24.0) * (5 - Math.Pow(Math.Tan(lat),2) + 9*e2*Math.Pow((Math.Cos(lat)),2) + 4*(Math.Pow(e2,2))*(Math.Pow(Math.Cos(lat),4)));
double y = K1 + K2*(Math.Pow(p,2)) + K3*(Math.Pow(p,4));

//計算X值
double K4 = k0*nu*Math.Cos(lat);
double K5 = (k0*nu*(Math.Pow(Math.Cos(lat),3))/6.0) * (1 - Math.Pow(Math.Tan(lat),2) + e2*(Math.Pow(Math.Cos(lat),2)));
double x = K4 * p + K5 * (Math.Pow(p, 3)) + dx;

TWD97 = x.ToString()+ "," + y.ToString();
return TWD97;
}


private string Cal_TWD97_To_lonlat(double x, double y)
{

double dy = 0;
double e = Math.Pow((1- Math.Pow(b,2)/Math.Pow(a,2)), 0.5);

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);

// 計算緯度
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);

// 計算經度
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);

lat = (lat * 180) / Math.PI; //緯
lon = (lon * 180) / Math.PI; //經


string lonlat = lon + "," + lat;
return lonlat;
}


}

最後,我用控制點去進行轉換也覺得誤差很小,如果要進行一般性的應用,我想是相當足夠。

題外話:剛剛又去Victor的部落格看了一看,我猜測他是資工的,沒想到這樣平凡、常用、實用的內容,竟然不是在所謂的測量或地理資訊相關網站找到,我真的不知道該說什麼。

13 則留言:

愛迪生。吳 提到...

我也碰到了這個問題
說真的,我不覺得這些東西很困難,只是對我這種沒有相關基礎知識的人來說,要先k完知識再寫那麼一隻程式出來實在很沒有意思
anyway, 謝謝你的分享 :)
你的程式碼很完整,比victor的更容易用了 :)

ola的家 提到...

坐標轉換就學理部份是一個很高深的學問啦,說真的也不是基礎知識讀一讀就可以用出來的,很多參數會算會帶也不見得了解其中的涵義,但是就應用面來說,一個精度尚可的轉換方式就足以應付絕大部分的場合了。

匿名 提到...

謝謝你的提供~
面對艱深的座標轉換的英文資料實在令人無奈...
在此看到你的無私真的太高興了!!

繆拉 提到...

板主您好

感謝您釋出原始碼,我最近也有相同的需求,在參考了您的程式之後,對於轉換函數本身做了一些修正,並放在我的blog上:

http://sask989.blogspot.com/2012/05/twd97tm2.html

經過測試後精度有再提升,在此供您以及大家參考 :)

ola的家 提到...

Oh ..... really great

@o@

Unknown 提到...

老師你好:
在private string Cal_lonlat_To_twd97(double lon ,double lat)這一段程式碼中
有一段
//計算X值
double K4 = k0*nu*Math.Cos(lat);
double K5 = (k0*nu*(Math.Pow(Math.Cos(lat),3))/6.0) * (1 - Math.Pow(Math.Tan(lat),2) + e2*(Math.Pow(Math.Cos(lat),2)));
請問為什麼計算x值是用lat? 因為前一段計算Y值已使用lat, 而計算x值是不是應該用lon?

Unknown 提到...

您好:
在您的程式撰寫 Cal_TWD97_To_lonlat中C1參數為 double C1 = Math.Pow(e2*Math.Cos(fp), 2);
而在Converting UTM to Latitude and Longitude (Or Vice Versa)網頁中C1參數為C1 = e'2cos2(fp)
不知道是哪一個才是正確的呢

ola的家 提到...

一個是算式,一個是程式碼

Unknown 提到...

有沒有辦法再同場加映
其他程式語言
javascript VBA
找了很久只有你的最接近我們入門款看得懂的

法拉魚 提到...

由衷感謝這篇,不然會吐血..

Kyle K. Lin 提到...

非常棒,轉換成 Java 也沒問題,感謝。

天亮damody 提到...

感謝 我也順便做了一個rust的版本
https://crates.io/crates/mercator

Unknown 提到...

超讚的
這是我看過,反算誤差最小的
謝謝大大無私地分享

張貼留言