dog@wutingan-All-Series:~/lake$ g++ i.c
dog@wutingan-All-Series:~/lake$ ./a.out
將WGS84經緯度座標轉TWD97,TWD67
台灣虎子山經緯度WGS84=120.982025,23.973875 其二度分帶座標TWD97=248170.787,2652129.936 TWD67=247341.925, 2652335.811
請輸入經度(例如虎子山lon=120.982025)=120.982025
請輸入緯度(例如虎子山lat=23.973875)=23.973875
轉換結果如下
WGS84(120.982025,23.973875)
TWD97(248170.825725,2652129.977372)
TWD67(247341.887019,2652335.877557)
dog@wutingan-All-Series:~/lake$ cat i.c
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
int lonlat2twd97(double lon,double lat,double *x, double *y)
{
double a,b,long0,k0,dx;
double e,e2,n,nu,p;
double A,B,C,D,E,S;
double K1,K2,K3;
double K4,K5;
a = 6378137.0;
b = 6356752.314245;
long0 = 121.0*M_PI/180.0;
k0 = 0.9999;
dx = 250000;
e = pow(1-(b*b)/(a*a),0.5);
e2 = pow(e,2)/(1-pow(e,2));
n = (a-b)/(a+b);
nu = a/pow(1-pow(e,2)*pow(sin(lat),2),0.5);
p = lon-long0;
A = a*(1 - n + (5/4.0)*(pow(n,2) - pow(n,3)) + (81/64.0)*(pow(n,4) - pow(n,5)));
B = (3*a*n/2.0)*(1 - n + (7/8.0)*(pow(n,2) - pow(n,3)) + (55/64.0)*(pow(n,4) - pow(n,5)));
C = (15*a*pow(n,2)/16.0)*(1 - n + (3/4.0)*(pow(n,2) - pow(n,3)));
D = (35*a*pow(n,3)/48.0)*(1 - n + (11/16.0)*(pow(n,2) - pow(n,3)));
E = (315*a*pow(n,4)/51.0)*(1 - n);
S = A*lat - B*sin(2*lat) + C*sin(4*lat) - D*sin(6*lat) + E*sin(8*lat);
K1 = S*k0;
K2 = k0*nu*sin(2*lat)/4.0;
K3 = (k0*nu*sin(lat)*pow(cos(lat),3)/24.0) * (5.0 - pow(tan(lat),2) + 9.0*e2*pow(cos(lat),2) + 4.0*pow(e2,2)*pow(cos(lat),4));
*y = K1 + K2*pow(p,2) + K3*pow(p,4);
K4 = k0*nu*cos(lat);
K5 = (k0*nu*pow(cos(lat),3)/6.0) * (1 - pow(tan(lat),2) + e2*pow(cos(lat),2));
*x = K4*p + K5*pow(p,3) + dx;
}
int lonlat2twd67(double lon,double lat,double *x, double *y)
{
double x97,y97;
lonlat2twd97(lon,lat,&x97,&y97);
*x=x97-807.8-0.00001549*x97-0.000006521*y97;
*y=y97+248.6-0.00001549*y97-0.000006521*x97;
}
int main()
{
double lon,lat,londeg,latdeg;
printf("將WGS84經緯度座標轉TWD97,TWD67\n");
printf("台灣虎子山經緯度WGS84=120.982025,23.973875 其二度分帶座標TWD97=248170.787,2652129.936 TWD67=247341.925, 2652335.811\n");
printf("請輸入經度(例如虎子山lon=120.982025)=");scanf("%lf",&londeg);
printf("請輸入緯度(例如虎子山lat=23.973875)=");scanf("%lf",&latdeg);
lon=londeg*M_PI/180.0; lat=latdeg*M_PI/180.0;
double x97,y97,x67,y67;
lonlat2twd97(lon,lat,&x97,&y97);
lonlat2twd67(lon,lat,&x67,&y67);
printf("轉換結果如下\n");
printf("WGS84(%f,%f)\n",londeg,latdeg);
printf("TWD97(%f,%f)\n",x97,y97);
printf("TWD67(%f,%f)\n",x67,y67);
}
//參考程式設計遇上小提琴BLOG: http://blog.ez2learn.com/2009/08/15/lat-lon-to-twd97/
http://www.kmvs.km.edu.tw/lf/index.php?op=ViewAlbum&albumId=247&blogId=70
|