Calculating Mercator coordinates from lat/lon

Question

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.

Answer 1

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.

Answer 2

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]

Answer 3

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)