Python | Tutorial: Intro to Cartopy

Introduction:

There are currently two main Python libraries for plotting geographic data on map: Cartopy and Basemap.  New users should use the Cartopy since support Cartopy will replace Basemap and support for Basemap is expected to wrap up in 2020.  Therefore, this tutorial will focus on Cartopy.  If you want to use Basemap (e.g., for plotting in "3D"), here are a couple of links to existing tutorials and examples to help get you started.  Furthermore, many concepts described here may be useful for understanding Basemap.

• Basemap Examples
• Basemap Tutorials 

Cartopy:

Cartopy is a geospatial plotting library built on top of Numpy and Matplotlib that makes plotting gridded data, shapefiles, and other geographic data on over 30 different map projections.  Furthermore, the Cartopy "transform" functionality makes it straightforward to convert data from one projection to another.  In this tutorial you will learn how to use Cartopy to:

• Plot GFS surface temperature data on a map
• Mask out land surfaces to plot a map of sea surface temperature and ice-cover from the GFS
• Plot a regional map of surface temperature with US states
• Use the transform function to plot GFS data on a different map projection

This tutorial accesses the NOMADS data server using the netCDF4 library, if you are unfamiliar with doing this, I recommend you see my tutorial on reading NetCDF data.

Getting Started - Grabbing GFS Data:
The following code is used to get started with this tutorial, essentially, it imports all of the required modules, pulls the GFS data from the NOMADS data server, and loads it into arrays.  Once the data is loaded, the remainder of the tutorial will focus on Cartopy. 

import numpy as np
from matplotlib import pyplot as plt
from cartopy import crs as ccrs

import cartopy.feature as cfeature

import netCDF4 as nc
import datetime as DT 

datetime=(DT.datetime.utcnow()-DT.timedelta(hours=6)).strftime('%Y%m%d')
fhour=20

nomads_path="https://nomads.ncep.noaa.gov:9090/dods/gfs_0p25/gfs{}/gfs_0p25_00z".format(datetime)

datafile=nc.Dataset(nomads_path)
lat=datafile.variables['lat'][:]
lon=datafile.variables['lon'][:] 

sfcT=datafile.variables['tmpsfc'] #Keep meta data for now!
sfcT_data=sfcT[fhour,:]

title=sfcT.long_name
time=datafile.variables['time']
tstart=DT.datetime.strptime(time.minimum,'%Hz%d%b%Y')
timestr=(tstart+DT.timedelta(hours=time.resolution*24.*fhour)).strftime('%c')

Now that your data is loaded, we initialize a figure, and subplot using matplotlibs "projection" keyword to link to the Cartopy library.  Once you've attached the Cartopy package to your subplot, there are a number of additional objects and classes you can add to your figure.  In the first example, I'll demonstrate by calling upon the "coastlines" attribute to plot continents on the map.

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())

Now, your figure and subplot are defined.  The Cartopy "PlateCarree()" projection is a basic cylindrical map projection that accepts coordinate information in the form of lat/lon coordinate pairs.  The lat/lon data can be either 1D or 2D, but essentially, if your coordinate information is lat/lon, your projection is PlateCarree().   You can pass 2 keywords into PlateCarree (central_longitude and globe), but in most instances, you won't need to.   Occasionally, I reset the central_longitude to 180, instead of it's default zero.

The code to plot the map is:

Z=plt.pcolormesh(lon,lat,sfcT_data,cmap='jet')

pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines()

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

And you'll get an image similar to this one:

GFS surface temperature from Cartopy

Now, lets say you wanted to mask out land areas, and add sea-ice to the map, to focus on the oceans.
Adding the sea-ice is easy, all it requires is for you to pull the icecsfc variable from the GFS data, and plot it overtop of the surface temperature map:

plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)

We use the numpy mask (ma) functionality to maskout gridcells where the ice cover is less than 10%.
To mask out land areas, we make use of the "Cartopy Feature" module, which handles shapefiles.  The Cartopy Feature module connects seemlessly to the Natural Earth Dataset, and allows for anyone to plot a host of different GIS datasets without needing to download individual shapefiles before hand.  The interface even has many pre-defined functions, including land.  This makes masking out the land areas, a breeze.

ax.add_feature(cfeature.LAND,zorder=4,color='gray')

In the above code, the "add_feature" function is used, and the pre-defined "LAND" feature from the Cartopy Feature package is supplied with a couple of basic keyword arguments.  Putting it all together:

fig=plt.figure(figsize=(11,5))
ax=plt.subplot(111,projection=ccrs.PlateCarree())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.)
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3)
ax.add_feature(cfeature.LAND,zorder=4,color='gray')
pos=ax.get_position()

plt.title(title,fontweight='bold',loc='left')
plt.title("Valid:{} UTC".format(timestr),loc='right')
ax.coastlines(color='k',zorder=5)

cbar_ax=fig.add_axes([pos.x1+0.01,pos.y0,0.015,pos.height])
cbar=plt.colorbar(Z,cax=cbar_ax)
cbar.set_label(title)

datafile.close()

plt.show()

Running the above code will give you an image like this:

GFS SST and Sea-ice cover

Now, lets zoom in over the United States, and add states to the map.  To clip the extent, the "set_extent" function is used:

ax.set_extent([-125, -65, 25, 55],crs=ccrs.PlateCarree())

The list corresponds to lon/lat boundaries from [western-most longitude, eastern-most longitude, southern-most latitude, northern-most latitude] and the crs argument indicates which transform is being applied.  To add the states, you need to define the states from Natural Earth database. the "scale" argument

states = cfeature.NaturalEarthFeature(category='cultural', scale='50m', facecolor='none',
                         name='admin_1_states_provinces_shp',edgecolor='k')

 By changing the extent, and adding the states, you'll get a map that looks like this:

GFS SST centered over the United States

You'll notice that the "LAND" feature masks out the Great Lakes, which is unfortunate, and to my knowledge the only way to get both the land mask and the SST over lakes using GFS surface temperature data, is to use a Lakes shapefile with shapefile clipping, a technique that is beyond this simple introduction, but may be the subject of a future tutorial.

Finally, if you want to use a different map projection to plot the data, e.g., a North Polar Sterographic projection, you simply define a different projection when setting up your subplot, and use the "transform" keyword in the plotting function to specify that you are transforming your lon/lat coordinates from the PlateCarree() projection:

ax=plt.subplot(111,projection=ccrs.NorthPolarStereo())
Z=plt.pcolormesh(lon,lat,sfcT_data-273.15,cmap='jet',vmin=-1,vmax=30.,transform=ccrs.PlateCarree())
plt.pcolormesh(lon,lat,np.ma.masked_less(ice,0.1),vmin=0,vmax=1,cmap='Greys_r',zorder=3,transform=ccrs.PlateCarree())

ax.set_extent([-180, 180, 45, 90],crs=ccrs.PlateCarree()) 

Note, that we also reset the extent to focus on the Northern Latitudes.  You should get something that looks like this:

GFS Sea-ice and SST Polar Plot

That is the extent of this introductory tutorial on how to use Cartopy to work with atmospheric data.  You can learn more from the Cartopy Website, and get a list of different map projections here Cartopy Map Projections. Future tutorials on Cartopy will focus on working with shapefiles, transforming individual points, how to work with data without lat/lon coordinates, and how to put fine details on your map.