Jingwen Zheng
  • Jan 10, 2021

How to convert projected coordinates to latitude longitude by Python?

This blog talks about how to convert projected coordinates to latitude longitude with python module "pyproj".

Last month I drew subway and tramway lines with folium according to the Île-de-France public transport network datasets. However, the coordinates are not what we use in GPS (EPSG:4326), but the Mercator projection (EPSG:3857).

input1

…

input2

Thus, before drawing the transport lines, we need to convert projected coordinates to latitude longitude. In this blog, I’ll talk about this with the following points:

  • What is EPSG:4326?
  • What is EPSG:3857?
  • How to convert projected coordinates to latitude longitude?

What is EPSG:4326?

The World Geodetic System (WGS) is a standard for use in cartography, geodesy, and satellite navigation including GPS. This standard includes the definition of the coordinate system’s fundamental and derived constants, the ellipsoidal (normal) Earth Gravitational Model (EGM), a description of the associated World Magnetic Model (WMM), and a current list of local datum transformations. The latest revision is WGS 84 (also known as WGS 1984, EPSG:4326), established and maintained by the United States National Geospatial-Intelligence Agency since 1984, and last revised in 2004. Earlier schemes included WGS 72, WGS 66, and WGS 60. WGS 84 is the reference coordinate system used by the Global Positioning System.

What is EPSG:3857?

Web Mercator, Google Web Mercator, Spherical Mercator, WGS 84 Web Mercator or WGS 84/Pseudo-Mercator is a variant of the Mercator projection and is the de facto standard for Web mapping applications. It rose to prominence when Google Maps adopted it in 2005. It is used by virtually all major online map providers, including Google Maps, Mapbox, Bing Maps, OpenStreetMap, Mapquest, Esri, and many others. Its official EPSG identifier is EPSG:3857, although others have been used historically.

How to convert projected coordinates to latitude longitude?

The Python module pyproj is a cartographic projections and coordinate transformations library, which converts from longitude, latitude to native map projection x,y coordinates and vice versa, its classes Proj and transform will help us to get the GPS coordinates.

A Proj class instance is initialized with proj map projection control parameter key/value pairs. The key/value pairs can either be passed in a dictionary, or as keyword arguments, or as a PROJ string (compatible with the proj command). performs cartographic transformations: converting from longitude, latitude to native map projection x,y coordinates and vice versa. The Transform.er class is also for facilitating re-using transforms without needing to re-create them. The goal is to make repeated transforms faster.

Method 1

import pyproj

wgs84 = pyproj.Proj(projparams = 'epsg:4326')
InputGrid = pyproj.Proj(projparams = 'epsg:3857')

x1, y1 = -11705274.6374,4826473.6922
pyproj.transform(InputGrid, wgs84, x1, y1)

We firstly set 2 variables wgs84 and InputGrid by specifying the projparams, then transfer the point (x1, y1) from the coordinate system “EPSG:3857” to the “EPSG:4326”.

method1

We got the results but with a warning message as above. After checking the recommendation, now we have another method.

Method 2

transformer = pyproj.Transformer.from_crs("epsg:3857", "epsg:4326")
transformer.transform(x1, y1)

For the second method, we create a transformer from a pyproj.crs.CRS, with crs_from="epsg:3857" and crs_to="epsg:4326", then transform the point (x1, y1), we can get the same result as the first method, without any warning message.

method2

Use case

So how to solve the problem at the beginning of this article? Idea is that we go through the projected coordinates, and transform each of them to “EPSG:4326”.

transformer = pyproj.Transformer.from_crs("epsg:3857", "epsg:4326")
for feature in subway_lines_df['features']:
    if feature['geometry']['type'] == 'LineString':
        # all coordinates
        coords = feature['geometry']['coordinates']
        # coordList is for each individual polygon
        for coordList in coords:
            lat_grid, lng_grid = coordList
            # do transformation
            # coordList[0],coordList[1] = pyproj.transform(InputGrid, wgs84, lat_grid, lng_grid) # methode 1
            coordList[0],coordList[1] = transformer.transform(lat_grid, lng_grid) # methode 2
    elif feature['geometry']['type'] == 'MultiLineString':
        len_coord = len(feature['geometry']['coordinates'])
        for c in range(len_coord):
            coords = feature['geometry']['coordinates'][c]
            for coordList in coords:
                lat_grid, lng_grid = coordList
                # coordList[0],coordList[1] = pyproj.transform(InputGrid, wgs84, lat_grid, lng_grid) # methode 1
                coordList[0],coordList[1] = transformer.transform(lat_grid, lng_grid) # methode 2

If you are curious about the scripts, you will find them here.

References

Search