GPS坐标转换代码

经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法收藏新一篇:C#的6种常用集合类大比拼|旧一篇:Silverlight中使用图片及常见问题一、经纬度BL换算到高斯平面直角坐标XY（高斯投影正算）的源码及算法//GaussBL2xy.cpp:Definestheentrypointfortheconsoleapplication.//#includestdafx.h#includemath.h#includeCoorTrans.h#includeiostreamusingnamespacestd;voidmain(intargc,char*argv[]){doubleMyL0=108;//中央子午线doubleMyB=33.44556666;//33du44fen55.6666miaodoubleMyL=109.22334444;//3度带，109d22m33.4444sPrjPoint_KrasovskyMyPrj;MyPrj.SetL0(MyL0);MyPrj.SetBL(MyB,MyL);doubleOutMyX;//结果应该大致是:3736714.783,627497.303doubleOutMyY;OutMyX=MyPrj.x;//正算结果:北坐标xOutMyY=MyPrj.y;//结果：东坐标y//////////////////反算////////////////////////////////////////doubleInputMyX=3736714.783;//如果是独立计算，应该给出中央子午线L0doubleInputMyY=627497.303;MyPrj.Setxy(InputMyX,InputMyY);MyPrj.GetBL(&MyPrj.B,&MyPrj.L);//把计算出的BL的弧度值换算为dms形式doubleOutputMyB;doubleOutputMyL;OutputMyB=MyPrj.B;//反算结果：BOutputMyL=MyPrj.L;//反算结果：L//分析表明，此程序的结果和Coord4.2的转换结果是一样的，只差到毫米级//原程序有几个问题，1.Pi的值不对。2.SetBL中多了两行错误代码}doubleDms2Rad(doubleDms){doubleDegree,Miniute;doubleSecond;intSign;doubleRad;if(Dms=0)Sign=1;elseSign=-1;Dms=fabs(Dms);Degree=floor(Dms);Miniute=floor(fmod(Dms*100.0,100.0));Second=fmod(Dms*10000.0,100.0);Rad=Sign*(Degree+Miniute/60.0+Second/3600.0)*PI/180.0;returnRad;}doubleRad2Dms(doubleRad){doubleDegree,Miniute;doubleSecond;intSign;doubleDms;if(Rad=0)Sign=1;elseSign=-1;Rad=fabs(Rad*180.0/PI);Degree=floor(Rad);Miniute=floor(fmod(Rad*60.0,60.0));Second=fmod(Rad*3600.0,60.0);Dms=Sign*(Degree+Miniute/100.0+Second/10000.0);returnDms;}/////////////////////////////////////////////////////DefinitionofPrjPoint///////////////////////////////////////////////////boolPrjPoint::BL2xy(){doubleX,N,t,t2,m,m2,ng2;doublesinB,cosB;X=A1*B*180.0/PI+A2*sin(2*B)+A3*sin(4*B)+A4*sin(6*B);sinB=sin(B);cosB=cos(B);t=tan(B);t2=t*t;N=a/sqrt(1-e2*sinB*sinB);m=cosB*(L-L0);m2=m*m;ng2=cosB*cosB*e2/(1-e2);//x,y的计算公式见孔祥元等主编武汉大学出版社2002年出版的《控制测量学》的第72页//书的的括号有问题，(和[应该交换x=X+N*t*((0.5+((5-t2+9*ng2+4*ng2*ng2)/24.0+(61-58*t2+t2*t2)*m2/720.0)*m2)*m2);y=N*m*(1+m2*((1-t2+ng2)/6.0+m2*(5-18*t2+t2*t2+14*ng2-58*ng2*t2)/120.0));y+=500000;returntrue;}boolPrjPoint::xy2BL(){doublesinB,cosB,t,t2,N,ng2,V,yN;doublepreB0,B0;doubleeta;y-=500000;B0=x/A1;do{preB0=B0;B0=B0*PI/180.0;B0=(x-(A2*sin(2*B0)+A3*sin(4*B0)+A4*sin(6*B0)))/A1;eta=fabs(B0-preB0);}while(eta0.000000001);B0=B0*PI/180.0;B=Rad2Dms(B0);sinB=sin(B0);cosB=cos(B0);t=tan(B0);t2=t*t;N=a/sqrt(1-e2*sinB*sinB);ng2=cosB*cosB*e2/(1-e2);V=sqrt(1+ng2);yN=y/N;B=B0-(yN*yN-(5+3*t2+ng2-9*ng2*t2)*yN*yN*yN*yN/12.0+(61+90*t2+45*t2*t2)*yN*yN*yN*yN*yN*yN/360.0)*V*V*t/2;L=L0+(yN-(1+2*t2+ng2)*yN*yN*yN/6.0+(5+28*t2+24*t2*t2+6*ng2+8*ng2*t2)*yN*yN*yN*yN*yN/120.0)/cosB;returntrue;}boolPrjPoint::SetL0(doubledL0){L0=Dms2Rad(dL0);returntrue;}boolPrjPoint::SetBL(doubledB,doubledL){B=Dms2Rad(dB);L=Dms2Rad(dL);//B=dB;//我靠，Iwanasayfuck//L=dL;//delit!BL2xy();returntrue;}boolPrjPoint::GetBL(double*dB,double*dL){*dB=Rad2Dms(B);*dL=Rad2Dms(L);returntrue;}boolPrjPoint::Setxy(doubledx,doubledy){x=dx;y=dy;xy2BL();returntrue;}boolPrjPoint::Getxy(double*dx,double*dy){*dx=x;*dy=y;returntrue;}/////////////////////////////////////////////////////DefinitionofPrjPoint_IUGG1975///////////////////////////////////////////////////PrjPoint_IUGG1975::PrjPoint_IUGG1975()//在类外定义构造成员函数，要加上类名和域限定符::{a=6378140;f=298.257;e2=1-((f-1)/f)*((f-1)/f);e12=(f/(f-1))*(f/(f-1))-1;A1=111133.0047;//这几个A是什么意思？A2=-16038.5282;A3=16.8326;A4=-0.0220;}PrjPoint_IUGG1975::~PrjPoint_IUGG1975()//析构函数，释放构造函数占用的内存{}/////////////////////////////////////////////////////DefinitionofPrjPoint_Krasovsky///////////////////////////////////////////////////PrjPoint_Krasovsky::PrjPoint_Krasovsky(){a=6378245;f=298.3;e2=1-((f-1)/f)*((f-1)/f);e12=(f/(f-1))*(f/(f-1))-1;A1=111134.8611;A2=-16036.4803;A3=16.8281;A4=-0.0220;}PrjPoint_Krasovsky::~PrjPoint_Krasovsky(){}//CoorTrans.h文件#ifndef_COORTRANS_H_INCLUDED#define_COORTRANS_H_INCLUDED#includestdafx.h#includemath.hconstdoublePI=3.141592653589793;doubleDms2Rad(doubleDms);doubleRad2Dms(doubleRad);classPrjPoint{public:doubleL0;//中央子午线经度doubleB,L;//大地坐标doublex,y;//高斯投影平面坐标public:boolBL2xy();boolxy2BL();protected:doublea,f,e2,e12;//基本椭球参数doubleA1,A2,A3,A4;//用于计算X的椭球参数public:boolSetL0(doubledL0);boolSetBL(doubledB,doubledL);boolGetBL(double*dB,double*dL);boolSetxy(doubledx,doubledy);boolGetxy(double*dx,double*dy);};classPrjPoint_Krasovsky:virtualpublicPrjPoint//类的继承，声明基类是PrjPoint，虚基类{public:PrjPoint_Krasovsky();//定义构造函数，用来初始化。(函数名和类名相同)~PrjPoint_Krasovsky();};classPrjPoint_IUGG1975:virtualpublicPrjPoint{public:PrjPoint_IUGG1975();~PrjPoint_IUGG1975();};#endif/*ndef_COORTRANS_H_INCLUDED*/