墨卡托投影

墨卡托投影 投影的封装于一个实实实实实实实实class MercatorProj中。实实实实实实实实实实实实实实实中定若干私有量,保存相参数 int __IterativeTimes; //反向程序中的迭代次数实实实实实实实实实实 double __IterativeValue; //反向程序中的迭代初始实实实实实实实实实实实 double __A; //实实实实球体半,米 double __B; //实实球体短半,米 double __B0; //实实实准度,弧度 double __L0; //原点度实实,弧度 以上参数的定由如下几个实实实实实实实public函数完成//实定__A与__B void MercatorProj::SetAB(double a, double b) { if(a<=0||b<=0) { return; } __A=a; __B=b; } //实定__B0 void MercatorProj::SetB0(double b0) { if(b0<-PI/2||b0>PI/2) { return; } __B0=b0; } //实定__L0 void MercatorProj::SetL0(double l0) { if(l0<-PI||l0>PI) { return; } __L0=l0; } //构造函数中予参数默实实实实实实实 MercatorProj::MercatorProj(){ __IterativeTimes=10; //迭代次数实10 __IterativeValue=0; //迭代初始实 __B0=0; __L0=0; __A=1; __B=1; } /******************************************* 投影正向程序实实实实 double B: 实度,弧度 double L: 实度,弧度 double& X: 实实向直角坐 double& Y: 横向直角坐实 *******************************************/ int MercatorProj::ToProj(double B, double L, double &X, double &Y) { double f/*扁率*/,e/*第一偏心率*/,e_/*第二偏心率*/,NB0/*卯酉圈曲率半径*/,K,dtemp; double E=exp(1); if(L<-PI||L>PI||B<-PI/2||B>PI/2) { return 1; } if(__A<=0||__B<=0) { return 1; } f =(__A-__B)/__A; dtemp=1-(__B/__A)*(__B/__A); if(dtemp<0) { return 1; } e= sqrt(dtemp); dtemp=(__A/__B)*(__A/__B)-1; if(dtemp<0) { return 1; } e_= sqrt(dtemp); NB0=((__A*__A)/__B)/sqrt(1+e_*e_*cos(__B0)*cos(__B0)); K=NB0*cos(__B0); Y=K*(L-__L0); X=K*log(tan(PI/4+B/2)*pow((1-e*sin(B))/(1+e*sin(B)),e/2)); return 0; } /*******************************************投影反向程序实实实实 double X: 实实向直角坐 double Y: 横向直角坐实 double& B: 实度,弧度 double& L: 实度,弧度 *******************************************/int MercatorProj::FromProj(double X, double Y, double& B, double& L) { double f/*扁率*/,e/*第一偏心率*/,e_/*第二偏心率*/,NB0/*卯酉圈曲率半径*/,K,dtemp; double E=exp(1); if(__A<=0||__B<=0) { return 1; } f =(__A-__B)/__A; dtemp=1-(__B/__A)*(__B/__A); if(dtemp<0) { return 1; } e= sqrt(dtemp); dtemp=(__A/__B)*(__A/__B)-1; if(dtemp<0) { return 1; } e_= sqrt(dtemp); NB0=((__A*__A)/__B)/sqrt(1+e_*e_*cos(__B0)*cos(__B0)); K=NB0*cos(__B0); L=Y/K+__L0; B=__IterativeValue; for(int i=0;i<__IterativeTimes;i++) { B=PI/2-2*atan(pow(E,(-X/K))*pow(E,(e/2)*log((1-e*sin(B))/(1+e*sin(B))))); } return 0; } 另需几个常量和函数: //实周率 const double PI=3.1415926535897932;//角度到弧度的实实 double DegreeToRad(double degree){ return PI*((double)degree/(double)180); } //弧度到角度的实实 double RadToDegree(double rad){ return (180*rad)/PI; } 实实主函数: double b0,l0; double latS,lgtS,latD,lgtD; b0=30; l0=0; latS=60; lgtS=120; m_mp.SetAB(6378137, 6378245,6378140); // WGS 84 m_mp.SetB0(DegreeToRad(b0)); m_mp.SetL0(DegreeToRad(l0)); m_mp.ToProj(DegreeToRad(latS),DegreeToRad(lgtS),latD,lgtD); cout<< latD<<”:”<< lgtD<