I've created a couple of functions using the math package from python 3.5 to transform GPS coordinates to x/y web mercator (pseudo mercator) as explained here: http://wiki.openstreetmap.org/wiki/Mercator#Java

I took the Java version and ported it to Python:

import math

# derived from the Java version explained here: http://wiki.openstreetmap.org/wiki/Mercator
RADIUS = 6378137.0 # in meters on the equator

def lat2y(a):
  return math.log(math.tan(math.pi / 4 + math.radians(a) / 2)) * RADIUS

def lon2x(a):
  return math.radians(a) * RADIUS

Then I can use them to get for example the projection of Trafalgar Square (taking the GPS coordinates from the Google Maps URL):

ts_gm = [51.50809,-0.1285907]
print('latitude web mercator y: {} longitude web mercator x: {}'.format(lat2y(ts_gm[0]), lon2x(ts_gm[1])))

Which gives back:

latitude web mercator y: 6711665.883938471 longitude web mercator x: -14314.651244750603

However I would like to achieve this using pyproj syntax, but I am lost in the details of the string I should pass to build a projection using Web Mercator EPSG:3857. I suspect the formulas above have some sort of precision issue and I would like to rely on a library like pyproj hoping to minimise projection errors.

Is there a minimum working example I could see to achieve the mapping in python using pyproj?

from pyproj import Proj, transform
print(transform(Proj(init='epsg:4326'), Proj(init='epsg:3857'), -0.1285907, 51.50809))  # longitude first, latitude second.
# output (meters east of 0, meters north of 0): (-14314.651244750548, 6711665.883938471)

The "trick" is to use these shortcuts for Web Mercator (EPSG 3857) and WGS 84 longitude and latitude (EPSG 4326). The init key means "initialize this projection by reading the definition for 3857 from the 'epsg' file."

  • Thanks! Is the zoom level influencing at all the precision of this transformation? For example in here bokeh.pydata.org/en/latest/docs/user_guide/… they are using a bounding box with side 2 * 20,000,000 meters for both x and y and I was wondering if the above transformation is the most precise that someone could use assuming the GPS/WGS84 coordinates are "in Europe" – TPPZ Jul 17 '17 at 13:53
  • I double checked a few pairs of coordinates, this transformation does not seem to depend on the scale / zoom level. The formulas I mentioned do not contain a hook to the zoom level, but somewhere I read the pseudo mercator renders better in some areas and for some reason I thought this was not because of the distortion due to latitude. – TPPZ Jul 17 '17 at 22:25

With the latest version of pyproj, the recommended way is

from pyproj import Transformer

TRAN_4326_TO_3857 = Transformer.from_crs("EPSG:4326", "EPSG:3857")

def mytransform(lon, lat):
  return TRAN_4326_TO_3857.transform(lon, lat)

This can be orders of magnitude faster, especially if you need to convert many coordinates

See: https://pyproj4.github.io/pyproj/dev/advanced_examples.html#advanced-examples

  • 1
    it should be TRAN_4326_TO_3857.transform(lon, lat) – fccoelho Dec 15 '20 at 12:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.