又稱正軸等角圓柱投影。圓柱投影的一種,由荷蘭地圖學家墨卡托(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 為自然對數底,緯度B 通過迭代計算很快就收斂了。
程序實現
投影的實現封裝於一個類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、百度百科
原文鏈接:墨卡托投影實現
