转GPS经纬度数据转换到以米单位的平面坐标系

Posted 2020-06-16

以前在网上找的，现在不知道出处是哪里了，呵呵...

转载出来。

  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 }