Taiwan Datum (TWD67, TWD97, TM2 Source)
台灣座標系統轉換原始碼 Taiwan Datum Convert Source(TWD67, TWD97, TM2 Source)
/*
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));
}
Related Site:
洪朝貴的首頁 (Chinese)
積丹尼 ＝ Dan Jacobson (Chinese/English)
Dan Jacobson (English)
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));
}
意見
#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轉換公式與相關變數的定義在這整份原始碼中應該說明得蠻清楚了, 不很確定你不瞭解的是哪一部份?
Most drawback playing councils are nonprofits funded by casinos, lotteries and sportsbooks. In almost each one of the 33 state legislatures that have legalized sports playing, the bills mandated the playing operators themselves present funding for drawback playing councils. Mike pushes for her opinion, so Brown sidesteps by mentioning some of the the} things would possibly be} 카지노사이트 commonest to drawback playing. Is playing fun and entertaining, or a compulsion that you just become preoccupied with? She'll often inform folks verify out|to take a look at} Gamblers Anonymous' 20-question survey, which covers some of the the} similar topics and features a statement that the majority all} drawback gamblers answer sure to seven or extra questions. Within 30 days after the event, the group shall submit a monetary report that lists all the receipts and expenditures for the Casino Night event.
張貼留言