Taiwan Datum (TWD67, TWD97, TM2 Source)

台灣座標系統轉換原始碼 Taiwan Datum Convert Source(TWD67, TWD97, TM2 Source)


Related Site:

洪朝貴的首頁 (Chinese)

積丹尼 = Dan Jacobson (Chinese/English)

Dan Jacobson (English)

Source Below:




/*
TWD67/TWD97/TM2 datum convert functions

based on: GFC (
http://www.samblackburn.com/gfc/)
proj (
http://www.remotesensing.org/proj)

rewrite: minstrel(http://www.minstrel.idv.tw), 2004/6/1

license:
GNU General Public License
*/

#include <math.h>

// from long-lat to tm2
void toTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y);

// from tm2 to long-lat
void fromTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y);

// used by toTM2 and fromTM2
double mercator(double y, double a, double ecc);

// from TWD97 to TWD67
void toTWD67(double* x, double* y, double* z);

// from TWD67 to TWD97
void toTWD97(double* x, double* y, double* z);

/*
Definition of math related value
*/

#define COS67_5 0.3826834323650897717284599840304e0
#define PI 3.14159265358979323e0
#define HALF_PI 1.570796326794896615e0
#define DEG_RAD 0.01745329251994329572e0
#define RAD_DEG 57.295779513082321031e0

/*
Definition of datum related value
*/

#define AD_C 1.0026000e0

#define TWD67_A 6378160.0e0
#define TWD67_B 6356774.7192e0
#define TWD67_ECC 0.00669454185458e0
#define TWD67_ECC2 0.00673966079586e0
#define TWD67_DX -752.32e0 // different from garmin and already knowned value, but those all value only
#define TWD67_DY -361.32e0 // got 5-15m accuracy. the real offical value is holded by somebody and not
#define TWD67_DZ -180.51e0 // release to public. if can got more enough twd67/twd97 control point coordinare,
#define TWD67_RX -0.00000117e0 // then we can calculate a better value than now.
#define TWD67_RY 0.00000184e0 //
#define TWD67_RZ 0.00000098e0 // and, also lack twd67/twd97 altitude convertion value...
#define TWD67_S 0.00002329e0 //


#define TWD97_A 6378137.0e0
#define TWD97_B 6356752.3141e0
#define TWD97_ECC 0.00669438002290e0
#define TWD97_ECC2 0.00673949677556e0

#define TWD67_TM2 0.9999e0 // TWD67->TM2 scale
#define TWD97_TM2 0.9999e0 // TWD97->TM2 scale

/*
datum convert function
*/

void toTWD97(double* x, double* y, double* z)
{
double r, pole, sin_lat, cos_lat;
double lat, lon, height;
double x1, y1, z1, x2, y2, z2;
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

lon = *x*DEG_RAD;
lat = *y*DEG_RAD;
height = *z*DEG_RAD;

if(lat<-HALF_PI&&lat>-1.001*HALF_PI)
lat = -HALF_PI;
else if(lat>HALF_PI&&lat<1.001*HALF_PI)
lat = HALF_PI;
else if((lat<-HALF_PI)||(lat>HALF_PI))
return;

if(lon>PI)
lon -= (2*PI);

sin_lat = sin(lat);
cos_lat = cos(lat);
r = TWD67_A/(sqrt(1.0-TWD67_ECC*sin_lat*sin_lat));
x1 = (r+height)*cos_lat*cos(lon);
y1 = (r+height)*cos_lat*sin(lon);
z1 = ((r*(1-TWD67_ECC))+height)*sin_lat;

x2 = x1+TWD67_DX+TWD67_S*(lon+TWD67_RZ*lat-TWD67_RY*height);
y2 = y1+TWD67_DY+TWD67_S*(-TWD67_RZ*lon+lat+TWD67_RX*height);
z2 = z1+TWD67_DZ+TWD67_S*(TWD67_RY*lon-TWD67_RX*lat+height);

pole = 0;

if(x2!=0.0)
{
lon=atan2(y2,x2);
}
else
{
if(y2>0)
{
lon = HALF_PI;
}
else if(y2<0)
{
lon=-HALF_PI;
}
else
{
pole=1;

lon=0;

if(z2>0)
{
lat=HALF_PI;
}
else if(z2<0)
{
lat=-HALF_PI;
}
else
{
lat=HALF_PI;

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = -TWD97_B;

return;
}
}
}

q2 = x2*x2+y2*y2;
q = sqrt(q2);
t = z2*AD_C;
s = sqrt(t*t+q2);
sin_b = t/s;
cos_b = q/s;
t1 = z2+TWD97_B*TWD97_ECC2*sin_b*sin_b*sin_b;
sum = q-TWD97_A*TWD97_ECC*cos_b*cos_b*cos_b;
s1 = sqrt(t1*t1+sum*sum);
sin_p = t1/s1;
cos_p = sum/s1;
r = TWD97_A/sqrt(1.0-TWD97_ECC*sin_p*sin_p);

if(cos_p>=COS67_5)
{
height=q/cos_p-r;
}
else if(cos_p<=-COS67_5)
{
height=q/-cos_p-r;
}
else
{
height=z2/sin_p+r*(TWD97_ECC-1.0);
}

if(!pole)
{
lat = atan(sin_p/cos_p);
}

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = height;
}

void toTWD67(double* x, double* y, double* z)
{
double r, pole, sin_lat, cos_lat;
double lat, lon, height;
double x1, y1, z1, x2, y2, z2;
double q, q2, t, t1, s, s1, sum, sin_b, cos_b, sin_p, cos_p;

lon = *x*DEG_RAD;
lat = *y*DEG_RAD;
height = *z*DEG_RAD;

if(lat<-HALF_PI&&lat>-1.001*HALF_PI)
lat = -HALF_PI;
else if(lat>HALF_PI&&lat<1.001*HALF_PI)
lat = HALF_PI;
else if((lat<-HALF_PI)||(lat>HALF_PI))
return;

if(lon>PI)
lon -= (2*PI);

sin_lat = sin(lat);
cos_lat = cos(lat);
r = TWD97_A/(sqrt(1.0-TWD97_ECC*sin_lat*sin_lat));
x1 = (r+height)*cos_lat*cos(lon);
y1 = (r+height)*cos_lat*sin(lon);
z1 = ((r*(1-TWD97_ECC))+height)*sin_lat;

x2 = x1-TWD67_DX-TWD67_S*(lon+TWD67_RZ*lat-TWD67_RY*height);
y2 = y1-TWD67_DY-TWD67_S*(-TWD67_RZ*lon+lat+TWD67_RX*height);
z2 = z1-TWD67_DZ-TWD67_S*(TWD67_RY*lon-TWD67_RX*lat+height);

pole = 0;

if(x2!=0.0)
{
lon=atan2(y2,x2);
}
else
{
if(y2>0)
{
lon = HALF_PI;
}
else if(y2<0)
{
lon=-HALF_PI;
}
else
{
pole=1;

lon=0;

if(z2>0)
{
lat=HALF_PI;
}
else if(z2<0)
{
lat=-HALF_PI;
}
else
{
lat=HALF_PI;

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = -TWD67_B;

return;
}
}
}

q2 = x2*x2+y2*y2;
q = sqrt(q2);
t = z2*AD_C;
s = sqrt(t*t+q2);
sin_b = t/s;
cos_b = q/s;
t1 = z2+TWD67_B*TWD67_ECC2*sin_b*sin_b*sin_b;
sum = q-TWD67_A*TWD67_ECC*cos_b*cos_b*cos_b;
s1 = sqrt(t1*t1+sum*sum);
sin_p = t1/s1;
cos_p = sum/s1;
r = TWD67_A/sqrt(1.0-TWD67_ECC*sin_p*sin_p);

if(cos_p>=COS67_5)
{
height=q/cos_p-r;
}
else if(cos_p<=-COS67_5)
{
height=q/-cos_p-r;
}
else
{
height=z2/sin_p+r*(TWD67_ECC-1.0);
}

if(!pole)
{
lat = atan(sin_p/cos_p);
}

*x = lon*RAD_DEG;
*y = lat*RAD_DEG;
*z = height;
}

void toTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y)
{
double x0, y0, x1, y1, m0, m1;
double n, t, c, A;

x0 = *x*DEG_RAD;
y0 = *y*DEG_RAD;

x1 = lon*DEG_RAD;
y1 = lat*DEG_RAD;

m0 = mercator(y1, a, ecc);
m1 = mercator(y0, a, ecc);

n = a/sqrt(1-ecc*pow(sin(y0), 2.0));
t = pow(tan(y0), 2.0);
c = ecc2*pow(cos(y0), 2.0);
A = (x0-x1)*cos(y0);

*x = scale*n*(A + (1.0 - t + c)*A*A*A/6.0 + (5.0 - 18.0*t + t*t + 72.0*c - 58.0*ecc2)*pow(A, 5.0)/120.0);
*y = scale*(m1 - m0 + n*tan(y0)*(A*A/2.0 + (5.0 - t + 9.0*c + 4*c*c)*pow(A, 4.0)/24.0
+ (61.0 - 58.0*t + t*t + 600.0*c - 330.0*ecc2)*pow(A, 6.0)/720.0));
}

void fromTM2(double a, double ecc, double ecc2, double lat, double lon, double scale, double* x, double* y)
{
double x0, y0, x1, y1, phi, m, m0, mu, e1;
double c1, t1, n1, r1, d;

x0 = *x;
y0 = *y;

x1 = lon*DEG_RAD;
y1 = lat*DEG_RAD;

m0 = mercator(y1, a, ecc);
m = m0 + y0/scale;

e1 = (1.0-sqrt(1.0-ecc))/(1.0+sqrt(1.0-ecc));
mu = m / (a * (1.0 - ecc/4.0 - 3.0*ecc*ecc/64.0 - 5.0*ecc*ecc*ecc/256.0));

phi = mu + (3.0*e1/2.0 - 27.0*pow(e1, 3.0)/32.0)*sin(2.0*mu)
+ (21.0*e1*e1/16.0 - 55.0*pow(e1,4.0)/32.0)*sin(4.0*mu)
+ 151.0*pow(e1,3.0)/96.0*sin(6.0*mu) + 1097.0*pow(e1,4.0)/512.0*sin(8.0*mu);

c1 = ecc2*pow(cos(phi),2.0);
t1 = pow(tan(phi),2.0);
n1 = a/sqrt(1-ecc*pow(sin(phi),2.0));
r1 = a*(1.0-ecc)/pow(1.0-ecc*pow(sin(phi),2.0), 1.5);
d = x0/(n1*scale);

*x = (x1+(d-(1.0+2.0*t1+c1)*pow(d,3.0)/6.0
+(5.0-2.0*c1+28.0*t1-3.0*c1*c1+8.0*ecc2+24.0*t1*t1)*pow(d,5.0)/120.0)/cos(phi))*RAD_DEG;
*y = (phi-n1*tan(phi)/r1*(d*d/2.0-(5.0+3.0*t1+10.0*c1-4.0*c1*c1-9.0*ecc2)
*pow(d,4.0)/24.0+(61.0+90.0*t1+298.0*c1+45.0*t1*t1-252.0*ecc2-3.0*c1*c1)*pow(d,6.0)/72.0))*RAD_DEG;
}

double mercator(double y, double a, double ecc)
{
if(y == 0.0)
{
return 0.0;
}
else
{
return a * (
( 1.0 - ecc/4.0 - 3.0*ecc*ecc/64.0 - 5.0*ecc*ecc*ecc/256.0 ) * y -
( 3.0*ecc/8.0 + 3.0*ecc*ecc/32.0 + 45.0*ecc*ecc*ecc/1024.0 ) * sin(2.0 * y) +
( 15.0*ecc*ecc/256.0 + 45.0*ecc*ecc*ecc/1024.0 ) * sin(4.0 * y) -
( 35.0*ecc*ecc*ecc/3072.0 ) * sin(6.0 * y) );
}
}

/*************************************************************************************

Sample code below, using coordinate in Dan Jacob's website.

*/

#include <stdio.h>
#include <stdlib.h>

int main()
{
double x1, y1, z1, x2, y2, z2;
double tx1, ty1, tx2, ty2;
double dx, dy, dz, dx1, dy1;

x1 = 120.85788004; // TWD67
y1 = 24.18347242;
z1 = 777;

x2 = 120.86603958; // TWD97
y2 = 24.18170479;
z2 = 777;

tx1 = 235561; // TWD67->TM2
ty1 = 2675359;

tx2 = 236389.849; // TWD97->TM2
ty2 = 2675153.168;

/////////////////////////////////////////////
//
//
// convert TWD67->TM2
//
//
/////////////////////////////////////////////

dx = x1;
dy = y1;

toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);
// center longitude of taiwan is 121, for penghu is 119

dx += 250000; // TM2 in Taiwan should add 250000

printf("TWD67->TM2\nTWD67 (%f, %f)\nConvert (%.3f, %.3f)\nOrigin (%.3f, %.3f)\n", x1, y1, dx, dy, tx1, ty1);
printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)\n\n", sqrt((dx-tx1)*(dx-tx1)+(dy-ty1)*(dy-ty1)), (dx-tx1), (dy-ty1));

/////////////////////////////////////////////
//
//
// convert TWD97->TM2
//
//
/////////////////////////////////////////////

dx = x2;
dy = y2;

toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);
// center longitude of taiwan is 121, for penghu is 119

dx += 250000; // TM2 in Taiwan should add 250000

printf("TWD97->TM2\nTWD97 (%f, %f)\nConvert (%.3f, %.3f)\nOrigin (%.3f, %.3f)\n", x2, y2, dx, dy, tx2, ty2);
printf("Acuuracy (%.3f, X:%.3f, Y:%.3f)\n\n", sqrt((dx-tx2)*(dx-tx2)+(dy-ty2)*(dy-ty2)), (dx-tx2), (dy-ty2));

/////////////////////////////////////////////
//
//
// convert TM2->TWD67
//
//
/////////////////////////////////////////////

dx = tx1-250000; // should minus 250000 first in Taiwan
dy = ty1;

fromTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx, &dy);

printf("TM2->TWD67\nTM2 (%f, %f)\nConvert (%.9f, %.9f)\nOrigin (%.9f, %.9f)\n", tx1, ty1, dx, dy, x1, y1);
printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)\n\n", sqrt((dx-x1)*(dx-x1)+(dy-y1)*(dy-y1)), (dx-x1), (dy-y1));

/////////////////////////////////////////////
//
//
// convert TM2->TWD97
//
//
/////////////////////////////////////////////

dx = tx2-250000; // should minus 250000 first in Taiwan
dy = ty2;

fromTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx, &dy);

printf("TM2->TWD97\nTM2 (%f, %f)\nConvert (%.9f, %.9f)\nOrigin (%.9f, %.9f)\n", tx2, ty2, dx, dy, x2, y2);
printf("Acuuracy (%.9f, X:%.9f, Y:%.9f)\n\n", sqrt((dx-x2)*(dx-x2)+(dy-y2)*(dy-y2)), (dx-x2), (dy-y2));

/////////////////////////////////////////////
//
//
// convert TWD67->TWD97
//
//
/////////////////////////////////////////////

dx = x1;
dy = y1;
dz = z1;

toTWD97(&dx, &dy, &dz);

dx1 = dx;
dy1 = dy;

toTM2(TWD97_A, TWD97_ECC, TWD97_ECC2, 0, 121, TWD97_TM2, &dx1, &dy1);

dx1 += 250000; // TM2 in Taiwan should add 250000

printf("TWD67->TWD97\nTWD67 (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x1, y1, z1, tx1, ty1);
printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", dx, dy, dz, dx1, dy1);
printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x2, y2, z2, tx2, ty2);
printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)\n\n", sqrt((dx1-tx2)*(dx1-tx2)+(dy1-ty2)*(dy1-ty2)), (dx1-tx2), (dy1-ty2));

/////////////////////////////////////////////
//
//
// convert TWD97->TWD67
//
//
/////////////////////////////////////////////

dx = x2;
dy = y2;
dz = z2;

toTWD67(&dx, &dy, &dz);

dx1 = dx;
dy1 = dy;

toTM2(TWD67_A, TWD67_ECC, TWD67_ECC2, 0, 121, TWD67_TM2, &dx1, &dy1);

dx1 += 250000; // TM2 in Taiwan should add 250000

printf("TWD97->TWD67\nTWD97 (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x2, y2, z2, tx2, ty2);
printf("Convert (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", dx, dy, dz, dx1, dy1);
printf("Origin (%.9f, %.9f, %6.2f) (%.3f, %.3f)\n", x1, y1, z1, tx1, ty1);
printf("Acuuracy (%.4f, X:%.4f, Y:%.4f)\n\n", sqrt((dx1-tx1)*(dx1-tx1)+(dy1-ty1)*(dy1-ty1)), (dx1-tx1), (dy1-ty1));
}

意見

Chauncey 提到...

#define COS67_5 0.3826834323650897717284599840304e0
請問這裡宣告常數,最後的e0有意義嗎?可以這樣宣告?在下要改成VB的寫法!

提到...

這個常數, 其實就是預先定義好的 cos(67.5°)

至於e , 代表這個數字是以科學記號來做表示的. e0是指10的零次方, 也就是數字要乘以1.

請參考維基百科中關於科學記號的相關說明.

匿名 提到...

Hello,
我把你的原始碼改寫成VB之後,用來轉換全國三角點的資料(TWD67TM2→TWD97DMS),但是發現越遠離虎子山,偏差量會越多,有沒有辦法修正?

詳情請見:
http://www.mobile01.com/topicdetail.php?f=228&t=506714&p=1

可從84樓開始:
http://www.mobile01.com/topicdetail.php?f=228&t=506714&p=6#5111185

提到...

這個問題主要是以下幾個常數的設定所造成的問題.

TWD67_DX
TWD67_DY
TWD67_DZ
TWD67_RX
TWD67_RY
TWD67_RZ
TWD67_S

目前所使用的數值, 是許久之前寫地圖程式時, 根據當時手
上的全台灣 TWD67 資料與 WGS84 資料比對而來, 確實是有
一定程度的變形. 不過當時內政部不肯開放數據, 只能自行
交叉比對.

你可以試試使用其他幾組不同的數據. 數據請參考
http://blog.minstrel.idv.tw/2004/06/taiwan-datum-parameter.html
的第 4,5,6 項. 或許結果會好一些.

小喵 提到...

你好~

我正在找TM2轉經緯度的轉換程式和公式~

剛好逛到這裡~有些問題想請教一下

TM2->TWD67和TM2->TWD97中的相關變數~

各代表的是什麼?!有沒有轉換公式?

不好意思~能不能幫我解個答呢?

謝謝噢^^

提到...

原始公式的由來可以參考 http://blog.minstrel.idv.tw/2004/05/datum-transformation.html 中七參數法的部份.

至於TM2/TWD67/TWD97轉換公式與相關變數的定義在這整份原始碼中應該說明得蠻清楚了, 不很確定你不瞭解的是哪一部份?