# Calculating Mercator coordinates from lat/lon

I am new to GIS but I am working on a project mapping a building using floor plans and JOSM.

Using the PicLayer plugin I added the .png floor plan as a layer.

I wrote a .cal file to set the scale, rotation, and center position of the image.

In the JOSM default Mercator projection, the image moved only slightly when the center was changed (using lat and lon decimal degrees). The image was still centered somewhere around zero. In WGS84 the image moved to the correct location but was sheared and stretched.

I found this reference on converting lat/lon to a Mercator projection: http://wiki.openstreetmap.org/wiki/Mercator

The Python code under Elliptical Mercator correctly converts my longitude, but the latitude is off by about 0.5% (a few metres).

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

My nodes and ways were lined up correctly with the image with this latitude and longitude in WGS84, so I think they are correct and the conversion is not. I am in the northern hemisphere and using an area near 49 N, -122 E.

• So I was looking at the elliptical Mercator code because the Python for spherical Mercator didn't show X and gave me the wrong Y value. In the ActionScript section (wiki.openstreetmap.org/wiki/…), it says that the X coordinate is the longitude multiplied by a chosen scale. I guessed that I could take the output X value from the elliptical Mercator conversion and divide it by the actual longitude to get the scale, and it worked. I then took this value and multiplied it by the result of lat2y to get the Mercator Y value. Will post as answer. Jul 28, 2015 at 6:12

## 3 Answers

Take a look at Proj4 https://github.com/OSGeo/proj.4/wiki python's implementation pyproj https://github.com/jswhit/pyproj

You could specify src and dest projections (like `EPSG:4326` and `EPSG:3857`) and convert coordinates between them. And pyproj will made all raw calculations.

Note: this solution targets array inputs. Inspired by the accepted @JoshVazquez's solution.

### If `lat` and `lon` are both arrays (or dataframe columns)...

The current accepted answer would work very well when `lat` and `lon` are scalar parameters.

I've however recently bumped into a scenario where `lat` and `lon` are numpy arrays (or pandas dataframe columns). To make this work for arrays / dataframe columns, I've essentially modified @JoshVazquez's solution above - essentially replacing all the `math` functions with `numpy` functions.

``````import numpy as np
def merc_from_arrays(lats, lons):
r_major = 6378137.000
x = r_major * np.radians(lons)
scale = x/lons
y = 180.0/np.pi * np.log(np.tan(np.pi/4.0 + lats * (np.pi/180.0)/2.0)) * scale
return (x, y)

lats = np.asarray([54.984105, 56])
lons = np.asarray([-3.193693, -2.2])

xs, ys = merc_from_arrays(lats, lons)

print('xs: {}'.format(xs))  # => [-355520.27851004 -244902.8797452 ]
print('ys: {}'.format(ys))  # => [ 7358781.82857011  7558415.65608178]
``````
• There is a bug in this code: in the mere_from_arrays function, you reference variables lon and lat, but those do not exist.Presumably you meant lats and lons, the arguments to the function. Feb 7, 2019 at 0:20
• thank you @ibrewster - typo corrected now. Feb 9, 2019 at 10:35

Longitude conversion from elliptical Mercator (same formula in both systems?:

``````import math
def merc_x(lon):
r_major=6378137.000
return r_major*math.radians(lon)
``````

Mercator X now known. Divide this by the longitude to get the scale value.

Latitude conversion from spherical Mercator:

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

Multiply the result by the scale value to get Mercator Y value.

In one function:

``````import math
def merc(lat, lon):
r_major = 6378137.000
x = r_major * math.radians(lon)
scale = x/lon
y = 180.0/math.pi * math.log(math.tan(math.pi/4.0 + lat * (math.pi/180.0)/2.0)) * scale
return (x, y)
``````
• I converted this to js, gist.github.com/anonymous/3c82d80547d3ed26baf5d87f85855198 , however, i think im getting wrong answer. I thought top left would be 0,0 `merc(90, -180) == 0,0` but im getting `238107693.26496768` Mar 19, 2018 at 4:28
• There is an error with the formula, the 'scale' variable will throw a division error if the lon is 0. It should be replaced by scale = r_major * π / 180 Sep 3, 2020 at 5:40
• @MuhammadUmer the projection cannot be used near the poles, anywhere above 85 lat or below -85 lat. See developers.google.com/maps/documentation/javascript/… Sep 3, 2020 at 5:42