Creating maps with Cartopy
DS Python for GIS and Geoscience
October, 2020© 2020, Joris Van den Bossche and Stijn Van Hoey. Licensed under CC BY-SA 4.0 Creative Commons Adapted from material from Phil Elson and Ryan Abernathey (see below).
Making maps is a fundamental part of GIS and geoscience research. Maps differ from regular figures in the following principle ways:
The maps we have made up to now, for example using the GeoPandas plot()
method, assume the data can be plotted as is on the figure. This works fine if your data is already in projected coordinates, and has a limited extent (small study area). When mapping data of a larger area of the full globe, properly taking into account the projection becomes much more important!
In this notebook, we will learn about Cartopy, one of the most common packages for making maps within python.
Lots of the material in this notebook was adopted from https://earth-env-data-science.github.io/intro.html by Ryan Abernathey, which itself was adopted from Phil Elson's excellent Cartopy Tutorial. Phil is the creator of Cartopy and published his tutorial under an open license, meaning that we can copy, adapt, and redistribute it as long as we give proper attribution. THANKS PHIL AND RYAN! 👏👏👏
Our two most common media are flat:
A map projection (or more commonly refered to as just "projection") is:
a systematic transformation of the latitudes and longitudes of locations from the surface of a sphere or an ellipsoid into locations on a plane. [Wikipedia: Map projection].
There are many different ways to make a projection, and we will not attempt to explain all of the choices and tradeoffs here. Instead, you can read Phil's original tutorial for a great overview of this topic. Instead, we will dive into the more practical sides of Caropy usage.
https://scitools.org.uk/cartopy/docs/latest/
Cartopy makes use of the powerful PROJ.4, numpy and shapely libraries and includes a programatic interface built on top of Matplotlib for the creation of publication quality maps.
Key features of cartopy are its object oriented projection definitions, and its ability to transform points, lines, vectors, polygons and images between those projections.
In Cartopy, each projection is a class. Most classes of projection can be configured in projection-specific ways, although Cartopy takes an opinionated stance on sensible defaults.
Let's create a Plate Carree projection instance.
To do so, we need cartopy's crs module. This is typically imported as ccrs
(Cartopy Coordinate Reference Systems).
import cartopy.crs as ccrs
import cartopy
Cartopy's projection list tells us that the Plate Carree projection is available with the ccrs.PlateCarree
class: https://scitools.org.uk/cartopy/docs/latest/crs/projections.html
Note: we need to instantiate the class in order to do anything projection-y with it!
ccrs.PlateCarree()
<cartopy.crs.PlateCarree object at 0x7fb34d267860>
Cartopy optionally depends upon matplotlib, and each projection knows how to create a matplotlib Axes (or AxesSubplot) that can represent itself.
The Axes that the projection creates is a cartopy.mpl.geoaxes.GeoAxes. This Axes subclass overrides some of matplotlib's existing methods, and adds a number of extremely useful ones for drawing maps:
import matplotlib.pyplot as plt
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
ax
<GeoAxesSubplot:>
That was a little underwhelming, but we can see that the Axes created is indeed one of those GeoAxes[Subplot] instances.
One of the most useful methods that this class adds on top of the standard matplotlib Axes class is the coastlines
method. With no arguments, it will add the Natural Earth 1:110,000,000
scale coastline data to the map.
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb2ec26f7f0>
Projection classes have options we can use to customize the map
ccrs.PlateCarree?
Init signature: ccrs.PlateCarree(central_longitude=0.0, globe=None) Docstring: The abstract class which denotes cylindrical projections where we want to allow x values to wrap around. Init docstring: Parameters ---------- proj4_params: iterable of key-value pairs The proj4 parameters required to define the desired CRS. The parameters should not describe the desired elliptic model, instead create an appropriate Globe instance. The ``proj4_params`` parameters will override any parameters that the Globe defines. globe: :class:`~cartopy.crs.Globe` instance, optional If omitted, the default Globe instance will be created. See :class:`~cartopy.crs.Globe` for details. File: ~/miniconda3/envs/DS-geospatial-test/lib/python3.8/site-packages/cartopy/crs.py Type: ABCMeta Subclasses:
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree(central_longitude=180)})
ax.coastlines()
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb2eee6bdf0>
The cartopy.mpl.geoaxes.GeoAxes class adds a number of useful methods.
Let's take a look at:
set_global - zoom the map out as much as possible
set_extent - zoom the map to the given bounding box
gridlines - add a graticule (and optionally labels) to the axes
coastlines - add Natural Earth coastlines to the axes
stock_img - add a low-resolution Natural Earth background image to the axes
imshow - add an image (numpy array) to the axes
add_geometries - add a collection of geometries (Shapely / GeoPandas) to the axes
projections = [ccrs.PlateCarree(),
ccrs.Robinson(),
ccrs.Mercator(),
ccrs.Orthographic(),
ccrs.InterruptedGoodeHomolosine()
]
for proj in projections:
fig, ax = plt.subplots(subplot_kw={'projection': proj})
ax.stock_img()
ax.coastlines()
ax.set_title(f'{type(proj)}')
To create a regional map, we use the set_extent
method of GeoAxis to limit the size of the region.
ax.set_extent?
Signature: ax.set_extent(extents, crs=None) Docstring: Set the extent (x0, x1, y0, y1) of the map in the given coordinate system. If no crs is given, the extents' coordinate system will be assumed to be the Geodetic version of this axes' projection. Parameters ---------- extents Tuple of floats representing the required extent (x0, x1, y0, y1). File: ~/miniconda3/envs/DS-geospatial-test/lib/python3.8/site-packages/cartopy/mpl/geoaxes.py Type: method
central_lon, central_lat = -10, 45
extent = [-40, 20, 30, 60]
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.Orthographic(central_lon, central_lat)})
ax.set_extent(extent)
ax.gridlines()
ax.coastlines(resolution='50m')
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb2ec07c130>
To give our map more styles and details, we add cartopy.feature
objects.
Many useful features are built in. These "default features" are at coarse (110m) resolution.
Name | Description |
---|---|
cartopy.feature.BORDERS |
Country boundaries |
cartopy.feature.COASTLINE |
Coastline, including major islands |
cartopy.feature.LAKES |
Natural and artificial lakes |
cartopy.feature.LAND |
Land polygons, including major islands |
cartopy.feature.OCEAN |
Ocean polygons |
cartopy.feature.RIVERS |
Single-line drainages, including lake centerlines |
cartopy.feature.STATES |
(limited to the United States at this scale) |
Below we illustrate these features in a customized map of North America.
import cartopy.feature as cfeature
import numpy as np
central_lat = 37.5
central_lon = -96
extent = [-120, -70, 24, 50.5]
central_lon = np.mean(extent[:2])
central_lat = np.mean(extent[2:])
plt.figure(figsize=(12, 6))
ax = plt.axes(projection=ccrs.AlbersEqualArea(central_lon, central_lat))
ax.set_extent(extent)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.LAND, edgecolor='black')
ax.add_feature(cartopy.feature.LAKES, edgecolor='black')
ax.add_feature(cartopy.feature.RIVERS)
ax.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x7fb2ddda6ee0>
If we want higher-resolution features, Cartopy can automatically download and create them from the Natural Earth Data database or the GSHHS dataset database.
You will typically have your own data to add to a map, of course. To add vector data, we can add single Shapely geometries with ax.add_geometries()
, or combine with GeoPandas.
The GeoPandas plot()
method works if:
ax=
keywordtransform=
in the plot()
method.import geopandas
cities = geopandas.read_file("zip://./data/ne_110m_populated_places.zip")
proj = ccrs.AlbersEqualArea()
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': proj})
# Reproject the cities dataset to the corresponding CRS of the cartopy projection
cities_albers = cities.to_crs(proj.proj4_init)
cities_albers.plot(ax=ax)
ax.coastlines()
ax.set_global()
ax.gridlines(draw_labels=True)
<cartopy.mpl.gridliner.Gridliner at 0x7fb2dddb4fd0>
proj = ccrs.AlbersEqualArea()
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': proj})
# cities is in geographic lat/lon data, so we can specify PlateCarree as the cartopy equivalent
cities.plot(ax=ax, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_global()
ax.gridlines(draw_labels=True)
<cartopy.mpl.gridliner.Gridliner at 0x7fb2ddeca9d0>
Or we can also plot specific points:
proj = ccrs.AlbersEqualArea()
fig, ax = plt.subplots(figsize=(10, 10), subplot_kw={'projection': proj})
ax.plot(0, 0, marker='o', color='red', markersize=12, transform=ccrs.PlateCarree())
ax.coastlines()
ax.set_global()
ax.gridlines()
<cartopy.mpl.gridliner.Gridliner at 0x7fb2ddde8fa0>
Because our map is a matplotlib axis, we can use all the familiar maptplotlib commands to make plots.
By default, the map extent will be adjusted to match the data. We can override this with the .set_global
or .set_extent
commands.
Another example to show the difference of adding a transform or not:
# create some test data
new_york = dict(lon=-74.0060, lat=40.7128)
honolulu = dict(lon=-157.8583, lat=21.3069)
lons = [new_york['lon'], honolulu['lon']]
lats = [new_york['lat'], honolulu['lat']]
Key point: the data also have to be transformed to the projection space.
This is done via the transform=
keyword in the plotting method. The argument is another cartopy.crs
object.
If you don't specify a transform, Cartopy assume that the data is using the same projection as the underlying GeoAxis.
From the Cartopy Documentation
The core concept is that the projection of your axes is independent of the coordinate system your data is defined in. The
projection
argument is used when creating plots and determines the projection of the resulting plot (i.e. what the plot looks like). Thetransform
argument to plotting functions tells Cartopy what coordinate system your data are defined in.
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
ax.plot(lons, lats, label='Equirectangular straight line')
ax.plot(lons, lats, label='Great Circle', transform=ccrs.Geodetic())
ax.coastlines()
ax.legend()
ax.set_global()
The same principles apply to 2D data. Below we create some example data defined in regular lat / lon coordinates.
import numpy as np
lon = np.linspace(-80, 80, 25)
lat = np.linspace(30, 70, 25)
lon2d, lat2d = np.meshgrid(lon, lat)
data = np.cos(np.deg2rad(lat2d) * 4) + np.sin(np.deg2rad(lon2d) * 4)
plt.contourf(lon2d, lat2d, data)
<matplotlib.contour.QuadContourSet at 0x7fb2ddd5f7c0>
Now we create a PlateCarree
projection and plot the data on it without any transform
keyword.
This happens to work because PlateCarree
is the simplest projection of lat / lon data.
fig, ax = plt.subplots(subplot_kw={'projection': ccrs.PlateCarree()})
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data)
<cartopy.mpl.contour.GeoContourSet at 0x7fb2ddd1da90>
However, if we try the same thing with a different projection, we get the wrong result.
projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)
fig, ax = plt.subplots(subplot_kw={'projection': projection})
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data)
<cartopy.mpl.contour.GeoContourSet at 0x7fb2ddcc7df0>
To fix this, we need to pass the correct transform argument to contourf
:
projection = ccrs.RotatedPole(pole_longitude=-177.5, pole_latitude=37.5)
fig, ax = plt.subplots(subplot_kw={'projection': projection})
ax.set_global()
ax.coastlines()
ax.contourf(lon, lat, data, transform=ccrs.PlateCarree())
<cartopy.mpl.contour.GeoContourSet at 0x7fb2ddc79970>
We can plot a satellite image easily on a map if we know its extent
! wget https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif
--2020-10-04 22:26:01-- https://github.com/mapbox/rasterio/raw/master/tests/data/RGB.byte.tif Resolving github.com (github.com)... 140.82.121.3 Connecting to github.com (github.com)|140.82.121.3|:443... connected. HTTP request sent, awaiting response... 302 Found Location: https://raw.githubusercontent.com/mapbox/rasterio/master/tests/data/RGB.byte.tif [following] --2020-10-04 22:26:02-- https://raw.githubusercontent.com/mapbox/rasterio/master/tests/data/RGB.byte.tif Resolving raw.githubusercontent.com (raw.githubusercontent.com)... 151.101.36.133 Connecting to raw.githubusercontent.com (raw.githubusercontent.com)|151.101.36.133|:443... connected. HTTP request sent, awaiting response... 200 OK Length: 1743350 (1,7M) [image/tiff] Saving to: ‘RGB.byte.tif.2’ RGB.byte.tif.2 100%[===================>] 1,66M 1,45MB/s in 1,1s 2020-10-04 22:26:03 (1,45 MB/s) - ‘RGB.byte.tif.2’ saved [1743350/1743350]
import rasterio
import rasterio.plot
with rasterio.open('RGB.byte.tif') as src:
img_extent = rasterio.plot.plotting_extent(src)
img = src.read()
print(src.meta)
img = rasterio.plot.reshape_as_image(img)
{'driver': 'GTiff', 'dtype': 'uint8', 'nodata': 0.0, 'width': 791, 'height': 718, 'count': 3, 'crs': CRS.from_epsg(32618), 'transform': Affine(300.0379266750948, 0.0, 101985.0, 0.0, -300.041782729805, 2826915.0)}
proj = ccrs.UTM(zone=18, southern_hemisphere=False)
proj
<cartopy.crs.UTM object at 0x7fb2ddcaaef0>
fig, ax = plt.subplots(figsize=(8, 12), subplot_kw={'projection': proj})
ax.imshow(img, extent=img_extent, transform=proj)
ax.coastlines(resolution='10m', color='red', linewidth=1)
<cartopy.mpl.feature_artist.FeatureArtist at 0x7fb2ec26fb80>
Cartopy transforms can be passed to xarray! This creates a very quick path for creating professional looking maps from netCDF data.
import xarray as xr
url = 'http://www.esrl.noaa.gov/psd/thredds/dodsC/Datasets/noaa.ersst.v5/sst.mnmean.nc'
ds = xr.open_dataset(url, drop_variables=['time_bnds'])
ds
<xarray.Dataset> Dimensions: (lat: 89, lon: 180, time: 2000) Coordinates: * lat (lat) float32 88.0 86.0 84.0 82.0 80.0 ... -82.0 -84.0 -86.0 -88.0 * lon (lon) float32 0.0 2.0 4.0 6.0 8.0 ... 350.0 352.0 354.0 356.0 358.0 * time (time) datetime64[ns] 1854-01-01 1854-02-01 ... 2020-08-01 Data variables: sst (time, lat, lon) float32 ... Attributes: climatology: Climatology is based on 1971-2000 SST, X... description: In situ data: ICOADS2.5 before 2007 and ... keywords_vocabulary: NASA Global Change Master Directory (GCM... keywords: Earth Science > Oceans > Ocean Temperatu... instrument: Conventional thermometers source_comment: SSTs were observed by conventional therm... geospatial_lon_min: -1.0 geospatial_lon_max: 359.0 geospatial_laty_max: 89.0 geospatial_laty_min: -89.0 geospatial_lat_max: 89.0 geospatial_lat_min: -89.0 geospatial_lat_units: degrees_north geospatial_lon_units: degrees_east cdm_data_type: Grid project: NOAA Extended Reconstructed Sea Surface ... original_publisher_url: http://www.ncdc.noaa.gov References: https://www.ncdc.noaa.gov/data-access/ma... source: In situ data: ICOADS R3.0 before 2015, N... title: NOAA ERSSTv5 (in situ only) history: created 07/2017 by PSD data using NCEI's... institution: This version written at NOAA/ESRL PSD: o... citation: Huang et al, 2017: Extended Reconstructe... platform: Ship and Buoy SSTs from ICOADS R3.0 and ... standard_name_vocabulary: CF Standard Name Table (v40, 25 January ... processing_level: NOAA Level 4 Conventions: CF-1.6, ACDD-1.3 metadata_link: :metadata_link = https://doi.org/10.7289... creator_name: Boyin Huang (original) date_created: 2017-06-30T12:18:00Z (original) product_version: Version 5 creator_url_original: https://www.ncei.noaa.gov license: No constraints on data access or use comment: SSTs were observed by conventional therm... summary: ERSST.v5 is developed based on v4 after ... dataset_title: NOAA Extended Reconstructed SST V5 data_modified: 2020-10-03 DODS_EXTRA.Unlimited_Dimension: time
array([ 88., 86., 84., 82., 80., 78., 76., 74., 72., 70., 68., 66., 64., 62., 60., 58., 56., 54., 52., 50., 48., 46., 44., 42., 40., 38., 36., 34., 32., 30., 28., 26., 24., 22., 20., 18., 16., 14., 12., 10., 8., 6., 4., 2., 0., -2., -4., -6., -8., -10., -12., -14., -16., -18., -20., -22., -24., -26., -28., -30., -32., -34., -36., -38., -40., -42., -44., -46., -48., -50., -52., -54., -56., -58., -60., -62., -64., -66., -68., -70., -72., -74., -76., -78., -80., -82., -84., -86., -88.], dtype=float32)
array([ 0., 2., 4., 6., 8., 10., 12., 14., 16., 18., 20., 22., 24., 26., 28., 30., 32., 34., 36., 38., 40., 42., 44., 46., 48., 50., 52., 54., 56., 58., 60., 62., 64., 66., 68., 70., 72., 74., 76., 78., 80., 82., 84., 86., 88., 90., 92., 94., 96., 98., 100., 102., 104., 106., 108., 110., 112., 114., 116., 118., 120., 122., 124., 126., 128., 130., 132., 134., 136., 138., 140., 142., 144., 146., 148., 150., 152., 154., 156., 158., 160., 162., 164., 166., 168., 170., 172., 174., 176., 178., 180., 182., 184., 186., 188., 190., 192., 194., 196., 198., 200., 202., 204., 206., 208., 210., 212., 214., 216., 218., 220., 222., 224., 226., 228., 230., 232., 234., 236., 238., 240., 242., 244., 246., 248., 250., 252., 254., 256., 258., 260., 262., 264., 266., 268., 270., 272., 274., 276., 278., 280., 282., 284., 286., 288., 290., 292., 294., 296., 298., 300., 302., 304., 306., 308., 310., 312., 314., 316., 318., 320., 322., 324., 326., 328., 330., 332., 334., 336., 338., 340., 342., 344., 346., 348., 350., 352., 354., 356., 358.], dtype=float32)
array(['1854-01-01T00:00:00.000000000', '1854-02-01T00:00:00.000000000', '1854-03-01T00:00:00.000000000', ..., '2020-06-01T00:00:00.000000000', '2020-07-01T00:00:00.000000000', '2020-08-01T00:00:00.000000000'], dtype='datetime64[ns]')
[32040000 values with dtype=float32]
sst = ds.sst.sel(time='2000-01-01', method='nearest')
fig = plt.figure(figsize=(9,6))
ax = plt.axes(projection=ccrs.Robinson())
ax.coastlines()
ax.gridlines()
sst.plot(ax=ax, transform=ccrs.PlateCarree(),
vmin=2, vmax=30, cbar_kwargs={'shrink': 0.4})
<matplotlib.collections.QuadMesh at 0x7fb2d738cd90>
Browse the Cartopy Gallery to learn about all the different types of data and plotting methods available!