【转】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 }
