In this notebook we'll have a look at some exemplary coordinate reference systems for geospatial data.

When working with geospatial data we usually want to simplify the analysis and work with data in a 2D setting with 2D maps. However, any point of the surface of the earth is actually a point on a 3D shape. Hence, in order to get from 3D to 2D we need to do some kind of transformation (mapping, projection).

I think a pretty helpful way to think about this problem is described nicely in an ESRI blog post:

"Imagine peeling an orange and trying to lay the peel flat on a table. You can get close, but only if you start tearing the peel apart. This is where map projections come in. They tell you how to distort the earth—how to tear and stretch that orange peel—so the parts that are most important to your map get the least distorted and are displayed b est on the flat surface of the map."

Image source: https://www.esri.com/arcgis-blog/products/arcgis-pro/mapping/gcs_vs_pcs/

There are infinitely many ways to peel an orange. Similarly, there is not just one way to map the 3D earth surface into 2D. Different transformations come with different advantages and disadvantages. Each transformation defines a coordinate reference system (CRS). They differ with regards to:

units: are points given in latitude / longitude or in meters
distances: some will distort distances, such that individual countries can appear larger than they really are
local / global: some CRSs are optimized for local applications, some for global applications
use cases: e.g. some early CRSs where optimized for navigation. The shortest way between two points could either show as a straight line on a map or as an arc

To get a better feeling for these differences, let's visualize some CRSs. For this, we will plot a world map based on country boundaries together with graticules (representing latitudes and longitudes) for each CRS. The data is taken from https://www.naturalearthdata.com/downloads/110m-physical-vectors/. The examples are somewhat adapted from https://www.earthdatascience.org/courses/use-data-open-source-python/intro-vector-data-python/spatial-data-vector-shapefiles/intro-to-coordinate-reference-systems-python/.

CRS comparison¶

First, let's load the data. We will load them from local paths. In order to replicate this, you will first need to download the files.

In [1]:

import geopandas as gpd
import matplotlib.pyplot as plt
import numpy as np
from shapely import Point
import contextily as cx

In [2]:

world_boundaries = gpd.read_file('../../data/external/ne_110m_land/ne_110m_land.shp')
graticule = gpd.read_file('../../data/external/ne_110m_graticules_all/ne_110m_graticules_15.shp')

Let's define some CRSs in a dictionary. The dictionary keys are name / description of the CRS, the values are projection strings that can be used with geopandas.

In [3]:

crs_list = {'WGS 84 - World Geodetic System 1984, used in GPS (EPSG:4326)': '+proj=longlat +datum=WGS84 +no_defs +type=crs',
            'Adams Square II': '+proj=adams_ws2 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +type=crs',
            'Azimuthal equidistant': '+proj=aeqd +lat_0=0 +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs',
            'Aitoff': '+proj=aitoff +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs',
            'Robinson': '+proj=robin +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs'}

Now let's visualize our data for each CRS:

In [4]:

for this_key in crs_list.keys():
    
    this_proj_str = crs_list[this_key]
    
    # reproject to crs
    world_boundaries_reproj = world_boundaries.to_crs(this_proj_str)
    graticule_reproj = graticule.to_crs(this_proj_str)
    
    # get crs type
    crs_type = world_boundaries_reproj.crs.type_name
    
    # visualize
    fig, ax = plt.subplots(1, 1, figsize=(15, 8))

    graticule_reproj.plot(ax=ax, color='lightgrey')
    world_boundaries_reproj.plot(ax=ax, color='black')

    # title and axes labels
    ax.set(title=f"{this_key} - {crs_type}");

Shape distortions through reprojection¶

Although we can easily transform points and geometries from one coordinate reference system to another, it is important to keep in mind that reprojections will modify shapes. Let's look at an example in order to get a better feeling for this. For this, we select an arbitrary location (Flamingo, Florida) and create a circular and a rectangular shape around it.

In [5]:

poi = Point([-80.9341253, 25.1418102]) # Flamingo, Florida
poi_circle = poi.buffer(10)
poi_rect = poi.buffer(20, cap_style=3, quad_segs=200, resolution=200)

roi_gpd = gpd.GeoDataFrame([poi, poi_circle, poi_rect], crs='epsg:4326', columns=['geometry'])

The following image shows all three geometries together with a base map:

In [6]:

ax = roi_gpd.plot(
    edgecolor="black",
    facecolor="None",
    categorical=True,
    aspect="equal",
    alpha=0.2,
    figsize=(12, 12),
    legend=True,
    legend_kwds={"loc": "upper left", "frameon": False, "ncol": 1},
)

cx.add_basemap(
    ax, crs=roi_gpd.crs.to_string(), source=cx.providers.CartoDB.Voyager
)

Zooming in, we can better verify that the chose point matches our expectation:

In [7]:

ax = roi_gpd.head(2).plot(
    edgecolor="black",
    facecolor="None",
    categorical=True,
    aspect="equal",
    alpha=0.2,
    figsize=(12, 12),
    legend=True,
    legend_kwds={"loc": "upper left", "frameon": False, "ncol": 1},
)

cx.add_basemap(
    ax, crs=roi_gpd.crs.to_string(), source=cx.providers.CartoDB.Voyager
)

Now let's reproject the shapes into a different CRS and visualize it again with a base map:

In [8]:

this_proj_str = '+proj=robin +lon_0=0 +x_0=0 +y_0=0 +R=6371000 +units=m +no_defs +type=crs'
roi_gpd_reproj = roi_gpd.to_crs(this_proj_str)

In [9]:

ax = roi_gpd_reproj.plot(
    edgecolor="black",
    facecolor="None",
    categorical=True,
    aspect="equal",
    alpha=0.2,
    figsize=(12, 12),
    legend=True,
    legend_kwds={"loc": "upper left", "frameon": False, "ncol": 1},
)

cx.add_basemap(
    ax, crs=roi_gpd_reproj.crs.to_string(), source=cx.providers.CartoDB.Voyager
)

As we can see, the circle is not a circle anymore, but it looks a bit like an Easter egg. And the rectangle is distorted as well. However, for the rectangle it is important to keep in mind that we only reprojected 4 points, but not the edges between the points. Hence, the edges still look like straight lines. In reality they would be curved.

The source of the notebook can be found here

You can also find a good overview blog post about the topic on GEOAWESOME.