樹心幽徑

C GPX Reader Project STEP1: 20170423將WGS84經緯度座標轉TWD97(TWD67)二度分帶座標
2017/04/17,14:00

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

 

 
Accessible and Valid XHTML 1.0 Strict and CSS Powered by LifeType