Mercator

From OpenStreetMap Wiki
Jump to navigation Jump to search
Mercator - Other languages
For the OSM editor program, see Merkaartor.

This article describes algorithms for Mercator projection.

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. Combined with the zoom level, the system is known as a Web Mercator on Wikipedia.

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)

The following functions take or return their results in degrees:

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 lon2x(lon) { return lon; }

The following take or return their results in Meters and Lat/Lon (similar to Python implementation below):

PI = Math.PI;
RAD2DEG = 180/PI;
DEG2RAD = PI/180;
R =  6378137.0

function y2lat(y) { return (2 * Math.atan( Math.exp(y / R)) - PI/2 ) * RAD2DEG }
function x2lon(x) { return RAD2DEG * (x/R); }

function lat2y(lat) { return Math.log(Math.tan(PI/4 + lat * DEG2RAD / 2)) * R}
function lon2x(lon) { return lon * DEG2RAD * R; }

C

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

/* 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); }
double x2lon_m(double x)   { return RAD2DEG(              x/EARTH_RADIUS           ); }

/* 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; }
double lon2x_m(double lon) { return          DEG2RAD(lon)                 * 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;
    }

D

double lon2x(double lon) { import std.math; return PI*(lon/180.0) * 6378137.0; }
double lat2y(double lat) { import std.math; return log(tan(PI_4+PI*(lat/360.0)))*6378137.0; }
double x2lon(double x) { import std.math; return 180*(x/6378137)/PI; }
double y2lat(double y) { import std.math; return 180*(2*atan(exp(y/6378137))-PI_2)/PI; }

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
A y lat: = ARCTAN(EXP(A1/180*PI())/PI()*360-90
B 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


Excel formulas for converting Lat/Lon to Meters

Lat to Y (expecting Lat in A1):

=LN(TAN((PI()/4)+(RADIANS(A1)/2)))*6378137

Lon to X (expecting Lon in B1):

=RADIANS(B1)*6378137

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) {
    return Math.toDegrees(aX / RADIUS);
  }

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

  public static double lat2y(double aLat) {
    return Math.log(Math.tan(Math.PI / 4 + Math.toRadians(aLat) / 2)) * RADIUS;
  }  
  public static double lon2x(double aLong) {
    return Math.toRadians(aLong) * RADIUS;
  }
}

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:

[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');

Python

import math
R =  6378137.0
def y2lat(y):
    return math.degrees(2 * math.atan(math.exp (y / R)) - math.pi / 2.0)
def lat2y(lat):
    return math.log(math.tan(math.pi / 4 + math.radians(lat) / 2)) * R
def x2lng(x):
    return math.degrees(x / R)
def lon2x(lon):
    return math.radians(lon) * R

Rust

// length of semi-major axis of the WGS84 ellipsoid, i.e. radius at equator
const EARTH_RADIUS_KM: f64 = 6_378.137;
    
pub fn lon2x(lon: f64) -> f64 {
    EARTH_RADIUS_KM * 1000. * lon.to_radians()
}

pub fn x2lon(x: f64) -> f64 {
    (x / (EARTH_RADIUS_KM * 1000.)).to_degrees()
}

pub fn lat2y(lat: f64) -> f64 {
    ((lat.to_radians() / 2. + std::f64::consts::PI / 4.).tan()).log(std::f64::consts::E)
        * EARTH_RADIUS_KM
        * 1000.
}

pub fn y2lat(y: f64) -> f64 {
    (2. * ((y / (EARTH_RADIUS_KM * 1000.)).exp()).atan() - std::f64::consts::PI / 2.)
        .to_degrees()
}

Swift

class SphericalMercator {
    private let radius: Double = 6378137.0; /* in meters on the equator */
    
    /* These functions take their length parameter in meters and return an angle in degrees */
    
    func y2lat(aY: Double) -> Double {
        (atan(exp(aY / radius)) * 2 - Double.pi/2).radiansToDegrees
    }
    
    func x2lon(aX: Double) -> Double {
        (aX / radius).radiansToDegrees;
    }
    
    /* These functions take their angle parameter in degrees and return a length in meters */
    
    func lat2y(aLat: Double) -> Double {
        log(tan(Double.pi / 4 + (aLat.degreesToRadians / 2))) * radius
    }
    
    func lon2x(aLong: Double) -> Double {
        (aLong).degreesToRadians * radius;
    }
}

extension FloatingPoint {
    var degreesToRadians: Self { self * .pi / 180 }
    var radiansToDegrees: Self { self * 180 / .pi }
}

Go (Golang)

package main

import "math"

const R = 6378137.0

func Radians(deg float64) float64 {
	return deg * math.Pi / 180
}

func Degrees(rad float64) float64 {
	return rad * 180 / math.Pi
}

func Y2lat(y float64) float64 {
	return Degrees(2*math.Atan(math.Exp(y/R)) - math.Pi/2)
}

func Lat2y(lat float64) float64 {
	return R * math.Log(math.Tan(math.Pi/4+Radians(lat)/2))
}

func X2lat(x float64) float64 {
	return Degrees(x / R)
}

func Lon2x(lon float64) float64 {
	return R * Radians(lon)
}

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 projection is not commonly 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;
    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# 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 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;
    }
}

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

PHP implementation

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

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

Retrieved from "https://wiki.openstreetmap.org/w/index.php?title=Mercator&oldid=2745139"