【轉】GPS經緯度資料轉換到以米單位的平面座標系統

標籤：

以前在網上找的，現在不知道出處是哪裡了，呵呵...

轉載出來。

  1 //笛卡爾座標系    2 typedef struct tagCRDCARTESIAN{    3 double x;    4 double y;    5 double z;    6 }CRDCARTESIAN,*PCRDCARTESIAN;    7 //typedef CRDCARTESIAN *PCRDCARTESIAN;    8 //大地座標系    9 typedef struct tagCRDGEODETIC{   10 double longitude; //經度   11 double latitude;  //緯度    12 double height;    //大地高,可設為0   13 }CRDGEODETIC;   14 typedef CRDGEODETIC *PCRDGEODETIC;   15    16 void CoordCovert::GeodeticToCartesian (PCRDCARTESIAN pcc, PCRDGEODETIC pcg,   17 double dSemiMajorAxis, double dFlattening)   18 {   19     double B;    //緯度度數   20     double L;    //經度度數   21     double L0;    //中央經線度數   22     double l;    //L-L0   23     double t;    //tanB   24     double m;    //ltanB   25     double N;    //卯酉圈曲率半徑    26     double q2;   27     double x;    //高斯平面縱座標   28     double y;    //高斯平面橫座標   29     double s;    //赤道至緯度B的經線弧長   30     double f;    //參考橢球體扁率   31     double e2;    //橢球第一偏心率   32     double a;    //參考橢球體長半軸   33     //double b;    //參考橢球體短半軸   34     double a1;   35     double a2;   36     double a3;   37     double a4;   38     double b1;   39     double b2;   40     double b3;   41     double b4;   42     double c0;   43     double c1;   44     double c2;   45     double c3;   46     int Datum=84;    //投影基準面類型：北京54基準面為54，西安80基準面為80，WGS84基準面為84   47     int prjno=0;    //投影帶號   48     int zonewide=3;       49     double IPI=0.0174532925199433333333;    //3.1415926535898/180.0   50     B=pcg->latitude ; //緯度   51     L=pcg->longitude ; //經度   52     if (zonewide==6)    53     {   54          prjno=(int)(L/zonewide)+1;   55          L0=prjno*zonewide-3;   56     }   57     else   58     {   59         prjno=(int)((L-1.5)/3)+1;   60         L0=prjno*3;   61     }   62        63     if(Datum==54)   64     {   65          a=6378245;   66          f=1/298.3;   67     }       68     else if(Datum==84)   69     {   70         a=6378137;   71         f=1/298.257223563;   72     }   73     L0=L0*IPI;   74     L=L*IPI;   75     B=B*IPI;   76     e2=2*f-f*f;//(a*a-b*b)/(a*a);   77     l=L-L0;   78     t=tan(B);   79     m=l * cos(B);   80     N=a/sqrt(1-e2* sin(B) * sin(B));   81     q2=e2/(1-e2)* cos(B)* cos(B);   82     a1=1+(double)3/4*e2+(double)45/64*e2*e2+(double)175/256*e2*e2*e2+(double)11025/16384*e2*e2*e2*e2+(double)43659/65536*e2*e2*e2*e2*e2;   83     a2=(double)3/4*e2+(double)15/16*e2*e2+(double)525/512*e2*e2*e2+(double)2205/2048*e2*e2*e2*e2+(double)72765/65536*e2*e2*e2*e2*e2;   84     a3=(double)15/64*e2*e2+(double)105/256*e2*e2*e2+(double)2205/4096*e2*e2*e2*e2+(double)10359/16384*e2*e2*e2*e2*e2;   85     a4=(double)35/512*e2*e2*e2+(double)315/2048*e2*e2*e2*e2+(double)31185/13072*e2*e2*e2*e2*e2;   86     b1=a1*a*(1-e2);   87     b2=(double)-1/2*a2*a*(1-e2);   88     b3=(double)1/4*a3*a*(1-e2);   89     b4=(double)-1/6*a4*a*(1-e2);   90     c0=b1;   91     c1=2*b2+4*b3+6*b4;   92     c2=-(8*b3+32*b4);   93     c3=32*b4;   94     s=c0*B+cos(B)*(c1*sin(B)+c2*sin(B)*sin(B)*sin(B)+c3*sin(B)*sin(B)*sin(B)*sin(B)*sin(B));   95     x=s+(double)1/2*N*t*m*m+(double)1/24*(5-t*t+9*q2+4*q2*q2)*N*t*m*m*m*m+(double)1/720*(61-58*t*t+t*t*t*t)*N*t*m*m*m*m*m*m;   96     y=N*m+(double)1/6*(1-t*t+q2)*N*m*m*m+(double)1/120*(5-18*t*t+t*t*t*t-14*q2-58*q2*t*t)*N*m*m*m*m*m;   97     y=y+1000000*prjno+500000;   98     pcc->x=x;   99     pcc->y=y-38000000;  100     pcc->z=0;  101 }  

 

【轉】GPS經緯度資料轉換到以米單位的平面座標系統