# Mercator

For the OSM editor program, see Merkaartor.

## Spherical Pseudo-Mercator projection

Most of OSM, including the main tiling system, uses a Pseudo-Mercator projection where the Earth is modelized as if it was a perfect a sphere.

This produces a fast approximation to the truer, but heavier elliptical projection, where the Earth would be projected on a more accurate ellipsoid (flattened on poles). As a consequence, direct mesurements of distances in this projection will be approximative, except on the Equator, and the aspect ratios on the rendered map for true squares measured on the surface on Earth will slightly change with latitude and angles not so precisely preserved by this spherical projection.

### JavaScript (or ActionScript)

```RAD2DEG = 180 / Math.PI;
PI_4 = Math.PI / 4;

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

function y2lat(y) { return (Math.atan(Math.exp(y / RAD2DEG)) / PI_4 - 1) * 90; }
function x2lon(x) { return x; }

function lat2y(lat) { return Math.log(Math.tan((lat / 90 + 1) * PI_4 )) * RAD2DEG; }
function lon2y(lon) { return lon; }
```

### C

```#include <math.h>
#define DEG2RAD(a)   ((a) / (180 / M_PI))
#define RAD2DEG(a)   ((a) * (180 / M_PI))

/* The following functions take their parameter and return their result in degrees */

double y2lat_d(double y)   { return RAD2DEG( atan(exp( DEG2RAD(y) )) * 2 - M_PI/2 ); }
double x2lon_d(double x)   { return x; }

double lat2y_d(double lat) { return RAD2DEG( log(tan( DEG2RAD(lat) / 2 +  M_PI/4 )) ); }
double lon2x_d(double lon) { return lon; }

/* The following functions take their parameter in something close to meters, along the equator, and return their result in degrees */

double y2lat_m(double y)   { return RAD2DEG(2 * atan(exp( y/EARTH_RADIUS)) - M_PI/2); }

/* The following functions take their parameter in degrees, and return their result in something close to meters, along the equator */

double lat2y_m(double lat) { return log(tan( DEG2RAD(lat) / 2 + M_PI/4 )) * EARTH_RADIUS; }
```

### C#

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

### Excel

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

 1 2 y lat: = ARCTAN(EXP(A1/180*PI())/PI()*360-90 y: = LOG(TAN((B2+90)/360*PI()))/PI()*180 lat

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

### Java

```import java.lang.Math;
public class SphericalMercator {
public static final double RADIUS = 6378137.0; /* in meters on the equator */

/* These functions take their length parameter in meters and return an angle in degrees */

public static double y2lat(double aY) {
return Math.toDegrees(Math.atan(Math.exp(aY / RADIUS)) * 2 - Math.PI/2);
}
public static double x2lon(double aX) {
}

/* These functions take their angle parameter in degrees and return a length in meters */

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

### 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
```

### Perl

```use Geo::Mercator;
my (\$x, \$y) = mercate(\$lat, \$long)
```

### 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); }
```

### 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');
```

### 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))
```

FIXME: add spherical code in other languages here

## Elliptical (true) Mercator Projection

This projection gives more accurate aspect ratios for objects anywhere on Earth, and respects their angles with higher precision. However this project is not used on most maps used on the OSM websites and in editors.

### JavaScript (or ActionScript) implementation

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 (elliptical) Mercator projection, in Javascript:

```function deg_rad(ang) {
return ang * (Math.PI/180.0)
}
function merc_x(lon) {
var r_major = 6378137.000;
}
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 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
{

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 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 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) {
}

double merc_y (double lat) {
lat = fmin (89.5, fmax (lat, -89.5));
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) {
}

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# implementation

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 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)
{
}

public static double latToY(double lat)
{
lat = Math.Min(89.5, Math.Max(lat, -89.5));
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)
{
}

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++;
}
}

{
}

{
}
}
```

### 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) {
}

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 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;
}
}
```

### PHP implementation

Php Code by Erhan Baris 19:19, 01.09.2007

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

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);
```

### 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) implementation

bash-script by Frank Baettermann, 18.05.2008

```#!/bin/bash

# Compute Mercator 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! Try to install it: sudo apt install 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 Office implementation

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
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)