Introduction to Cartopy

---

Overview

The concepts covered in this section include:

1. Learning core Cartopy concepts: map projections and GeoAxes

2. Exploring some of Cartopy’s map projections

3. Creating regional maps

This tutorial will lead you through some basics of creating maps with specified projections using Cartopy, and adding geographical features (like coastlines and borders) to those maps.

Plotting data on map projections will be covered in later tutorials.

Prerequisites

Concepts	Importance	Notes
Matplotlib	Necessary

• Time to learn: 30 minutes

---

Imports

Here, we import the main libraries of Cartopy: crs and feature. In addition, we import numpy, as well as matplotlib’s pyplot interface. Finally, we import a library called warnings, and use it to remove extraneous warnings that Cartopy produces in later examples.

import warnings

import matplotlib.pyplot as plt
import numpy as np
from cartopy import crs as ccrs, feature as cfeature

#  Suppress warnings issued by Cartopy when downloading data files
warnings.filterwarnings('ignore')

---

Basic concepts: map projections and GeoAxes

Extend Matplotlib’s axes into georeferenced GeoAxes

Recall from earlier tutorials that a figure in Matplotlib has two elements: a Figure object, and a list of one or more Axes objects (subplots).

Since we imported cartopy.crs, we now have access to Cartopy’s Coordinate Reference System, which contains many geographical projections. We can specify one of these projections for an Axes object to convert it into a GeoAxes object. This will effectively georeference the subplot. Examples of converting Axes objects into GeoAxes objects can be found later in this section.

Create a map with a specified projection

In this example, we’ll create a GeoAxes object that uses the PlateCarree projection. PlateCarree is a global lat-lon map projection in which each point is evenly spaced in terms of degrees. The name “Plate Carree” is French for “flat square”.

fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=ccrs.PlateCarree(central_longitude=-75))
ax.set_title("A Geo-referenced subplot, Plate Carree projection");

Although the figure seems empty, it has, in fact, been georeferenced using a map projection; this projection is provided by Cartopy’s crs (coordinate reference system) class. We can now add in cartographic features, in the form of shapefiles, to our subplot. One such cartographic feature is coastlines, which can be added to our subplot using the callable GeoAxes method simply called coastlines.

ax.coastlines()

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f237b39dc40>

Info

To get the figure to display again with the features that we’ve added since the original display, just type the name of the Figure object in its own cell.

fig

Add cartographic features to the map

Cartopy provides other cartographic features via its features class, which was imported at the beginning of this page, under the name cfeature. These cartographic features are laid out as data in shapefiles. The shapefiles are downloaded when their cartographic features are used for the first time in a script or notebook, and they are downloaded from https://www.naturalearthdata.com/. Once downloaded, they “live” in your ~/.local/share/cartopy directory (note the ~ represents your home directory).

We can add these features to our subplot via the add_feature method; this method allows the definition of attributes using arguments, similar to Matplotlib’s plot method. A list of the various Natural Earth shapefiles can be found at https://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html. In this example, we add borders and U. S. state lines to our subplot:

ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='black')
ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')

<cartopy.mpl.feature_artist.FeatureArtist at 0x7f232ff54260>

Once again, referencing the Figure object will re-render the figure in the notebook, now including the two features.

fig

Explore some of Cartopy’s map projections

Info

You can find a list of supported projections in Cartopy, with examples, at https://scitools.org.uk/cartopy/docs/latest/reference/crs.html

Mollweide Projection (often used with global satellite mosaics)

To save typing later, we can define a projection object to store the definition of the map projection. We can then use this object in the projection kwarg of the subplot method when creating a GeoAxes object. This allows us to use this exact projection in later scripts or Jupyter Notebook cells using simply the object name, instead of repeating the same call to ccrs.

fig = plt.figure(figsize=(11, 8.5))
projMoll = ccrs.Mollweide(central_longitude=0)
ax = plt.subplot(1, 1, 1, projection=projMoll)
ax.set_title("Mollweide Projection")

Text(0.5, 1.0, 'Mollweide Projection')

Add in the cartographic shapefiles

This example shows how to add cartographic features to the Mollweide projection defined earlier:

ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue')
fig

Add a fancy background image to the map.

We can also use the stock_img method to add a pre-created background to a Mollweide-projection plot:

ax.stock_img()
fig

Lambert Azimuthal Equal Area Projection

This example is similar to the above example set, except it uses a Lambert azimuthal equal-area projection instead:

fig = plt.figure(figsize=(11, 8.5))
projLae = ccrs.LambertAzimuthalEqualArea(central_longitude=0.0, central_latitude=0.0)
ax = plt.subplot(1, 1, 1, projection=projLae)
ax.set_title("Lambert Azimuthal Equal Area Projection")
ax.coastlines()
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');

Create regional maps

Cartopy’s set_extent method

For this example, let’s create another PlateCarree projection, but this time, we’ll use Cartopy’s set_extent method to restrict the map coverage to a North American view. Let’s also choose a lower resolution for coastlines, just to illustrate how one can specify that. In addition, let’s also plot the latitude and longitude lines.

Natural Earth defines three resolutions for cartographic features, specified as the strings “10m”, “50m”, and “110m”. Only one resolution can be used at a time, and the higher the number, the less detailed the feature becomes. You can view the documentation for this functionality at the following reference link: https://www.naturalearthdata.com/downloads/

projPC = ccrs.PlateCarree()
lonW = -140
lonE = -40
latS = 15
latN = 65
cLat = (latN + latS) / 2
cLon = (lonW + lonE) / 2
res = '110m'

fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=projPC)
ax.set_title('Plate Carree')
gl = ax.gridlines(
    draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'
)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.coastlines(resolution=res, color='black')
ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');

Info

Please note, even though the calls to the subplot method use different projections, the calls to set_extent use PlateCarree. This ensures that the values we passed into set_extent will be transformed from degrees into the values appropriate for the projection we use for the map.

The PlateCarree projection exaggerates the spatial extent of regions closer to the poles. In the following examples, we use set_extent with stereographic and Lambert-conformal projections, which display polar regions more accurately.

projStr = ccrs.Stereographic(central_longitude=cLon, central_latitude=cLat)
fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=projStr)
ax.set_title('Stereographic')
gl = ax.gridlines(
    draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'
)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.coastlines(resolution=res, color='black')
ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');

projLcc = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)
fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=projLcc)
ax.set_title('Lambert Conformal')
gl = ax.gridlines(
    draw_labels=True, linewidth=2, color='gray', alpha=0.5, linestyle='--'
)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.coastlines(resolution='110m', color='black')
ax.add_feature(cfeature.STATES, linewidth=0.3, edgecolor='brown')
# End last line with a semicolon to suppress text output to the screen
ax.add_feature(cfeature.BORDERS, linewidth=0.5, edgecolor='blue');

Info

Lat/lon labeling for projections other than Mercator and PlateCarree is a recent addition to Cartopy. As you can see, work still needs to be done to improve the placement of labels.

Create a regional map centered over New York State

Here we set the domain, which defines the geographical region to be plotted. (This is used in the next section in a set_extent call.) Since these coordinates are expressed in degrees, they correspond to a PlateCarree projection, even though the map projection is set to LambertConformal.

Warning

Be patient; when plotting a small geographical area, the high-resolution “10m” shapefiles are used by default. As a result, these plots take longer to create, especially if the shapefiles are not yet downloaded from Natural Earth. Similar issues can occur whenever a GeoAxes object is transformed from one coordinate system to another. (This will be covered in more detail in a subsequent page.)

latN = 45.2
latS = 40.2
lonW = -80.0
lonE = -71.5
cLat = (latN + latS) / 2
cLon = (lonW + lonE) / 2
projLccNY = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)

Add some predefined features

Some cartographical features are predefined as constants in the cartopy.feature package. The resolution of these features depends on the amount of geographical area in your map, specified by set_extent.

fig = plt.figure(figsize=(15, 10))
ax = plt.subplot(1, 1, 1, projection=projLccNY)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='--')
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.STATES)
ax.add_feature(cfeature.RIVERS)
ax.set_title('New York and Vicinity');

Note:

For high-resolution Natural Earth shapefiles such as this, while we could add Cartopy’s OCEAN feature, it currently takes much longer to render on the plot. You can create your own version of this example, with the OCEAN feature added, to see for yourself how much more rendering time is added. Instead, we take the strategy of first setting the facecolor of the entire subplot to match that of water bodies in Cartopy. When we then layer on the LAND feature, pixels that are not part of the LAND shapefile remain in the water facecolor, which is the same color as the OCEAN.

Use lower-resolution shapefiles from Natural Earth

In this example, we create a new map. This map uses lower-resolution shapefiles from Natural Earth, and also eliminates the plotting of country borders.

This example requires much more code than previous examples on this page. First, we must create new objects associated with lower-resolution shapefiles. This is performed by the NaturalEarthFeature method, which is part of the Cartopy feature class. Second, we call add_feature to add the new objects to our new map.

fig = plt.figure(figsize=(15, 10))
ax = plt.subplot(1, 1, 1, projection=projLccNY)
ax.set_extent((lonW, lonE, latS, latN), crs=projPC)

# The features with names such as cfeature.LAND, cfeature.OCEAN, are higher-resolution (10m)
# shapefiles from the Naturalearth repository.  Lower resolution shapefiles (50m, 110m) can be
# used by using the cfeature.NaturalEarthFeature method as illustrated below.

resolution = '110m'

land_mask = cfeature.NaturalEarthFeature(
    'physical',
    'land',
    scale=resolution,
    edgecolor='face',
    facecolor=cfeature.COLORS['land'],
)
sea_mask = cfeature.NaturalEarthFeature(
    'physical',
    'ocean',
    scale=resolution,
    edgecolor='face',
    facecolor=cfeature.COLORS['water'],
)
lake_mask = cfeature.NaturalEarthFeature(
    'physical',
    'lakes',
    scale=resolution,
    edgecolor='face',
    facecolor=cfeature.COLORS['water'],
)
state_borders = cfeature.NaturalEarthFeature(
    category='cultural',
    name='admin_1_states_provinces_lakes',
    scale=resolution,
    facecolor='none',
)

ax.add_feature(land_mask)
ax.add_feature(sea_mask)
ax.add_feature(lake_mask)
ax.add_feature(state_borders, linestyle='solid', edgecolor='black')
ax.set_title('New York and Vicinity; lower resolution');

A figure with two different regional maps

Finally, let’s create a figure with two subplots. On the first subplot, we’ll repeat the high-resolution New York State map created earlier; on the second, we’ll plot over a different part of the world.

# Create the figure object
fig = plt.figure(
    figsize=(30, 24)
)  # Notice we need a bigger "canvas" so these two maps will be of a decent size

# First subplot
ax = plt.subplot(2, 1, 1, projection=projLccNY)
ax.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax.set_facecolor(cfeature.COLORS['water'])
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linestyle='--')
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.STATES)
ax.set_title('New York and Vicinity')

# Set the domain for defining the second plot region.
latN = 70
latS = 30.2
lonW = -10
lonE = 50
cLat = (latN + latS) / 2
cLon = (lonW + lonE) / 2

projLccEur = ccrs.LambertConformal(central_longitude=cLon, central_latitude=cLat)

# Second subplot
ax2 = plt.subplot(2, 1, 2, projection=projLccEur)
ax2.set_extent([lonW, lonE, latS, latN], crs=projPC)
ax2.set_facecolor(cfeature.COLORS['water'])
ax2.add_feature(cfeature.LAND)
ax2.add_feature(cfeature.COASTLINE)
ax2.add_feature(cfeature.BORDERS, linestyle='--')
ax2.add_feature(cfeature.LAKES, alpha=0.5)
ax2.add_feature(cfeature.STATES)
ax2.set_title('Europe');

An example of plotting data

First, we’ll create a lat-lon grid and define some data on it.

lon, lat = np.mgrid[-180:181, -90:91]
data = 2 * np.sin(3 * np.deg2rad(lon)) + 3 * np.cos(4 * np.deg2rad(lat))
plt.contourf(lon, lat, data)
plt.colorbar();

Plotting data on a Cartesian grid is equivalent to plotting data in the PlateCarree projection, where meridians and parallels are all straight lines with constant spacing. As a result of this simplicity, the global datasets we use often begin in the PlateCarree projection.

Once we create our map again, we can plot these data values as a contour map. We must also specify the transform keyword argument. This is an argument to a contour-plotting method that specifies the projection type currently used by our data. The projection type specified by this argument will be transformed into the projection type specified in the subplot method. Let’s plot our data in the Mollweide projection to see how shapes change under a transformation.

fig = plt.figure(figsize=(11, 8.5))
ax = plt.subplot(1, 1, 1, projection=projMoll)
ax.coastlines()
dataplot = ax.contourf(lon, lat, data, transform=ccrs.PlateCarree())
plt.colorbar(dataplot, orientation='horizontal');

---

Summary

• Cartopy allows for the georeferencing of Matplotlib Axes objects.

• Cartopy’s crs class supports a variety of map projections.

• Cartopy’s feature class allows for a variety of cartographic features to be overlaid on a georeferenced plot or subplot.

---

What’s Next?

In the next notebook, we will delve further into how one can transform data that is defined in one coordinate reference system (crs) so it displays properly on a map that uses a different crs.

Resources and References

1. Cartopy Documentation

2. Full list of projections in Cartopy

3. Maps with Cartopy (Ryan Abernathey)

4. Map Projections (GeoCAT)

5. NCAR xdev Cartopy Tutorial Video