CartoPy

Unidata Logo

Plotting on a Map with CartoPy

Unidata Python Workshop


CartoPy

Questions

  1. How do we plot on a map in Python?
  2. How do I specify a map projection?
  3. How do I tell CartoPy how to reference my data?
  4. How do I add map features to a CartoPy plot?

Objectives

  1. Create a basic figure using CartoPy
  2. Add maps to the figure
  3. Plot georeferenced data on the figure

1. Basic CartoPy Plotting

  • High level API for dealing with maps
  • CartoPy allows you to plot data on a 2D map.
  • Support many different map projections
  • Support for shapefiles from the GIS world
In [1]:
# Set things up
%matplotlib inline

# Importing CartoPy
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt

The simplest plot we can make sets a projection with no parameters. The one below uses the Robinson projection:

In [2]:
# Works with matplotlib's built-in transform support.
fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.Robinson())

# Sets the extent to cover the whole globe
ax.set_global()

# Adds standard background map
ax.stock_img()
Out[2]:
<matplotlib.image.AxesImage at 0x7fa46237c350>

We also have fine-tuned control over the globe used in the projection as well as lots of standard parameters, which depend on individual projections:

In [3]:
# Set up a globe with a specific radius
globe = ccrs.Globe(semimajor_axis=6371000.)

# Set up a Lambert Conformal projection
proj = ccrs.LambertConformal(standard_parallels=[25.0], globe=globe)

fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=proj)

# Sets the extent using a lon/lat box
ax.set_extent([-130, -60, 20, 55])

ax.stock_img()
Out[3]:
<matplotlib.image.AxesImage at 0x7fa46228c150>

Top


2. Adding maps to CartoPy

CartoPy provides a couple helper methods for adding maps to the plot:

In [4]:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

ax.stock_img()

ax.add_feature(cfeature.COASTLINE)

ax.set_extent([-130, -60, 20, 55])

Cartopy also has a lot of built-in support for a variety of map features:

In [5]:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

# Add variety of features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)

# Can also supply matplotlib kwargs
ax.add_feature(cfeature.BORDERS, linestyle=':')
ax.add_feature(cfeature.STATES, linestyle=':')
ax.add_feature(cfeature.LAKES, alpha=0.5)
ax.add_feature(cfeature.RIVERS, edgecolor='tab:green')

ax.set_extent([-130, -60, 20, 55])
/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/50m/physical/ne_50m_ocean.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)
/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/50m/physical/ne_50m_lakes.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)
/home/travis/miniconda/envs/unidata/lib/python3.7/site-packages/cartopy/io/__init__.py:260: DownloadWarning: Downloading: https://naciscdn.org/naturalearth/50m/physical/ne_50m_rivers_lake_centerlines.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)

The map features are available at several different scales depending on how large the area you are covering is. The scales can be accessed using the with_scale method. Natural Earth features are available at 110m, 50m and 10m.

In [6]:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

# Add variety of features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.OCEAN)
ax.add_feature(cfeature.COASTLINE)

# Can also supply matplotlib kwargs
ax.add_feature(cfeature.BORDERS.with_scale('50m'), linestyle=':')
ax.add_feature(cfeature.STATES.with_scale('50m'), linestyle=':')
ax.add_feature(cfeature.LAKES.with_scale('50m'), alpha=0.5)
ax.add_feature(cfeature.RIVERS.with_scale('50m'), edgecolor='tab:green')

ax.set_extent([-130, -60, 20, 55])

You can also grab other features from the Natural Earth project: http://www.naturalearthdata.com/

US Counties

MetPy has US Counties built in at the 20m, 5m, and 500k resolutions.

In [7]:
from metpy.plots import USCOUNTIES
In [8]:
proj = ccrs.LambertConformal(central_longitude=-85.0, central_latitude=45.0)

fig = plt.figure(figsize=(12, 9))
ax1 = fig.add_subplot(1, 3, 1, projection=proj)
ax2 = fig.add_subplot(1, 3, 2, projection=proj)
ax3 = fig.add_subplot(1, 3, 3, projection=proj)

for scale, axis in zip(['20m', '5m', '500k'], [ax1, ax2, ax3]):
    axis.set_extent([270.25, 270.9, 38.15, 38.75], ccrs.Geodetic())
    axis.add_feature(USCOUNTIES.with_scale(scale), edgecolor='black')
Downloading file 'us_counties_20m.dbf' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_20m.dbf' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_20m.shx' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_20m.shx' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_20m.shp' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_20m.shp' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_5m.dbf' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_5m.dbf' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_5m.shx' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_5m.shx' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_5m.shp' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_5m.shp' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_500k.dbf' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_500k.dbf' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_500k.shx' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_500k.shx' to '/home/travis/.cache/metpy/v0.12.2'.
Downloading file 'us_counties_500k.shp' from 'https://github.com/Unidata/MetPy/raw/v0.12.2/staticdata/us_counties_500k.shp' to '/home/travis/.cache/metpy/v0.12.2'.

Top


3. Plotting Data

CartoPy supports all of the matplotlib plotting options you would expect on a map. It handles transforming your data between different coordinate systems transparently, provided you provide the correct information. (More on this later...). To start, let's put a marker at -105, 40:

In [9]:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=2)
ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black')

ax.plot(-105, 40, marker='o', color='tab:red')

ax.set_extent([-130, -60, 20, 55])

So that did not succeed at putting a marker at -105 longitude, 40 latitude (Boulder, CO). Instead, what actually happened is that it put the marker at (-105, 40) in the map projection coordinate system; in this case that's a Lambert Conformal projection, and x,y are assumed in meters relative to the origin of that coordinate system. To get CartoPy to treat it as longitude/latitude, we need to tell it that's what we're doing. We do this through the use of the transform argument to all of the plotting functions.

In [10]:
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=ccrs.LambertConformal())

ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=2)
ax.add_feature(cfeature.STATES, linestyle='--', edgecolor='black')

data_projection = ccrs.PlateCarree()
ax.plot(-105, 40, marker='o', color='tab:red', transform=data_projection)

ax.set_extent([-130, -60, 20, 55])