大胖子與小個子的部落格: WGS84經緯度與TWD97(TM2)投影坐標轉換程式

最近處理 GIS 資料時，遇到了 WGS84 經緯度與 TWD97 投影坐標互相轉換的需求，於是在網路上找到了程式設計遇上小提琴以及 ola的家兩個網頁，分別提供了 Python 與 C# 的原始碼。不過實際測試後卻發現經緯度轉 TWD97 有 2 公尺， TWD97 轉經緯度有 3 公厘左右的誤差值。這對精密測量來說，還是稍嫌不足。

因此，我又尋找了其他的資料，最後在中國的代码发芽网找到了經緯度與 UTM 互轉的程式碼，由於 UTM 以及 TWD97 都是同一種投影方式（橫麥卡托投影，Transverse Mercator Projection），只有 X 軸平移量、中央子午線以及投影尺度因子等常數的不同而已（尤瑞哲，2003）。因此我參考了 ola的家所提供的程式框架，並將經緯度轉 TWD97 的函式重新改寫成 USGS (U.S. Geological Survey) 提供的轉換公式； TWD97 轉經緯度的轉換函式則是修正了其中一個算式。

經過測試，將 TWD97 投影坐標轉成 WGS84 經緯度後再轉回 TWD97 ，其轉換前後的誤差值在 0.1 公厘以下，大幅提升轉換的正確性，在此將 C# 原始碼提供出來給有需要的人參考。

使用方法和 ola的家提供的程式相同。

經緯度轉 TWD97 ：

 CoordinateTransform CoordinateTransform = new CoordinateTransform();

string twd97 = CoordinateTransform.lonlat_To_twd97(121.5, 23.5);

TWD97_E.Text = twd97.Substring(0, twd97.IndexOf(","));

TWD97_N.Text = twd97.Substring(twd97.IndexOf(",") + 1, twd97.Length - (twd97.IndexOf(",") + 1));

TWD97 轉經緯度：

 CoordinateTransform CoordinateTransform_Test = new CoordinateTransform();

string lonlat = CoordinateTransform_Test.TWD97_To_lonlat(System.Convert.ToDouble(TWD97_E.Text), System.Convert.ToDouble(TWD97_N.Text), 2);

Long.Text = lonlat.Substring(0, lonlat.IndexOf(","));

Lat.Text = lonlat.Substring(lonlat.IndexOf(",") + 1, lonlat.Length - (lonlat.IndexOf(",") + 1));

原始碼：

public class CoordinateTransform

{

  private static double a = 6378137.0;

  private static double b = 6356752.3142451;

  private double lon0 = 121 * Math.PI / 180;

  private double k0 = 0.9999;

  private int dx = 250000;

  private int dy = 0;

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

  private double e2 = (1 - Math.Pow(b, 2) / Math.Pow(a, 2)) / (Math.Pow(b, 2) / Math.Pow(a, 2));



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

      double LonSValue = (((double.Parse(Answer[0]) - LonDValue) * 60) - LonMValue) * 60;



      int LatDValue = (int)double.Parse(Answer[1]);

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

      double LatSValue = (((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 - Math.Floor((lon + 180) / 360) * 360) * Math.PI / 180;

    lat = lat * Math.PI / 180;



    double V = a / Math.Sqrt(1 - e * Math.Pow(Math.Sin(lat), 2));

    double T = Math.Pow(Math.Tan(lat), 2);

    double C = e2 * Math.Pow(Math.Cos(lat), 2);

    double A = Math.Cos(lat) * (lon - lon0);

    double M = a *((1.0 - e / 4.0 - 3.0 * Math.Pow(e, 2) / 64.0 - 5.0 * Math.Pow(e, 3) / 256.0) * lat -

          (3.0 * e / 8.0 + 3.0 * Math.Pow(e, 2) / 32.0 + 45.0 * Math.Pow(e, 3) / 1024.0) *

          Math.Sin(2.0 * lat) + (15.0 * Math.Pow(e, 2) / 256.0 + 45.0 * Math.Pow(e, 3) / 1024.0) * 

          Math.Sin(4.0 * lat) - (35.0 * Math.Pow(e, 3) / 3072.0) * Math.Sin(6.0 * lat));

    // x

    double x = dx + k0 * V * (A + (1 - T + C) * Math.Pow(A, 3) / 6 + (5 - 18 * T + Math.Pow(T, 2) + 72 * C - 58 * e2) * Math.Pow(A, 5) / 120);

    // y

    double y = dy + k0 * (M + V * Math.Tan(lat) * (Math.Pow(A, 2) / 2 + (5 - T + 9 * C + 4 * Math.Pow(C, 2)) * Math.Pow(A, 4) / 24 + ( 61 - 58 * T + Math.Pow(T, 2) + 600 * C - 330 * e2) * Math.Pow(A, 6) / 720));



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

    return TWD97;

  }



  private string Cal_TWD97_To_lonlat(double x, double y)

  {

    x -= dx;

    y -= dy;



    // Calculate the Meridional Arc

    double M = y / k0;



    // Calculate Footprint Latitude

    double mu = M / (a * (1.0 - e / 4.0 - 3 * Math.Pow(e, 2) / 64.0 - 5 * Math.Pow(e, 3) / 256.0));

    double e1 = (1.0 - Math.Sqrt(1.0 - e)) / (1.0 + Math.Sqrt(1.0 - e));



    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 C1 = e2 * Math.Pow(Math.Cos(fp), 2);

    double T1 = Math.Pow(Math.Tan(fp), 2);

    double R1 = a * (1 - e) / Math.Pow((1 - e * Math.Pow(Math.Sin(fp), 2)), (3.0 / 2.0));

    double N1 = a / Math.Pow((1 - e * 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;

  }

}

原始碼下載： SkyDrive

參考資料：
尤瑞哲，2003。測量坐標系統第二版，國立成功大學，台南，第64頁。

大胖子與小個子的部落格

2012年5月9日星期三

WGS84經緯度與TWD97(TM2)投影坐標轉換程式

2 則留言:

2012年5月9日 星期三

WGS84經緯度與TWD97(TM2)投影坐標轉換程式

2 則留言:

2012年5月9日星期三