Mercator
From OpenStreetMap
This article describes algorithms for Mercator projection.
Contents |
Spherical Mercator
Most of OSM, including the main tiling system, uses a spherical Mercator projection. This produces a fast approximation to the truer ellptical 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))
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)];
}
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 = 0;
while ((fabs(dphi) > 0.000000001) && (i < 15)) {
double con = ECCENT * sin (phi);
dphi = M_PI_2 - 2 * atan (ts * pow((1.0 - con) / (1.0 + con), COM)) - phi;
phi += dphi;
i++;
}
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):
// 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
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 double[] toPixel(double lon, double lat){
return new double[]{this.lonToX(lon), this.latToY(lat)};
}
public double[] toGeoCoord(double x, double y){
return new double[] { this.XToLon(x), this.YToLat(y) };
}
public double lonToX(double lon){
return R_MAJOR * this.DegToRad(lon);
}
public double latToY(double lat){
lat = Math.Min(89.5, Math.Max(lat, -89.5));
double phi = this.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 double xToLon(double x){
return this.RadToDeg(x) / R_MAJOR;
}
public 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 this.RadToDeg(phi);
}
private double RadToDeg(double rad){
return rad * RAD2Deg;
}
private double DegToRad(double deg){
return deg * DEG2RAD;
}
}
PHP Code
Php Code by Erhan Baris 19:19, 01.09.2007
function deg_rad($ang)
{
return (float)((float)$ang * (float)(M_PI / 180.0));
}
function merc_x($lon)
{
$r_major = 6378137.000;
return (float)($r_major * deg_rad($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 = deg_rad($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 = 0 - $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 deg_rad(ang): return ang*(math.pi/180.0) def merc_x(lon): r_major=6378137.000 return r_major*deg_rad(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 es=1-(temp*temp) eccent=math.sqrt(es) phi=(lat*math.pi)/180 sinphi=math.sin(phi) con=eccent*sinphi com=.5*eccent con=math.pow(((1.0-con)/(1.0+con)),com) ts=math.tan(.5*((math.pi*0.5)-phi))/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`"
