墨卡托投影实现
烤鱼片@eii.dlmu)
cleverysm@163.com
又称正轴等角圆柱投影。圆柱投影的一种，由荷兰地图学家墨卡托（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::SetABdouble a, double b)
{
ifa<=0||b<=0)
{
return;
}
__A=a;
__B=b;
}
//设定__B0
void MercatorProj::SetB0double b0)
{
ifb0<-PI/2||b0>PI/2)
{
return;
}
__B0=b0;
}
//设定__L0
void MercatorProj::SetL0double l0)
{
ifl0<-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::ToProjdouble B, double L, double &X, double &Y)
{
double f/*扁率*/,e/*第一偏心率*/,e_/*第二偏心率*/,NB0/*卯酉圈曲率半径*/,K,dtemp;
double E=exp1);
ifL<-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);
ifdtemp<0)
{
return 1;
}
e= sqrtdtemp);
dtemp=__A/__B)*__A/__B)-1;
ifdtemp<0)
{
return 1;
}
e_= sqrtdtemp);
NB0=__A*__A)/__B)/sqrt1+e_*e_*cos__B0)*cos__B0));
K=NB0*cos__B0);
Y=K*L-__L0);
X=K*logtanPI/4+B/2)*pow1-e*sinB))/1+e*sinB)),e/2));
return 0;
}
/*******************************************
投影反向转换程序
double X: 纵向直角坐标
double Y: 横向直角坐标
double& B: 维度,弧度
double& L: 经度,弧度
*******************************************/
int MercatorProj::FromProjdouble X, double Y, double& B, double& L)
{
double f/*扁率*/,e/*第一偏心率*/,e_/*第二偏心率*/,NB0/*卯酉圈曲率半径*/,K,dtemp;
double E=exp1);
if__A<=0||__B<=0)
{
return 1;
}
f =__A-__B)/__A;
dtemp=1-__B/__A)*__B/__A);
ifdtemp<0)
{
return 1;
}
e= sqrtdtemp);
dtemp=__A/__B)*__A/__B)-1;
ifdtemp<0)
{
return 1;
}
e_= sqrtdtemp);
NB0=__A*__A)/__B)/sqrt1+e_*e_*cos__B0)*cos__B0));
K=NB0*cos__B0);
L=Y/K+__L0;
B=__IterativeValue;
forint i=0;i<__IterativeTimes;i++)
{
B=PI/2-2*atanpowE,-X/K))*powE,e/2)*log1-e*sinB))/1+e*sinB)))));
}
return 0;
}
另需几个常量和函数：
//圆周率
const double PI=3.1415926535897932;
//角度到弧度的转换
double DegreeToRaddouble degree)
{
return PI*double)degree/double)180);
}
//弧度到角度的转换
double RadToDegreedouble rad)
{
return 180*rad)/PI;
}
测试主函数：
double b0,l0;
double latS,lgtS,latD,lgtD;
b0=30;
l0=0;
latS=60;
lgtS=120;
m_mp.SetAB6378137, 6378245,6378140); // WGS 84
m_mp.SetB0DegreeToRadb0));
m_mp.SetL0DegreeToRadl0));
m_mp.ToProjDegreeToRadlatS),DegreeToRadlgtS),latD,lgtD);
cout<< latD<<”:”<< lgtD<<endl;
// 7248377.351067:11578353.630128
latS=123456;
lgtS=654321;
m_mp.FromProjlatS,lgtS,latD,lgtD);
latD=RadToDegreelatD);
lgtD=RadToDegreelgtD);
cout<< latD<<”:”<< lgtD<<endl;
// 1.288032: 6.781493
参考材料：
《常用地图投影转换公式》 青岛海洋地质研究所 戴勤奋
百度百科