2012年5月9日 星期三

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頁。

1 則留言: