Mercator

From OpenStreetMap Wiki
Jump to: navigation, search

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>
#define deg2rad(d) (((d)*M_PI)/180)
#define rad2deg(d) (((d)*180)/M_PI)
#define earth_radius 6378137
 
/* 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:

[900913.sql]

 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]],
 PROJECTION["Mercator_1SP_Google"],PARAMETER["latitude_of_origin", 0.0],
 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) {
    return Math.toDegrees(Math.log(Math.tan(Math.PI/4+Math.toRadians(aLat)/2)));
  }
}

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_major:6378137.0,//Equatorial Radius, WGS84
	r_minor:6356752.314245179,//defined as constant
	f:298.257223563,//1/f=(a-b)/a , a=r_major, b=r_minor
	deg2rad:function(d)
	{
		var r=d*(Math.PI/180.0);
		return r;
	},
	rad2deg:function(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 lon=this.rad2deg((x/this.r_major));
 
		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;
        }
        return rad_deg (phi);
}

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++;
        }
        return RadToDeg(phi);
    }
 
    private static double RadToDeg(double rad)
    {
        return rad * RAD2Deg;
    }
 
    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);
    $phi = deg2rad($lat);
    $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
  return r_major*math.radians(lon)
 
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)
  phi=math.radians(lat)
  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
  deg_rad=tang*(pi/180.0)
  end function
 
function merc_x(lon):
  tlon=lon
  r_major=6378137.0
  merc_x=r_major*deg_rad(tlon)
  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; \
    phi_rad=phi*2*pi/360; \
    mercator=l((1+s(phi_rad))/(1-s(phi_rad)))/(2); \
    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)
 phi = deg_rad(lat)
 sinphi = Math.Sin(phi)
 con = eccent * sinphi
 com = 0.5 * eccent
 con = (1 - con / 1 + con) ^ com
 ts = Math.Tan(0.5 * ((Excel.WorksheetFunction.Pi * 0.5) - phi)) / con
 merc_y = 0 - r_major * Math.Log(ts)
End Function

Links

Lat/Lon to Mercator online converter