Mercator Projection Implementation

Also known as the positive axis conformal cylindrical projection. A cylindrical projection of a kind, by the Dutch map biologist Mercator (G. Mercator) was created in 1569. Is the most influential in the map projection method.
Imagine a cylinder aligned with the axis of the Earth cut or cut in the world, according to the isometric conditions of the graticule projection to the cylindrical surface, the cylindrical surface to show as a plane, the plane longitude net. The post-projection meridian is a set of vertical equidistant parallel lines, and parallels are a set of parallel lines perpendicular to the meridians. The interval between the adjacent parallels increases from the equator to the poles. The length ratio of any direction is equal, that is, there is no angular deformation, and the area deformation is significant, and increases with the distance from the standard parallels. The projection has the characteristics of the conformal route being expressed as a straight line, so it is widely used in compiling nautical charts and aerial diagrams.
The first and most commonly used projection of the Mercator projection in the tangent cylindrical projection and the cylindrical projection is the tangent cylindrical projections.

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

Formula Parameters

A--ellipsoid body length half axis

B--ellipsoid short half axis

F--flat rate (A-B)/A

E--First eccentricity 

E '--second eccentricity 

N--The curvature radius of the unitary circle of the Mao

R--meridional curvature radius

B--Latitude, L--longitude, Unit radian (RAD)

--longitudinal coordinate,--transverse rectangular coordinate, unit m (m)

parameters of ellipsoidal body

The 3 ellipsoid parameters commonly used in China are as follows (from "Global Positioning System measurement standard GB/T 18314-2001"):

Ellipsoidal body	Long half shaft A (m)	Short half axle B (m)
Krassovsky (Beijing 54 adopted)	6378245	6356863.0188
IAG 75 (XI ' an 80 adoption)	6378140	6356755.2882
WGS 84	6378137	6356752.3142

The formula of positive and negative solution of Mercator projection

Mercator Positive Solution formula:(b,l) → (x, y), standard latitude B0, origin latitude 0, origin longitude L0

Mercator projection Inverse formula: (x, y) → (b,l), standard latitude B0, origin latitude 0, origin longitude L0

Exp is a natural logarithmic base in the formula, and latitude B converges quickly by iterative computation.

Program Implementation

The implementation of the projection is encapsulated in a class mercatorproj.

Defines several private variables in a class, saving related parameters

int __iterativetimes; Number of iterations in the reverse conversion program

Double __iterativevalue; Iterative initial values in a reverse translator

Double __a; ellipsoid body length half axis, m

Double __b; ellipsoid Short half axis, m

Double __b0; Standard Latitude, radians

Double __l0; Origin longitude, radians

The setting of the above parameters is done by the following public functions

Setting __a and __b

void Mercatorproj::setab (Double A, double b)

{

if (a<=0| | b<=0)

{

Return

}

__a=a;

__b=b;

}

Set __b0

void Mercatorproj::setb0 (double b0)

{

if (b0<-pi/2| | B0>PI/2)

{

Return

}

__b0=b0;

}

Set __l0

void Mercatorproj::setl0 (double l0)

{

if (l0<-pi| | L0>PI)

{

Return

}

__l0=l0;

}

Give parameter default value in constructor

Mercatorproj::mercatorproj ()

{

__iterativetimes=10; Number of iterations is 10

__iterativevalue=0; Iteration Initial value

__b0=0;

__l0=0;

__a=1;

__b=1;

}

/*******************************************

Projected forward conversion program

Double B: Dimension, radians

Double L: Longitude, Radian

double& X: longitudinal rectangular coordinates

double& Y: Transverse rectangular coordinates

*******************************************/

int Mercatorproj::toproj (double B, double L, double &x, double &y)

{

Double f/* flat rate */,e/* first eccentricity */,e_/* second eccentricity */,nb0/* Mao unitary circle curvature radius */,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;

}

/*******************************************

Projection Reverse Conversion Program

Double X: longitudinal rectangular coordinates

Double Y: Transverse rectangular coordinates

double& B: Dimensions, radians

double& L: Longitude, Radian

*******************************************/

int Mercatorproj::fromproj (double X, double Y, double& B, double& L)

{

Double f/* flat rate */,e/* first eccentricity */,e_/* second eccentricity */,nb0/* Mao unitary circle curvature radius */,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;

}

several constants and functions are required:

Pi

const double pi=3.1415926535897932;

Angle-to-radian conversion

Double Degreetorad (double degree)

{

Return pi* (Double) degree/(double) 180);

}

ARC-to-angle conversion

Double Radtodegree (Double rad)

{

Return (180*rad)/pi;

}

To test the main function:

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;

72,483,77.,351,067:115,783,53.,630,128

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

 

Reference Material:

"Common map Projection Transformation Formula" Qingdao Institute of Marine Geology Dai Diligence