# Mercator

This article describes algorithms for Mercator projection. For the OSM editor program, see Merkaartor.

## Spherical Mercator

Most of OSM, including the main tiling system, uses a spherical Mercator projection. This produces a fast approximation to the truer, elliptical projection.

### ActionScript and JavaScript

The x co-ordinate is simply the longitude multiplied up by your chosen scale.

```function y2lat(a) { return 180/Math.PI * (2 * Math.atan(Math.exp(a*Math.PI/180)) - Math.PI/2); }
function lat2y(a) { return 180/Math.PI * Math.log(Math.tan(Math.PI/4+a*(Math.PI/180)/2)); }```

### Excel

I have simply converted the above script formulas into Excel formulas (value in A1): --Krza 23:36, 19 July 2008 (UTC)

```y2lat: =180/PI() * (2 * ARCTAN(EXP(A1*PI()/180)) - PI()/2)
lat2y: =180/PI() * LOG(TAN(PI()/4 + A1*(PI()/180)/2))
Note: this may work for older versions of Excel,
but in Excel 2010, ARCTAN is called ATAN, and you need to use LN instead of LOG to get the natural logarithm
```

### C

```#include <math.h>

/* The following functions take or return there results in degrees */

double y2lat_d(double y) { return rad2deg(2 * atan(exp(  deg2rad(y) ) ) - M_PI/2); }
double x2lon_d(double x) { return x; }
double lat2y_d(double lat) { return rad2deg(log(tan(M_PI/4+ deg2rad(lat)/2))); }
double lon2x_d(double lon) { return lon; }

/* The following functions take or return there results in something close to meters, along the equator */

double y2lat_m(double y) { return rad2deg(2 * atan(exp( (y / earth_radius ) )) - M_PI/2); }
double x2lon_m(double x) { return rad2deg(x / earth_radius); }
double lat2y_m(double lat) { return earth_radius * log(tan(M_PI/4+ deg2rad(lat)/2)); }
double lon2x_m(double lon) { return deg2rad(lon) * earth_radius; }```

### PostGIS / SQL

The SRID 900913 is commonly used for our projection but is not formally assigned. It is not included in the default projections supplied with PostGIS. It can be added into your database via:

``` INSERT INTO spatial_ref_sys (srid, auth_name, auth_srid, srtext, proj4text) VALUES
(900913,'EPSG',900913,'PROJCS["WGS84 / Simple Mercator",GEOGCS["WGS 84",
DATUM["WGS_1984",SPHEROID["WGS_1984", 6378137.0, 298.257223563]],PRIMEM["Greenwich", 0.0],
UNIT["degree", 0.017453292519943295],AXIS["Longitude", EAST],AXIS["Latitude", NORTH]],
PARAMETER["central_meridian", 0.0],PARAMETER["scale_factor", 1.0],PARAMETER["false_easting", 0.0],
PARAMETER["false_northing", 0.0],UNIT["m", 1.0],AXIS["x", EAST],
AXIS["y", NORTH],AUTHORITY["EPSG","900913"]]',
'+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +no_defs');```

### Java

```import java.lang.Math;

public class SphericalMercator {
public static double y2lat(double aY) {
return Math.toDegrees(2* Math.atan(Math.exp(Math.toRadians(aY))) - Math.PI/2);
}

public static double lat2y(double aLat) {
}
}```

### PHP

```function lon2x(\$lon) { return deg2rad(\$lon) * 6378137.0; }
function lat2y(\$lat) { return log(tan(M_PI_4 + deg2rad(\$lat) / 2.0)) * 6378137.0; }
function x2lon(\$x) { return rad2deg(\$x / 6378137.0); }
function y2lat(\$y) { return rad2deg(2.0 * atan(exp(\$y / 6378137.0)) - M_PI_2); }```

### LibreOffice Calc (3.5)

lon value in A1, lat value in A2

``` lon2x: = A1 * (PI()/180) * 6378137
x2lon: = A1 / (PI()/180) / 6378137
lat2y: = LN(TAN(A2 * (PI()/180)/2 + PI()/4)) * 6378137
y2lat: = (2*ATAN(EXP(A2/6378137)) - PI()/2) / (PI()/180```

### Python

```import math
def y2lat(a):
return 180.0/math.pi*(2.0*math.atan(math.exp(a*math.pi/180.0))-math.pi/2.0)
def lat2y(a):
return 180.0/math.pi*math.log(math.tan(math.pi/4.0+a*(math.pi/180.0)/2.0))```

### C# implementation

```		public double YToLatitude(double y)
{
return 180.0/System.Math.PI *
(2 *
System.Math.Atan(
System.Math.Exp(y*System.Math.PI/180)) - System.Math.PI/2);
}
public double LatitudeToY (double latitude)
{
return 180.0/System.Math.PI *
System.Math.Log(
System.Math.Tan(
System.Math.PI/4.0+latitude*(System.Math.PI/180.0)/2));
}```

FIXME: add spherical code in other languages here

## Elliptical Mercator

### JavaScript

From a posting by Christopher Schmidt to the dev list on 2nd December 2006:

So, for everyone's elucidation, here is the way to convert from latitude and longitude to a simple Mercator projection, in Javascript:

```function deg_rad(ang) {
return ang * (Math.PI/180.0)
}
function merc_x(lon) {
var r_major = 6378137.000;
return r_major * deg_rad(lon);
}
function merc_y(lat) {
if (lat > 89.5)
lat = 89.5;
if (lat < -89.5)
lat = -89.5;
var r_major = 6378137.000;
var r_minor = 6356752.3142;
var temp = r_minor / r_major;
var es = 1.0 - (temp * temp);
var eccent = Math.sqrt(es);
var phi = deg_rad(lat);
var sinphi = Math.sin(phi);
var con = eccent * sinphi;
var com = .5 * eccent;
con = Math.pow((1.0-con)/(1.0+con), com);
var ts = Math.tan(.5 * (Math.PI*0.5 - phi))/con;
var y = 0 - r_major * Math.log(ts);
return y;
}
function merc(x,y) {
return [merc_x(x),merc_y(y)];
}```

LatLon to/from Mercator converting class based on script above, and proj4 implementation

```var Conv=({
r_minor:6356752.314245179,//defined as constant
f:298.257223563,//1/f=(a-b)/a , a=r_major, b=r_minor
{
var r=d*(Math.PI/180.0);
return r;
},
{
var d=r/(Math.PI/180.0);
return d;
},
ll2m:function(lon,lat) //lat lon to mercator
{
//lat, lon in rad
var x=this.r_major * this.deg2rad(lon);

if (lat > 89.5) lat = 89.5;
if (lat < -89.5) lat = -89.5;

var temp = this.r_minor / this.r_major;
var es = 1.0 - (temp * temp);
var eccent = Math.sqrt(es);

var phi = this.deg2rad(lat);

var sinphi = Math.sin(phi);

var con = eccent * sinphi;
var com = .5 * eccent;
var con2 = Math.pow((1.0-con)/(1.0+con), com);
var ts = Math.tan(.5 * (Math.PI*0.5 - phi))/con2;
var y = 0 - this.r_major * Math.log(ts);
var ret={'x':x,'y':y};
return ret;
},
m2ll:function(x,y) //mercator to lat lon
{

var temp = this.r_minor / this.r_major;
var e = Math.sqrt(1.0 - (temp * temp));
var lat=this.rad2deg(this.pj_phi2( Math.exp( 0-(y/this.r_major)), e));

var ret={'lon':lon,'lat':lat};
return ret;
},
pj_phi2:function(ts, e)
{
var N_ITER=15;
var HALFPI=Math.PI/2;

var TOL=0.0000000001;
var eccnth, Phi, con, dphi;
var i;
var eccnth = .5 * e;
Phi = HALFPI - 2. * Math.atan (ts);
i = N_ITER;
do
{
con = e * Math.sin (Phi);
dphi = HALFPI - 2. * Math.atan (ts * Math.pow((1. - con) / (1. + con), eccnth)) - Phi;
Phi += dphi;

}
while ( Math.abs(dphi)>TOL && --i);
return Phi;
}
});
//usage
var mercator=Conv.ll2m(47.6035525, 9.770602);//output mercator.x, mercator.y
var latlon=Conv.m2ll(5299424.36041, 1085840.05328);//output latlon.lat, latlon.lon```

### C implementation

```#include <math.h>

/*
* Mercator transformation
* accounts for the fact that the earth is not a sphere, but a spheroid
*/
#define D_R (M_PI / 180.0)
#define R_D (180.0 / M_PI)
#define R_MAJOR 6378137.0
#define R_MINOR 6356752.3142
#define RATIO (R_MINOR/R_MAJOR)
#define ECCENT (sqrt(1.0 - (RATIO * RATIO)))
#define COM (0.5 * ECCENT)

static double deg_rad (double ang) {
return ang * D_R;
}

double merc_x (double lon) {
return R_MAJOR * deg_rad (lon);
}

double merc_y (double lat) {
lat = fmin (89.5, fmax (lat, -89.5));
double phi = deg_rad(lat);
double sinphi = sin(phi);
double con = ECCENT * sinphi;
con = pow((1.0 - con) / (1.0 + con), COM);
double ts = tan(0.5 * (M_PI * 0.5 - phi)) / con;
return 0 - R_MAJOR * log(ts);
}

static double rad_deg (double ang) {
return ang * R_D;
}

double merc_lon (double x) {
return rad_deg(x) / R_MAJOR;
}

double merc_lat (double y) {
double ts = exp ( -y / R_MAJOR);
double phi = M_PI_2 - 2 * atan(ts);
double dphi = 1.0;
int i;
for (i = 0; fabs(dphi) > 0.000000001 && i < 15; i++) {
double con = ECCENT * sin (phi);
dphi = M_PI_2 - 2 * atan (ts * pow((1.0 - con) / (1.0 + con), COM)) - phi;
phi += dphi;
}
}```

To compile in Visual Studio / MS Windows OS, I had to add these definitions MikeCollinson 14:17, 20 January 2007 (UTC):

```// Add this line before including math.h:
#define _USE_MATH_DEFINES
// Additions for MS Windows compilation:
#ifndef M_PI
#define M_PI acos(-1.0)
#endif
#ifndef M_PI_2
#define M_PI_2 1.57079632679489661922
#endif
inline double fmin(double x, double y) { return(x < y ? x : y); }
inline double fmax(double x, double y) { return(x > y ? x : y); }```

### C#

C# Implementation by Florian Müller, based on the C code published above, 14:50, 20.6.2008; updated to static functions by David Schmitt, 23.4.2010

```using System;

public static class MercatorProjection
{
private static readonly double R_MAJOR = 6378137.0;
private static readonly double R_MINOR = 6356752.3142;
private static readonly double RATIO = R_MINOR / R_MAJOR;
private static readonly double ECCENT = Math.Sqrt(1.0 - (RATIO * RATIO));
private static readonly double COM = 0.5 * ECCENT;

private static readonly double DEG2RAD = Math.PI / 180.0;
private static readonly double RAD2Deg = 180.0 / Math.PI;
private static readonly double PI_2 = Math.PI / 2.0;

public static double[] toPixel(double lon, double lat)
{
return new double[] { lonToX(lon), latToY(lat) };
}

public static double[] toGeoCoord(double x, double y)
{
return new double[] { xToLon(x), yToLat(y) };
}

public static double lonToX(double lon)
{
return R_MAJOR * DegToRad(lon);
}

public static double latToY(double lat)
{
lat = Math.Min(89.5, Math.Max(lat, -89.5));
double phi = DegToRad(lat);
double sinphi = Math.Sin(phi);
double con = ECCENT * sinphi;
con = Math.Pow(((1.0 - con) / (1.0 + con)), COM);
double ts = Math.Tan(0.5 * ((Math.PI * 0.5) - phi)) / con;
return 0 - R_MAJOR * Math.Log(ts);
}

public static double xToLon(double x)
{
return RadToDeg(x) / R_MAJOR;
}

public static double yToLat(double y)
{
double ts = Math.Exp(-y / R_MAJOR);
double phi = PI_2 - 2 * Math.Atan(ts);
double dphi = 1.0;
int i = 0;
while ((Math.Abs(dphi) > 0.000000001) && (i < 15))
{
double con = ECCENT * Math.Sin(phi);
dphi = PI_2 - 2 * Math.Atan(ts * Math.Pow((1.0 - con) / (1.0 + con), COM)) - phi;
phi += dphi;
i++;
}
}

{
}

private static double DegToRad(double deg)
{
return deg * DEG2RAD;
}
}```

### PHP Code

Php Code by Erhan Baris 19:19, 01.09.2007

```function merc_x(\$lon)
{
\$r_major = 6378137.000;
return \$r_major * deg2rad(\$lon);
}

function merc_y(\$lat)
{
if (\$lat > 89.5) \$lat = 89.5;
if (\$lat < -89.5) \$lat = -89.5;
\$r_major = 6378137.000;
\$r_minor = 6356752.3142;
\$temp = \$r_minor / \$r_major;
\$es = 1.0 - (\$temp * \$temp);
\$eccent = sqrt(\$es);
\$sinphi = sin(\$phi);
\$con = \$eccent * \$sinphi;
\$com = 0.5 * \$eccent;
\$con = pow((1.0-\$con)/(1.0+\$con), \$com);
\$ts = tan(0.5 * ((M_PI*0.5) - \$phi))/\$con;
\$y = - \$r_major * log(\$ts);
return \$y;
}

function merc(\$x,\$y) {
return array('x'=>merc_x(\$x),'y'=>merc_y(\$y));
}

\$array = merc(122,11);```

### Java Implementation

Java Implementation by Moshe Sayag, based on the JavaScript code published above, 17:11, 15.1.2008

```public class Mercator {
final private static double R_MAJOR = 6378137.0;
final private static double R_MINOR = 6356752.3142;

public double[] merc(double x, double y) {
return new double[] {mercX(x), mercY(y)};
}

private double  mercX(double lon) {
return R_MAJOR * Math.toRadians(lon);
}

private double mercY(double lat) {
if (lat > 89.5) {
lat = 89.5;
}
if (lat < -89.5) {
lat = -89.5;
}
double temp = R_MINOR / R_MAJOR;
double es = 1.0 - (temp * temp);
double eccent = Math.sqrt(es);
double phi = Math.toRadians(lat);
double sinphi = Math.sin(phi);
double con = eccent * sinphi;
double com = 0.5 * eccent;
con = Math.pow(((1.0-con)/(1.0+con)), com);
double ts = Math.tan(0.5 * ((Math.PI*0.5) - phi))/con;
double y = 0 - R_MAJOR * Math.log(ts);
return y;
}
}```

### Python Implementation

Python Implementation by Paulo Silva, based on all code published above, 13:32, 15.2.2008

```import math

def merc_x(lon):
r_major=6378137.000

def merc_y(lat):
if lat>89.5:lat=89.5
if lat<-89.5:lat=-89.5
r_major=6378137.000
r_minor=6356752.3142
temp=r_minor/r_major
eccent=math.sqrt(1-temp**2)
sinphi=math.sin(phi)
con=eccent*sinphi
com=eccent/2
con=((1.0-con)/(1.0+con))**com
ts=math.tan((math.pi/2-phi)/2)/con
y=0-r_major*math.log(ts)
return y```

### sdlBasic Implementation

sdlBasic Implementation by Paulo Silva, based on all code published above, 12:33, 18.2.2008

```function deg_rad(ang):
tang=ang
pi=3.14159265359
end function

function merc_x(lon):
tlon=lon
r_major=6378137.0
end function

function merc_y(lat):
tlat=lat
pi=3.14159265359
if tlat>89.5 then:tlat=89.5: end if
if tlat<-89.5 then:tlat=-89.5: end if
r_major=6378137.000
r_minor=6356752.3142
temp=r_minor/r_major
es=1-(temp*temp)
eccent=sqr(es)
phi=(tlat*pi)/180
sinphi=sin(phi)
con=eccent*sinphi
com=.5*eccent
con=((1.0-con)/(1.0+con))^com
ts=tan(.5*((pi*0.5)-phi))/con
y=0-r_major*log(ts)
merc_y=y
end function```

### bash-script (bc)

bash-script by Frank Baettermann, 18.05.2008

```#!/bin/bash

# Compute mercartor tile-coordinates

if [ \$# -ne 3 ]; then
echo "Usage: \$0 LATITUDE LONGITUDE ZOOM"
exit 1
fi

if [ -z "`which bc 2> /dev/null`" ]; then
echo "ERROR: Could not find bc!"
exit 1
fi

LATITUDE=\$1
LONGITUDE=\$2
ZOOM=\$3

# total width and height in tiles
MAP_SIZE=\$((2**\$ZOOM))
echo "map_size=\$MAP_SIZE"

# longitude -> x
echo "tile_x=`echo "(\$LONGITUDE + 180) * \$MAP_SIZE / 360" | bc`"

# latitude -> y
echo "tile_y=`echo "\
scale=10; \
phi=\$LATITUDE; \
mapsize=\$MAP_SIZE; \
pi=3.1415926535; \
tile=(1-mercator/pi)*(mapsize/2); \
scale=0; tile/1" \
| bc -l`"```

### VBA for use in Excel

Copy into a module to use as an Excel Function

Based on JavaScript Code above by Benjamin Judd

```Function deg_rad(deg As Double) As Double
deg_rad = deg * Excel.WorksheetFunction.Pi / 180
End Function

Function merc_x(lon As Double) As Double
Dim r_major As Double
r_major = 6378137
merc_x = r_major * deg_rad(lon)
End Function

Function merc_y(lat As Double) As Double
If (lat > 89.5) Then lat = 89.5
If (lat < -89.5) Then lat = -89.5
Dim r_major As Double, r_minor As Double, temp As Double, es As Double, eccent As Double
Dim phi As Double, sinphi As Double, con As Double, com As Double, ts As Double
r_major = 6378137
r_minor = 6356752.3142
temp = r_minor / r_major
es = 1 - temp ^ 2
eccent = Math.Sqr(es)