墨卡托投影實現


又稱正軸等角圓柱投影。圓柱投影的一種,由荷蘭地圖學家墨卡托(G. Mercator)於1569年創擬。為地圖投影方法中影響最大的。
設想一個與地軸方向一致的圓柱切於或割於地球,按等角條件將經緯網投影到圓柱面上,將圓柱面展為平面后,得平面經緯線網。投影后經線是一組豎直的等距離平 行直線,緯線是垂直於經線的一組平行直線。各相鄰緯線間隔由赤道向兩極增大。一點上任何方向的長度比均相等,即沒有角度變形,而面積變形顯著,隨遠離標准 緯線而增大。該投影具有等角航線被表示成直線的特性,故廣泛用於編制航海圖和航空圖等。
墨卡托投影在切圓柱投影與割圓柱投影中,最早也是最常用的是切圓柱投影。

 

(http://baike.baidu.com/view/301981.htm)

 

公式參數

 

a -- 橢球體長半軸

b -- 橢球體短半軸

f -- 扁率 (a-b) /a

e -- 第一偏心率  

e’-- 第二偏心率  

N -- 卯酉圈曲率半徑

 

R -- 子午圈曲率半徑

 

B -- 緯度,L -- 經度,單位弧度(RAD)

 -- 縱直角坐標,  -- 橫直角坐標,單位米(M)

 

橢球體參數

我國常用的3 個橢球體參數如下(源自“全球定位系統測量規范 GB/T 18314-2001”):

 

橢球體

長半軸 a(米)

短半軸b(米)

Krassovsky(北京54 采用)

6378245

6356863.0188

IAG 75(西安80 采用)

6378140

6356755.2882

WGS 84

6378137

6356752.3142

 

 

墨卡托投影正反解公式

墨卡托投影正解公式:(B,L)→(X,Y),標准緯度B0,原點緯度 0,原點經度L0

 

墨卡托投影反解公式:(X,Y) →(B,L),標准緯度B0,原點緯度 0,原點經度L0

 

 

公式中EXP 為自然對數底,緯度通過迭代計算很快就收斂了。

 

程序實現

 

投影的實現封裝於一個類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<<endl;
// 7248377.351067:11578353.630128
 
       latS=123456;
       lgtS=654321;
 
m_mp.FromProj(latS,lgtS,latD,lgtD);
       latD=RadToDegree(latD);
       lgtD=RadToDegree(lgtD);
 
cout<< latD<<”:”<< lgtD<<endl;
// 1.288032: 6.781493

參考材料: 

1、《常用地圖投影轉換公式》 青島海洋地質研究所 戴勤奮

2、百度百科

 

 原文鏈接:墨卡托投影實現


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM