“Maps are like campfires – everyone gathers around them, because they allow people to understand complex issues at a glance, and find agreement about how to help the land.” - Sonoma Ecology Center

Prerequisites

In Anaconda Powershell Prompt, install cartopy:

# Install cartopy
conda install -c conda-forge cartopy

Once cartopy is installed, load the following modules:

# Import modules
import numpy as np
import xarray as xr
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.ticker as mticker
import cartopy.crs as ccrs
import cartopy.feature as cfeature
%matplotlib inline

In this section, we will use the global monthly near surface temperature (CESM2_200001-201412.nc from the Community Earth System Model (CESM) again for demonstration. To read the data with xarray:

# Load the netCDF4 file
ds = xr.open_dataset("CESM2_200001-201412.nc", engine="netcdf4")
# Check the data
ds
# Use the latest monthly data
surface_T = ds.tas.isel(time=-1)
surface_T

Making maps with `cartopy`

Making maps is a fundamental part of geoscience and environmental research. Maps differ from regular figures in the following principal ways:

Maps require a projection of geographic coordinates on the 3D Earth to the 2D space of your figure.
Maps often include extra decorations besides the data (e.g. continents, country borders, etc.)

Making mapping is a hard and complicated problem, mostly due to the complexities of projection. In this lecture, we will learn about cartopy, one of the most common packages for making maps within python. Another popular and powerful library is basemap. However, basemap is retiring and being replaced by cartopy.

Projections

A map projection (or more commonly referred 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 (source: map projection wikipedia).

The major problems with map projection are:

The surface of a sphere is topologically different from a 2D surface, therefore we have to cut the sphere somewhere.
A sphere’s surface cannot be represented on a plane without distortion.

There are many different ways to make a projection, and they are beyond the scope of this course. You are encouraged to read Phil Elson’s projections tutorial for a great overview. Instead, we will dive into the more practical sides of caropy usage.

Introduction to `cartopy`

Cartopy makes use of the powerful PROJ, numpy, and shapely modules and includes a programmatic 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.

Let’s start with the simple Plate Carrée (French, for flat square) projection (a special case of equirectangular projection) with cartopy’s crs (coordinate reference systems) module. This is typically imported as ccrs (cartopy crs).

# Plate Carree projection
ccrs.PlateCarree()

# Show more details of the Plate Carree projection
print(ccrs.PlateCarree())

# Check the projection
ccrs.PlateCarree?

For other available cartopy projections, check cartopy coordinate reference systems.

Drawing a simple map

Cartopy depends upon matplotlib, and each projection knows how to create a matplotlib Axes that can represent itself. This Axes overrides some of matplotlib’s existing methods, and adds a number of extremely useful ones for drawing maps.

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with an basic PlateCarree projection style
proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# Add LAND feature to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey', zorder=0)

As you can see, ccrs.PlateCarree() creates an axes with a certain projection. We then use add_feature() to add land (LAND) feature to the axes. Here the pre-defined features are all small-scale (1:110 million) Natural Earth datasets:

BORDERS Country boundaries
COASTLINE Coastline, including major islands
LAKES Natural and artificial lakes
LAND Land polygons, including major islands
OCEAN Ocean polygons
RIVERS Single-line drainages, including lake centerlines
STATES Limited to the United States

Now try:

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with an basic PlateCarree projection style
proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.OCEAN, zorder=0)
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey', zorder=1)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)

You will soon figure out zorder means the order of each layer. Try to change those to different values, see what happens.

In fact, Natural Earth dataset can easily be used by creating an instance of NaturalEarthFeature(). For example:

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with an basic PlateCarree projection style
proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.NaturalEarthFeature(category='physical',
                                           name='land',
                                           scale='110m',
                                           facecolor='gray',
                                           edgecolor='black',
                                           linewidth=1))

And add country border lines (admin_0_countries) to the previous map:

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with an basic PlateCarree projection style
proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.NaturalEarthFeature(category='physical',
                                           name='land',
                                           scale='110m',
                                           facecolor='gray',
                                           edgecolor='black',
                                           linewidth=1))

# Add border lines over countries 
ax.add_feature(cfeature.NaturalEarthFeature(category='cultural',
                                           name='admin_0_countries',
                                           scale='110m',
                                           facecolor='none',
                                           edgecolor='black',
                                           linewidth=0.5))

For more about NaturalEarthFeature(), use help(cfeature.NaturalEarthFeature) or check the cartopy feature interface. For a list of available features in NaturalEarthFeature() visit Natural Earth Features.

Now let’s put all features together:

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with an basic PlateCarree projection style
proj = ccrs.PlateCarree()
ax = plt.axes(projection=proj)

# Add natural features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.OCEAN, zorder=0)
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey', zorder=1)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)

# Add border lines over countries 
ax.add_feature(cfeature.NaturalEarthFeature(category='cultural',
                                           name='admin_0_countries',
                                           scale='110m',
                                           facecolor='none',
                                           edgecolor='black',
                                           linewidth=0.5))

# Add lat/lon gridlines, draw gridlines
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5)

# Manipulate latitude and longitude gridline numbers and spacing
gl.ylocator = mticker.FixedLocator(np.arange(-90,91,30))
gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 30))

Here we use gridlines() to add a graticule (and optionally labels) to the axes, and use xlocator() and ylocator() to manipulate latitude and longitude gridline numbers and spacing.

Using different global projections

Simply change proj from ccrs.PlateCarree() to ccrs.Orthographic(), you will have a map with an Orthographic projection.

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with Orthographic projection style
proj = ccrs.Orthographic()         
ax = plt.axes(projection=proj)

# Add natural features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.OCEAN, zorder=0)
ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey', zorder=1)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)

# Add border lines over countries 
ax.add_feature(cfeature.NaturalEarthFeature(category='cultural',
                                           name='admin_0_countries',
                                           scale='110m',
                                           facecolor='none',
                                           edgecolor='black',
                                           linewidth=0.5))

# Add lat/lon gridlines, draw gridlines
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5)

# Manipulate latitude and longitude gridline numbers and spacing
gl.ylocator = mticker.FixedLocator(np.arange(-90,91,30))
gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 30))

Or play with different projections from a list:

# Define plotting function
def plot_map(my_projection):
    # Create and define the size of a figure object 
    plt.figure(figsize=(5,5), dpi=100)
    
    # Create an axes with Orthographic projection style
    ax = plt.axes(projection=my_projection)

    # Add natural features to axes using cartopy.feature (cfeature)
    ax.add_feature(cfeature.OCEAN, zorder=0)
    ax.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey', zorder=1)
    ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)

    # Add border lines over countries 
    ax.add_feature(cfeature.NaturalEarthFeature(category='cultural',
                                           name='admin_0_countries',
                                           scale='110m',
                                           facecolor='none',
                                           edgecolor='black',
                                           linewidth=0.5))

    # Add lat/lon gridlines, draw gridlines
    gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5)

    # Manipulate latitude and longitude gridline numbers and spacing
    gl.ylocator = mticker.FixedLocator(np.arange(-90,91,30))
    gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 30))
    
    # Add title
    ax.set_title(f'{type(my_projection)}')
    
# Set a list of projections
projections = [ccrs.PlateCarree(),
               ccrs.Robinson(),
               ccrs.Mercator(),
               ccrs.Orthographic(),
               ccrs.InterruptedGoodeHomolosine(),
               ccrs.Gnomonic()
              ]

# Loop the projections and call the plotting function
for proj in projections:
    plot_map(proj)

Check cartopy projection list for a list of available projections supported by cartopy.

Making regional maps

To create a regional map, we use the set_extent method to limit the size of the region. Use ax.set_extent? or help(ax.set_extent) to learn how to use it.

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Set Orthographic projection style
central_lon, central_lat = 114.06, 22.54 # Shenzhen
proj = ccrs.Orthographic(central_lon, central_lat) 

# Create an axes with Orthographic projection style
ax = plt.axes(projection=proj)

# Set a region and plot
extent = [central_lon-10, central_lon+10, central_lat-10, central_lat+10]
ax.set_extent(extent)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)
ax.add_feature(cfeature.RIVERS, edgecolor='blue', zorder=3)

# Add features to axes using methods
ax.coastlines(resolution='10m', linewidth=0.5)
ax.gridlines()

Here we create a regional map with Orthographic projection over Southern China. And we use the coastlines() (same as add_feature() as we previously used) method to add coastlines.

For higher-resolution features, cartopy can automatically download and create them from the Natural Earth Data or the Global Self-consistent, Hierarchical, High-resolution Geography Database database.

# Now use higher-resolution features
rivers_10m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m')

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Set Orthographic projection style
central_lon, central_lat = 114.06, 22.54 # Shenzhen
proj = ccrs.Orthographic(central_lon, central_lat) 

# Create an axes with Orthographic projection style
ax = plt.axes(projection=proj)

# Set a region and plot
extent = [central_lon-10, central_lon+10, central_lat-10, central_lat+10]
ax.set_extent(extent)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)
ax.add_feature(rivers_10m, facecolor='None', edgecolor='blue', linewidth=0.5)

# Add features to axes using coastlines method
ax.coastlines(resolution='10m')
ax.gridlines()

Adding data to the map

Now that we know how to create a map, let’s add our data to it.

Because our map is a matplotlib axes, we can use all the familiar maptplotlib commands to make plots. By default, the map extent will be adjusted to match the data.

# Create some test data
Shenzhen = dict(lon=114.06, lat=22.54)
Guangzhou = dict(lon=113.25, lat=23.13)
lons = [Shenzhen['lon'], Guangzhou['lon']]
lats = [Shenzhen['lat'], Guangzhou['lat']]

# Now use higher-resolution features
rivers_10m = cfeature.NaturalEarthFeature('physical', 'rivers_lake_centerlines', '10m')

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Set Orthographic projection style
central_lon, central_lat = 114.06, 22.54 # Shenzhen
proj = ccrs.Orthographic(central_lon, central_lat) 

# Create an axes with Orthographic projection style
ax = plt.axes(projection=proj)

# Set a region and plot
extent = [central_lon-1, central_lon+1, central_lat-1, central_lat+1]
ax.set_extent(extent)

# Add features to axes using cartopy.feature (cfeature)
ax.add_feature(cfeature.LAKES, edgecolor='blue', facecolor='blue', zorder=2)
ax.add_feature(rivers_10m, facecolor='None', edgecolor='blue', linewidth=0.5)

# Add features to axes using coastlines method
ax.coastlines(resolution='10m')
ax.gridlines()

# Add two points
ax.plot(lons, lats, 'go-', transform=ccrs.PlateCarree())

Notice, 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 assumes that the data is using the same projection as the underlying axes.

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). The transform argument to plotting functions tells cartopy what coordinate system your data are defined in. In our case, lat and lon are defined in degrees and in PlateCarree projection. For more, check understanding the transform and projection keywords.

`Xarray` integration

Cartopy transforms can be passed to xarray, this creates a very quick path for creating professional-looking maps from netCDF data.

# Create and define the size of a figure object 
plt.figure(figsize=(5,5), dpi=100)

# Create an axes with Orthographic projection style
proj = ccrs.Orthographic(central_lon, central_lat) 
ax = plt.axes(projection=proj)

# Plot the surface temperature
surface_T.plot(ax=ax, transform=ccrs.PlateCarree(),
         vmin=250, vmax=300, cbar_kwargs={'shrink': 0.4})

# Add border lines over countries 
ax.add_feature(cfeature.NaturalEarthFeature(category='cultural',
                                           name='admin_0_countries',
                                           scale='110m',
                                           facecolor='none',
                                           edgecolor='black',
                                           linewidth=0.5))

# Add lat/lon gridlines, draw gridlines
gl = ax.gridlines(crs=ccrs.PlateCarree(), linewidth=1, color='black', alpha=0.5)

# Manipulate latitude and longitude gridline numbers and spacing
gl.ylocator = mticker.FixedLocator(np.arange(-90,91,30))
gl.xlocator = mticker.FixedLocator(np.arange(-180, 181, 30))

Masking a map projection

One can also “mask” a cartopy map projection with data from the dataset. This allows us to, essentially, block out specific areas of data on the projection.

# Create and define the size of a figure object 
plt.figure(figsize=(10,10), dpi=100)

# Create an axes with Orthographic projection style
proj = ccrs.Orthographic(central_lon, central_lat) 

# Plot first panel
ax1 = plt.subplot(2, 2, 1, projection=proj)
ax1.coastlines(linewidths=0.5)
# Mask land, lakes, and oceans
ax1.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey',linewidths=0.5)
ax1.add_feature(cfeature.OCEAN)
plt.title("Many features mask")

# Plot second plot
ax2 = plt.subplot(2, 2, 2, projection=proj)
ax2.coastlines(linewidths=0.5)
# Contourf-plot data
surface_T.plot.contourf(ax=ax2, transform=ccrs.PlateCarree(),
               vmin=250, vmax=300, levels=11, cmap='magma',
               add_colorbar=False)
# Mask ocean data by adding ocean feature and changing its zorder
ax2.add_feature(cfeature.OCEAN, zorder=1)
plt.title("Applying ocean mask")

# Plot third plot
ax3 = plt.subplot(2, 2, 3, projection=proj)
ax3.coastlines(linewidths=0.5)
# Contourf-plot data
surface_T.plot.contourf(ax=ax3, transform=ccrs.PlateCarree(),
               vmin=250, vmax=300, levels=11, cmap='magma',
               add_colorbar=False)
# Mask land data by adding land feature and changing its zorder
ax3.add_feature(cfeature.LAND, edgecolor='black', facecolor='grey',
                linewidths=0.5, zorder=1)
plt.title("Applying land mask")

# Plot fourth plot
ax4 = plt.subplot(2, 2, 4, projection=proj)
ax4.coastlines(linewidths=0.5)
# Contourf-plot data
surface_T.plot.contourf(ax=ax4, transform=ccrs.PlateCarree(),
               vmin=250, vmax=300, levels=11, cmap='magma',
               add_colorbar=False)
plt.title("No mask data only")

The notes are modified from the excellent Maps in Scientific Python and Python Tutorial Seminar Series: Cartopy, both available freely online.

In-class exercises

Exercise #1

Go through the notes, make sure you understand the scripts.

Exercise #2

Applying the LambertConformal projection, plot the mean surface temperature in June from 2000 to 2010 near your hometown. Limit your figure to a 5 degrees by 5 degrees domain.

Section 09 Making maps with `cartopy`

Prerequisites

Making maps with `cartopy`

Projections

Introduction to `cartopy`

Drawing a simple map

Using different global projections

Making regional maps

Adding data to the map

`Xarray` integration

Masking a map projection

In-class exercises

Exercise #1

Exercise #2

Further reading

Section 09 Making maps with cartopy

Prerequisites

Making maps with cartopy

Projections

Introduction to cartopy

Drawing a simple map

Using different global projections

Making regional maps

Adding data to the map

Xarray integration

Masking a map projection

In-class exercises

Exercise #1

Exercise #2

Further reading

Section 09 Making maps with `cartopy`

Making maps with `cartopy`

Introduction to `cartopy`

`Xarray` integration