.ipynb .pdf

repository open issue suggest edit

Contents

Introduction to Cartopy

---

Overview

1. Basic concepts: map projections and GeoAxes

2. Explore some of Cartopy’s map projections

3. Create regional maps

This tutorial will lead you through some basics of creating maps with specified projections with Cartopy, and adding geographic features like coastlines and borders.

Later tutorials will focus on how to plot data on map projections.

Prerequisites

Concepts	Importance	Notes
Matplotlib	Necessary

• Time to learn: 30 minutes

---

Imports

import warnings

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

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

ERROR 1: PROJ: proj_create_from_database: Open of /knight/anaconda_aug22/envs/aug22_env/share/proj failed

---

Motivation: consider the following DataFrame:

loc= pd.read_csv('/spare11/atm350/common/data/locs.csv')

Examine the first few lines of the DataFrame

loc

	lat	lon
0	42.040360	-77.237260
1	42.182270	-74.801390
2	43.019940	-78.135660
3	41.528750	-73.945270
4	42.223220	-75.668520
...	...	...
121	44.393236	-73.858829
122	43.532408	-75.158597
123	43.485073	-73.423071
124	43.228680	-76.842610
125	42.855040	-77.847760

126 rows × 2 columns

Exercise

Create two Pandas Series objects: one for latitude; the second for longitude. Name them lats and lons, respectively.

# Write your code below. Uncomment the following cell only to check your code!

# %load /spare11/atm350/common/mar02/cartopy_1.py
lats = loc['lat']
lons = loc['lon']

Use Matplotlib’s ax.scatter function to plot these points on a figure.

fig = plt.figure(figsize=(11,8.5))
ax = fig.add_subplot(1,1,1)
ax.scatter(lons, lats)
ax.set_title("What does this make you think of??")
ax.set_xlabel('Longitudes')
ax.set_ylabel('Latitudes');

We can get a sense that this looks like a map! However, all we have here are x- and y- coordinate pairs, which just so happen to have values that resemble longitudes and latitudes, in degrees. Matplotlib, by itself, does not provide mapmaking (i.e., cartographic) functionality. We will use a separate package, Cartopy, which extends Matplotlib’s two-dimensional Axes into a cartographic context … i.e., geo-referenced GeoAxes.

Recall that in Matplotlib, what we might tradtionally term a figure consists of two key components: a figure and an associated subplot axes instance.

By virtue of importing Cartopy, we can now convert the axes into a GeoAxes by specifying a map projection that we have imported from Cartopy’s Coordinate Reference System class (ccrs). This will effectively georeference the subplot.

Create a map with a specified projection

Here we’ll create a GeoAxes that uses the PlateCarree projection (basically a global lat-lon map projection, which translates from French to “flat square” in English, where each point is equally spaced in terms of degrees).

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 one of Cartopy’s map projections that 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 of them is coastlines, which is a callable GeoAxes method that can be plotted directly on our subplot.

ax.coastlines()

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

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 we’ve imported as cfeature. These are also shapefiles, downloaded on initial request from https://www.naturalearthdata.com/. Once downloaded, they typically “live” in your ~/.local/share/cartopy directory (note the ~ represents your home directory).

We add them to our subplot via the add_feature method. We can define attributes for them in a manner similar to Matplotlib’s plot method. A list of the various Natural Earth shapefiles appears in https://scitools.org.uk/cartopy/docs/latest/matplotlib/feature_interface.html .

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 0x1532091be680>

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)

This time, we’ll define an object to store our projection definition. Any time we wish to use this particular projection later in the notebook, we can use the object name rather than 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

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

Add a fancy background image to the map.

ax.stock_img()
fig

Lambert Azimuthal Equal Area Projection

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

Now, let’s go back to PlateCarree, but let’s 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. Plot lat/lon lines as well.

Reference for Natural Earth’s three resolutions (10m, 50m, 110m; higher is coarser): 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

Note the in the set_extent call, we specified 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. Let’s try a couple different projections.

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

Set the domain for defining the plot region. We will use this in the set_extent line below. Since these coordinates are expressed in degrees, they correspond to the PlateCarree projection.

Warning

Be patient: with a limited regional extent as specified here, the highest resolution (10m) shapefiles are used; as a result (as with any GeoAxes object that must be transformed from one coordinate system to another, a subject for a subsequent notebook), this will take longer to plot, particularly if you haven’t previously retrieved these features from the Natural Earth shapefile server.

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 pre-defined Features

Some pre-defined Features exist as cartopy.feature constants. The resolution of these pre-defined Features will depend on the areal extent of your map, which you specify via 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'])
# Ocean feature is commented out ... see note below
# ax.add_feature(cfeature.OCEAN)
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 from Natural Earth, it currently takes much longer to render on the plot (try it yourself to see!). Instead, we take the strategy of first setting the facecolor of the entire Axes 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

Let’s create a new map, but this time use lower-resolution shapefiles from Natural Earth, and also eliminate plotting the country borders.

Notice this is a bit more involved. First we create objects for our lower-resolution shapefiles via the NaturalEarthFeature method from Cartopy’s feature class, and then we add them to the map with add_feature.

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 one, we’ll repeat our hi-res NYS map; 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');

Exercise

Create a figure with multiple subplots for different parts of the world.

---

Summary

• Cartopy allows for georeferencing 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 the figure.

---

What’s Next?

In the next notebook, we will plot data from the NYS Mesonet. Additionally, 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

previous

Datetime notebook, ATM350 Spring 2023

next

Week 8: 3/7, 3/9