Learning how to use iris cubes and cartopy plotting.

I'm experimenting with iris cubes and cartopy plotting. Both libraries are part of the SciTools package developed by the UK Met Office:

In [2]:

HTML('<iframe src=http://scitools.org.uk width=700 height=350></iframe>')

Out[2]:

They are great for loading, exploring, analyzing and visualizing gridded data like model outputs. iris supports netCDF, GRIB, and PP formats so far.

Here is a quick example to load and plot data from a dap url.

In [3]:

from oceans.colormaps import cm
from iris.fileformats.netcdf import load_cubes

url = 'http://data.nodc.noaa.gov/thredds/dodsC/woa/WOA09/NetCDFdata/silicate_annual_5deg.nc'
cubes = dict()
for cube in load_cubes([url]):
    cubes.update({cube.long_name: cube})

print(cubes.keys())

[u'Number of Observations', u'Standard Deviation from Statistical Mean', u'Statistical Mean', u'Standard Error of the Statistical Mean']

/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/iris/fileformats/cf.py:1139: UserWarning: NetCDF default loading behaviour currently does not expose variables which define reference surfaces for dimensionless vertical coordinates as independent Cubes. This behaviour is deprecated in favour of automatic promotion to Cubes. To switch to the new behaviour, set iris.FUTURE.netcdf_promote to True.
  warnings.warn(msg)

Now let's perform a horizontal slices at the Statistical Mean variable. .slices return an iterator along the depth dimension, so when we use .next() we extract the first depth level only.

In [4]:

surface = cubes['Statistical Mean'].slices(['latitude', 'longitude']).next()

iris provides its own plotting modules. Let's test the iris.quickplot:

In [5]:

import iris.quickplot as qplt

_ = qplt.contourf(surface, 25, extend='both', cmap=cm.odv)

A nice feature of iris.quickplot is that it "knows" about the cube's metadata. Therefore, we automagically got a title for the figure and a unit label for the colorbar. Let's try the iris.plot for more control over the plot. Here we will plot the data on a cartopy axes.

In [6]:

import numpy as np
import matplotlib.pyplot as plt
import iris.plot as iplt
import cartopy.crs as ccrs

levels = np.arange(0, 15, 1)
fig = plt.figure(figsize=(6, 4))
ax = plt.axes(projection=ccrs.PlateCarree())
cs = iplt.contourf(surface, cmap=cm.odv, levels=levels, extend='both')
ax.stock_img()
ax.coastlines()
gl = ax.gridlines()
dx, dy = 40, 20
_ = ax.set_xticks(range(-180, 180+dx, dx))
_ = ax.set_yticks(range(-90, 90+dy, dy))
cbar = fig.colorbar(cs, extend='both', orientation='horizontal', shrink=0.6)
_ = cbar.ax.set_xlabel(r'$\mu$mol l$^{-1}$')

/home/filipe/miniconda3/envs/Blog/lib/python2.7/site-packages/matplotlib/artist.py:221: MatplotlibDeprecationWarning: This has been deprecated in mpl 1.5, please use the
axes property.  A removal date has not been set.
  warnings.warn(_get_axes_msg, mplDeprecation, stacklevel=1)

We went from 2 to 14 lines of code.

In conclusion, iris.quickplot is definitely the way to go when you want a quick way to explorer the data. However, when you need more customization you will probably do better with iris.plot.

Take a look at the docs, and the following talk (from Scipy 2013) to learn more about these two packages.

In [7]:

from datetime import timedelta
from IPython.display import YouTubeVideo
start=int(timedelta(seconds=34).total_seconds())
YouTubeVideo("MW9wmGsscrs", start=start, autoplay=0, theme="light", color="red")

Out[7]:

In [8]:

HTML(html)

Out[8]:

This post was written as an IPython notebook. It is available for download or as a static html.

python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.

python4oceanographers

Turning ripples into waves

Using iris to read gridded files

Comments