7. Working with CRS and Projections in Python#

Download this compressed folder to get the notebook for this chapter and the corresponding environment.yml file with necessary packages to install.

The goal of this lecture is to work with PyProj, Cartopy and Shapely packages in Python and explore different crs formats. We will also transform data from one projection to another one.

Attribution: Parts of this notebook is developed using great resources from CLEX CMS Blog and Projection Tradeoffs and Distortion - Tissot

import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import pyproj
import shapely

Cartopy raises warning anytime you convert CRS from one format to another. Particulalrly, something like this: UserWarning: You will likely lose important projection information when converting to a PROJ string from another format. The following cell turns off all warnings to keep your notebook clean.

import warnings
warnings.filterwarnings('ignore')

7.1. Convert CRS from one format to the other#

You can generate a crs using caropy’s crs module. For example:

[ccrs.Mercator()]

[<Projected CRS: +proj=merc +ellps=WGS84 +lon_0=0.0 +x_0=0.0 +y_0=0 ...>
 Name: unknown
 Axis Info [cartesian]:
 - E[east]: Easting (metre)
 - N[north]: Northing (metre)
 Area of Use:
 - undefined
 Coordinate Operation:
 - name: unknown
 - method: Mercator (variant A)
 Datum: Unknown based on WGS 84 ellipsoid
 - Ellipsoid: WGS 84
 - Prime Meridian: Greenwich]

Similarly, you can use pyproj to generate a crs

pyproj.CRS(4326)

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

You can use built in function in both pyproj and cartopy to conver the crs between different formats

pyproj.CRS(ccrs.Mercator().to_proj4()).to_epsg()

3395

pyproj.CRS(ccrs.Mercator().to_proj4()).to_wkt()

'PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unknown based on WGS 84 ellipsoid",ELLIPSOID["WGS 84",6378137,298.257223563,LENGTHUNIT["metre",1],ID["EPSG",7030]]],PRIMEM["Greenwich",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8901]]],CONVERSION["unknown",METHOD["Mercator (variant A)",ID["EPSG",9804]],PARAMETER["Latitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8801]],PARAMETER["Longitude of natural origin",0,ANGLEUNIT["degree",0.0174532925199433],ID["EPSG",8802]],PARAMETER["Scale factor at natural origin",1,SCALEUNIT["unity",1],ID["EPSG",8805]],PARAMETER["False easting",0,LENGTHUNIT["metre",1],ID["EPSG",8806]],PARAMETER["False northing",0,LENGTHUNIT["metre",1],ID["EPSG",8807]]],CS[Cartesian,2],AXIS["(E)",east,ORDER[1],LENGTHUNIT["metre",1,ID["EPSG",9001]]],AXIS["(N)",north,ORDER[2],LENGTHUNIT["metre",1,ID["EPSG",9001]]]]'

pyproj.CRS("World_Mercator")

<Projected CRS: ESRI:54004>
Name: World_Mercator
Axis Info [cartesian]:
- E[east]: Easting (metre)
- N[north]: Northing (metre)
Area of Use:
- name: World.
- bounds: (-180.0, -90.0, 180.0, 90.0)
Coordinate Operation:
- name: World_Mercator
- method: Mercator (variant B)
Datum: World Geodetic System 1984 ensemble
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

7.2. Tissot’s indicatrix#

Next, we will plot Tissot’s indicatrix for several projections and study map distortions:

#Function to create cartopy plot with Tissot circles for input crs
def tissotplot(crs):
    
    #Create figure and axes
    fig, ax = plt.subplots(subplot_kw={'projection': crs}, figsize=(8,8))
    
    #Define positions of the circles, spaced every 20 degrees
    lons = range(-180, 180, 20)
    lats = range(-80, 81, 20)
    
    #Draw coastlines
    ax.coastlines()
    
    #Add gridlines
    gl = ax.gridlines(draw_labels=True, auto_inline=True, alpha=0.5, lw=0.5, linestyle=':')
    gl.xlocator = mticker.FixedLocator(lons)
    gl.ylocator = mticker.FixedLocator(lats)
    
    #Add tissot circles with 500 km radius
    ax.tissot(facecolor='orange', alpha=0.4, rad_km=500, lons=lons, lats=lats)
    
    #Title including projection name and proj string
    crs_name = str(crs).split(' ')[0].split('+')[-1]
    ax.set_title('%s EPSG:%s\nproj string: "%s"' % (crs_name, crs.to_epsg(), crs.to_proj4()), fontsize=10)

#Define projecitons to plot
#Full list is here: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html
crs_list = [ccrs.NorthPolarStereo(), ccrs.LambertAzimuthalEqualArea(), ccrs.LambertCylindrical(), ccrs.Robinson(), \
            ccrs.EqualEarth(), ccrs.Sinusoidal(), ccrs.Orthographic(), ccrs.Mercator(), \
            ccrs.TransverseMercator(), ccrs.TransverseMercator(central_longitude=-71), ccrs.UTM(19)]

for crs in crs_list:
    tissotplot(crs);
../_images/98aad2a36393da51276c80afeabea393d69a667b878787ca98e5c6626cb4d4ac.png ../_images/da0aada1c314e4dfc416a9f222a0c9716766a9bb58a332f62f737a13e69ec924.png ../_images/72e1156f99c0f89fb3aed3f94d60b87d53726e1c36ea888e59b8ce47f0630346.png ../_images/aa5e51752b044c595d6c1338c0e66151e14c05c04fbf01eb3061563a1793b8b4.png ../_images/8f214b8b49ef0b28d4ba5247501cb0783ba2a751961f72dfcf9a495ae6db4513.png ../_images/6464aaacbe27b4ea5fa89c5d14aa0f9a7bbe6787b758af6a1808745cb17c1b7b.png ../_images/ba0e22e10b8d989d5c7df2d82ae0d50255a040f742d364170edb3daf42f31c1f.png ../_images/4a379a5d3920b1e6014b3eda94811c52105d4c15c28734190450860eb4a1ff8c.png ../_images/d29baa099ae959344e7b514ba52fc862417614b66f48cf8743b5bfdb7dbc65dd.png ../_images/d38e4e93cbe14880df6928078f71681c84f6e81024d6895cb6cab73247974447.png ../_images/f628cbdcdb4a261bf937a1f9984085af0cd93c64153a8ce7b6102577253660f9.png

7.3. Converting coordinates with Pyproj#

In this part, we will donwload a sample scene of Sentinel-2, and work with its coordinates. To do so, we will use pystac_clinent to retrieve the metadata for a Sentinel-2 scene. We will cover pystac_client later in the semester, here you should just know that the following two cells will 1) import the Client, and 2) retrieve the metadata for a specific scene of Sentinel-2.

from pystac_client import Client

client = Client.open("https://earth-search.aws.element84.com/v1")
search_results = client.search(
    ids = ["S2B_37DFA_20230224_0_L2A"]
)
items = search_results.item_collection()
item = items[0]

item

Let’s find the CRS of this scene. The crs is recorded as EPSG numbers in the properties of the item:

item.properties["proj:epsg"]

32737

item

The STAC metadata contains a geometry for each item. This geometry is recorded in EPSG:4326 (check the specification here). Now, we will convert this from 4326 to the crs of the scene.

To transfrom coordinates from one crs to another one, we will creates a transformer using pyproj, and then input each point to transform it.

source_crs = 'epsg:4326' # Global lat-lon coordinate system used by `geometry` in the metadata
target_crs = f'epsg:{item.properties["proj:epsg"]}' # Coordinate system of the S-2 scene

latlon_to_s2_transformer = pyproj.Transformer.from_crs(source_crs, target_crs)

Let’s retrieve the geometry from the item metadata:

geometry = item.geometry

geometry

{'type': 'Polygon',
 'coordinates': [[[41.77852765935361, -71.18173043681364],
   [44.81460295895618, -71.11191009560223],
   [45.122130956889094, -72.09098682450592],
   [41.926326307987814, -72.16491345575483],
   [41.77852765935361, -71.18173043681364]]]}

Now, create an empty geometry for our projected geometry:

projected_geometry = {
    "type": "Polygon",
    "coordinates" : [[]]
}

projected_coordinates = []
for point in geometry["coordinates"][0]:
    projected_x, projected_y = latlon_to_s2_transformer.transform(point[0], point[1])
    projected_coordinates.append([projected_x, projected_y])

projected_geometry["coordinates"] = [projected_coordinates]

projected_geometry

{'type': 'Polygon',
 'coordinates': [[[-5034706.161102917, 22347375.738133498],
   [-4631227.676135911, 22123291.71715717],
   [-4542728.576578871, 22192903.107910063],
   [-4958974.741806454, 22435202.22502178],
   [-5034706.161102917, 22347375.738133498]]]}

7.4. Use Shapely to work with geometries#

Another useful Python package, is Shapely for analysis and manipulation of geometric features. Here, we will use it to investigate our geometries. Note that Shapely doesn’t have a notion of crs, and only look at the object as a geomtric object.

Shapely has a reach set of objects, and here we will use the Polygon one. You can create a Shapely Polygon using the following:

shapely_geo = shapely.geometry.shape(geometry)

shapely_geo
../_images/30439791a464a6ccd7c1eb231902a3fa5530fc5f4d9baf2f45bc452c2a8b4a8f.svg

Let’s retrieve a property of the this Polygon, such as area:

shapely_geo.area

3.0735047508188207

We can do the same using the projected geometry:

shapely_projected_geo = shapely.geometry.shape(projected_geometry)

shapely_projected_geo
../_images/e622634551cc7842050aa2c5064909e844fcabefd7fd5785d78824ba0e53d7df.svg
shapely_projected_geo.area

51412513482.37867

As you can see, the areas are very different. In this case, geometry is in lat/lon coordinates, and projected_geometry is in UTM coordinates (meter). Therefore, their corresponding areas are in their native units.

But we can confirm that these two polygons are indeed the same:

source_crs = f'epsg:{item.properties["proj:epsg"]}' # Coordinate system of the S-2 scene
target_crs = 'epsg:4326' # Global lat-lon coordinate system used by `geometry` in the metadata

s2_to_latlon_transformer = pyproj.Transformer.from_crs(source_crs, target_crs)

reprojected_geometry = {
    "type": "Polygon",
    "coordinates" : [[]]
}

reprojected_coordinates = []
for point in projected_geometry["coordinates"][0]:
    projected_x, projected_y = s2_to_latlon_transformer.transform(point[0], point[1])
    reprojected_coordinates.append([projected_x, projected_y])

reprojected_geometry["coordinates"] = [reprojected_coordinates]

reprojected_geometry

{'type': 'Polygon',
 'coordinates': [[[41.77852765935368, -71.18173043681335],
   [44.814602958956264, -71.1119100956021],
   [45.12213095688916, -72.0909868245058],
   [41.926326307987836, -72.16491345575456],
   [41.77852765935368, -71.18173043681335]]]}

geometry

{'type': 'Polygon',
 'coordinates': [[[41.77852765935361, -71.18173043681364],
   [44.81460295895618, -71.11191009560223],
   [45.122130956889094, -72.09098682450592],
   [41.926326307987814, -72.16491345575483],
   [41.77852765935361, -71.18173043681364]]]}

While there are minor differences in the reprojected geometry vs the original geometry, they look overall the same. These differences are due to transformation errors.