Geoscience visualization with Matplotlib and CartoPy

Sophisticated 2D and 3D plotting library.
We'll give you a taste of the possibilities, but covering < 1%.
Mature and can produce publication quality graphics.
Static images, nothing interactive yet.
We'll focus on geoscience, but matplotlib can likely meet most of your scientific plotting requirements.
The best is to find an example of what you want in a gallery and emulate that.

Some CartoPy examples to whet your appetite

Reference: http://scitools.org.uk/cartopy/docs/latest/gallery.html

Plotting netCDF data

Basic Plot

Plotting temperature as a function of depth as predicted by the RTOFS model



In [ ]:

    
# Importing libraries we will need.
import netCDF4
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np



In [ ]:

    
# Open the netCDF file
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')

# Getting the n-dimensional data
tempv = f.variables['temperature']
depth = f.variables['Depth']

print("The temperature variable...\n")
# Note the temperature variable has time, depth, y and x dimensions
print(tempv)
print("The dimensions...\n")
print(tempv.dimensions)



In [ ]:

    
# Continued from previous cell..

# Our goal is temperature as a function of depth so slicing along the depth axis
# at a specific time and place
temp = tempv[0, :, 123, 486]

# Masked arrays are arrays that have bad values idenitifed by the mask array.
print("The masked array containing the temperature data...")
print(repr(temp))

# trick for filtering out good values
x = temp[~temp.mask] 
y = depth[~temp.mask]

print("Plotting...")
# plot and show data
plt.plot(y, x)

# close netCDF
f.close()

Let's build upon the previous plot into something that is ready for publication.

Title
Axis Legends
Markers

Peruse matplotlib gallery and see and emulate what you like.



In [ ]:

    
# Adding text, adjusting borders, and figure size
# Get figure hook to manipulate our plot
fig = plt.figure()
desc ='Figure 1. Temperature as a function of ocean depth as\n'\
    'predicted by the RTOFS model'
# Adding our description
plt.figtext(.5,.15,desc,fontsize=12,ha='center')
#adjust margin
fig.subplots_adjust(bottom=0.35)
#adjust figure size
fig.set_size_inches(7,7)

# Improve axes
# Get axis hook to manipulate our plot
ax = fig.add_subplot(111)
# Add axis labels
ax.set_xlabel('Depth (m)', fontweight='bold')
ax.set_ylabel('Temperature ($^\circ$C)', fontweight='bold')
# Don't show top and right axis
ax.spines["right"].set_visible(False)
ax.spines["top"].set_visible(False)
# Define ticks
ax.tick_params(axis='both', direction='out')
ax.get_xaxis().tick_bottom() # remove unneeded ticks 
ax.get_yaxis().tick_left()

# Getting the data as we did before
f = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = f.variables['temperature']
depth = f.variables['Depth']
temp = tempv[0,:,123,486]
x = temp[~temp.mask] #trick for getting data
y = depth[~temp.mask]

# Plotting line with triangle markers, and red line color.
plt.plot(y, x, marker=r'^', color='r', markersize=10, clip_on=False, linewidth=2.0)
plt.show()
f.close()

Exercise

Create a new notebook cell here.
Based on cell above, plot salinity as a function of depth.
Try using a different marker.
Try using a different line color.

CartoPy

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 [ ]:

    
# Importing CartoPy
import cartopy



In [ ]:

    
# Works with matplotlib's built-in transform support.
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Mercator())

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

# Adds standard background map
ax.stock_img()

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



In [ ]:

    
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())

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

# Can also supply matplotlib kwargs
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')
ax.add_feature(cartopy.feature.LAKES, alpha=0.5)
ax.add_feature(cartopy.feature.RIVERS)

# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])

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



In [ ]:

    
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.PlateCarree())

ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.add_feature(cartopy.feature.COASTLINE)
ax.add_feature(cartopy.feature.BORDERS, linestyle=':')

# Grab state borders
state_borders = cartopy.feature.NaturalEarthFeature(category='cultural', name='admin_1_states_provinces_lines', scale='50m',
                                                    facecolor='none') 

ax.add_feature(state_borders, linestyle="--", edgecolor='blue')

# Set limits in lat/lon space
ax.set_extent([-140, -60, 25, 50])

Plotting data works with the set projection by specifying the transform for the data when plotting. CartoPy will take care of any needed interpolation and transformation. The Plate Carree projection (Equi-rectangular projection) specifies the data in lat/lon, but not on the earth sphere. The Geodetic projection uses lat/lon, but are considered on the sphere. This example shows the difference.



In [ ]:

    
fig = plt.figure(figsize=(12,8))
ax = fig.add_subplot(1, 1, 1, projection=cartopy.crs.Robinson())
ax.set_global()
ax.stock_img()
ax.coastlines()

# Plot circle with and without specifying the transform
ax.plot(-105.295175, 40.013176, 'o')
ax.plot(-105.295175, 40.013176, 'ro', transform=cartopy.crs.Geodetic())

# Plot line between two lat/lon points. One using equi-rectangular and one geodetic
plt.plot([-0.08, 132.], [51.53, 43.17], 'r', transform=cartopy.crs.PlateCarree())
plt.plot([-0.08, 132.], [51.53, 43.17], 'g', transform=cartopy.crs.Geodetic())

Exercise

Try it yourself: Let's look at those curved flight paths on a map

Pick any projection
Plot a map of the United States with borders, coastlines, and state borders.
Plot lines between Seattle, WA and Orlando, FL
- Plot one straight line on the map (assuming lat/lon are rectangular coordinates)
- Plot same points assuming they come from a sphere

Let's plot some data from RTOFS model on a map

Data are coming from the local file system



In [ ]:

    
# Open and read netCDF variables
nc = netCDF4.Dataset('data/rtofs_glo_3dz_f006_6hrly_reg3.nc', 'r')
tempv = nc.variables['temperature']
lat = nc.variables['Latitude'][:]
lon = nc.variables['Longitude'][:]
data = tempv[0, 0, :, :]

# Set up a stereographic projection
proj = cartopy.crs.Stereographic(central_latitude=60., central_longitude=-120.)

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

# Setting the plot size and text
plt.figtext(.5, .15, 'Figure 1. Sea surface temperature as predicted by the RTOFS model', fontsize=12, ha='center')

# define color map
cmap = plt.cm.hsv

# Nice high-level, human-readable abstractions for dealing with maps.
ax.add_feature(cartopy.feature.LAND)
ax.add_feature(cartopy.feature.OCEAN)
ax.coastlines(zorder=3)
ax.gridlines()

#Color-filled contour plot
cs = ax.contourf(lon, lat, data, 50, cmap=cmap, transform=cartopy.crs.PlateCarree(), zorder=2)

# Color bar
cbar = plt.colorbar(cs)
cbar.set_label('$^{o}C$')

nc.close()

Exercise

Create a new notebook cell here.
Based on cell above, create a plot showing salinity on map.
Try a different projection.
Try a different color map (examples).

Accessing and plotting data from the TDS

SST and ICE



In [ ]:

    
from datetime import datetime
date = datetime(2014, 9, 15, 0) # date to plot.

# open dataset.
dataset = netCDF4.Dataset('http://www.ncdc.noaa.gov/thredds/dodsC/OISST-V2-AVHRR_agg')
timevar = dataset.variables['time']
timeindex = netCDF4.date2index(date, timevar) # find time index for desired date.

# read sst.  Will automatically create a masked array using
# missing_value variable attribute. 'squeeze out' singleton dimensions.
sst = dataset.variables['sst'][timeindex, :].squeeze()

# read lats and lons (representing centers of grid boxes).
lats = dataset.variables['lat'][:]
lons = dataset.variables['lon'][:]

# create figure, axes instances.
fig = plt.figure(figsize=(20, 10))
proj = cartopy.crs.Robinson(central_longitude=-105.)
ax = fig.add_axes([0.05, 0.05, 0.9, 0.9], projection=proj)

# plot sst, then ice with imshow--this relies on evenly spaced points
bounds = (lons[0], lons[-1], lats[0], lats[-1])
im1 = ax.imshow(sst, extent=bounds, interpolation='nearest',
                cmap=plt.cm.jet, transform=cartopy.crs.PlateCarree())

# Add gridlines and coastlines
ax.gridlines()
ax.coastlines('50m')

# add colorbar
plt.colorbar(im1)

# add a title.
ax.set_title('Analysis for %s' % date)

Exercise

Modify the last example to add the sea ice analysis to the plot
Try a different projection

Winds and Pressure



In [ ]:

    
# specify date to plot.
yyyy=1993; mm=3; dd=14; hh=0

# set OpenDAP server URL.
URLbase = "http://nomads.ncdc.noaa.gov/thredds/dodsC/modeldata/cmd_pgbh/"
URL = URLbase + "{year}/{year}{month:02}/{year}{month:02}{day:02}/pgbh00.gdas.{year}{month:02}{day:02}{hour:02}.grb2".format(
    year=yyyy, month=mm, day=dd, hour=hh)
data = netCDF4.Dataset(URL)

# read lats,lons
latitudes = data.variables['lat'][:]
longitudes = data.variables['lon'][:]

# get sea level pressure and 10-m wind data.
# mult slp by 0.01 to put in units of hPa.
slpin = 0.01 * data.variables['Pressure_msl'][:].squeeze()
uin = data.variables['U-component_of_wind_height_above_ground'][:].squeeze()
vin = data.variables['V-component_of_wind_height_above_ground'][:].squeeze()

# make 2-d grid of lons, lats
lons, lats = np.meshgrid(longitudes, latitudes)

# make orthographic basemap.
proj = cartopy.crs.Orthographic(central_longitude=-60., central_latitude=60.)

# create figure, add axes
fig = plt.figure(figsize=(20, 10))
ax = fig.add_axes([0.1, 0.1, 0.8, 0.8], projection=proj)
ax.set_global()

# set desired contour levels.
clevs = np.arange(960., 1061., 8.)

# plot SLP contours.
ax.contour(lons, lats, slpin, clevs, linewidths=0.5, colors='k', transform=cartopy.crs.PlateCarree())
ax.contourf(lons, lats, slpin, clevs, cmap=plt.cm.RdBu_r, transform=cartopy.crs.PlateCarree())

stride = 15
ax.quiver(lons[::stride, ::stride], lats[::stride, ::stride],
          uin[::stride, ::stride], vin[::stride, ::stride],
          transform=cartopy.crs.PlateCarree())

# draw coastlines and gridlines
ax.coastlines(linewidth=1.5)
ax.gridlines()

date = datetime(yyyy, mm, dd, hh)
ax.set_title('SLP and Wind Vectors ' + str(date))