9

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?

Improve this question

3 Answers 3

Reset to default
16 
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."

Improve this answer
3
  • 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" Commented Jul 17, 2017 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. Commented Jul 17, 2017 at 22:25
  • Note that the recommended way of doing this with newer version of pyproj is shown in the answer from Jacopofar. Commented Mar 9, 2022 at 19:06
5

With the latest version of pyproj, the recommended way is

from pyproj import Transformer # See other answer about the always_xy=True parameter TRAN_4326_TO_3857 = Transformer.from_crs("EPSG:4326", "EPSG:3857") def mytransform(lat, lon): return TRAN_4326_TO_3857.transform(lat, lon)

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

Improve this answer
3
  • 1
    it should be TRAN_4326_TO_3857.transform(lon, lat) Commented Dec 15, 2020 at 12:10
  • 1
    I'm by no means an expert on geo systems, but in the experiments that I'm personally running the reversal of (lat, lon) that was applied here in a recent edit gives me the wrong results whereas the original gives me the correct results. Commented Nov 29, 2021 at 8:18
  • Since I was also confused about the order of x, y or lon, lat I added the relevant info to a separate answer. Commented Mar 9, 2022 at 19:04
3

Beware! Order of x, y or lon, lat might not be intuitive

Since there is some confusion about the order of x, y or lon, lat in pyproj I want to summarize the info that I just found after also messing up the order when reusing the solution from Jacopofar's answer.

From the FAQ - Why is the axis ordering in PROJ not consistent?

PROJ respects the axis ordering as it was defined by the authority in charge of a given coordinate reference system. This is in accordance to the ISO19111 standard

Warning from the doc of Transformer

The axis order may be swapped if the source and destination CRS’s are defined as having the first coordinate component point in a northerly direction (See PROJ FAQ on axis order). You can check the axis order with the pyproj.crs.CRS class. If you prefer to keep your axis order as always x,y, you can use the always_xy option when creating the pyproj.transformer.Transformer.

So, to make things clear in the code, my suggested variation of Jacopofar's solution using always_xy=True is

 from pyproj import Transformer lonlat_to_webmercator = Transformer.from_crs("EPSG:4326", "EPSG:3857", always_xy=True) def mytransform(lon, lat): x, y = lonlat_to_webmercator.transform(lon, lat) return x, y # lon-lat of center of Munich (Marienplatz) mytransform(11.576124, 48.137154)

correct result (which I verified on a webmercator map)

(1288648.2290397931, 6129702.78025012)
Improve this answer

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.