In [1]:
%load_ext load_style
%load_style talk.css


Iris and Cartopy

iris and cartopy are developed by the UK Met. Office.

  • iris is a Python package for analysing and visualising meteorological and oceanographic data sets
  • cartopy is a Python package for advanced map generation with a simple matplotlib interface.


In [2]:
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import datetime as dt

In [8]:
%matplotlib inline

In [9]:
matplotlib.rcParams['figure.figsize'] = (10.0, 8.0)

In [10]:
import iris

In [11]:
import cartopy.crs as ccrs

In [12]:
import iris.quickplot as qplt

loading one day of Tropical Rainfall Measurement Mission (TRMM) rainfall (mm)


In [13]:
date = dt.datetime(2014,07,23)

In [17]:
fname = '../data/3B42RT_daily.{}.nc'.format(date.strftime("%Y.%m.%d"))

In [19]:
trmm = iris.load_cube(fname)

In [20]:
print(trmm)


trmm / (unknown)                    (time: 1; latitude: 240; longitude: 420)
     Dimension coordinates:
          time                           x            -               -
          latitude                       -            x               -
          longitude                      -            -               x
     Attributes:
          title: saved netcdf variable

In [21]:
lats = trmm.coord('latitude').points
lons = trmm.coord('longitude').points

In [22]:
lats


Out[22]:
array([-49.875, -49.625, -49.375, -49.125, -48.875, -48.625, -48.375,
       -48.125, -47.875, -47.625, -47.375, -47.125, -46.875, -46.625,
       -46.375, -46.125, -45.875, -45.625, -45.375, -45.125, -44.875,
       -44.625, -44.375, -44.125, -43.875, -43.625, -43.375, -43.125,
       -42.875, -42.625, -42.375, -42.125, -41.875, -41.625, -41.375,
       -41.125, -40.875, -40.625, -40.375, -40.125, -39.875, -39.625,
       -39.375, -39.125, -38.875, -38.625, -38.375, -38.125, -37.875,
       -37.625, -37.375, -37.125, -36.875, -36.625, -36.375, -36.125,
       -35.875, -35.625, -35.375, -35.125, -34.875, -34.625, -34.375,
       -34.125, -33.875, -33.625, -33.375, -33.125, -32.875, -32.625,
       -32.375, -32.125, -31.875, -31.625, -31.375, -31.125, -30.875,
       -30.625, -30.375, -30.125, -29.875, -29.625, -29.375, -29.125,
       -28.875, -28.625, -28.375, -28.125, -27.875, -27.625, -27.375,
       -27.125, -26.875, -26.625, -26.375, -26.125, -25.875, -25.625,
       -25.375, -25.125, -24.875, -24.625, -24.375, -24.125, -23.875,
       -23.625, -23.375, -23.125, -22.875, -22.625, -22.375, -22.125,
       -21.875, -21.625, -21.375, -21.125, -20.875, -20.625, -20.375,
       -20.125, -19.875, -19.625, -19.375, -19.125, -18.875, -18.625,
       -18.375, -18.125, -17.875, -17.625, -17.375, -17.125, -16.875,
       -16.625, -16.375, -16.125, -15.875, -15.625, -15.375, -15.125,
       -14.875, -14.625, -14.375, -14.125, -13.875, -13.625, -13.375,
       -13.125, -12.875, -12.625, -12.375, -12.125, -11.875, -11.625,
       -11.375, -11.125, -10.875, -10.625, -10.375, -10.125,  -9.875,
        -9.625,  -9.375,  -9.125,  -8.875,  -8.625,  -8.375,  -8.125,
        -7.875,  -7.625,  -7.375,  -7.125,  -6.875,  -6.625,  -6.375,
        -6.125,  -5.875,  -5.625,  -5.375,  -5.125,  -4.875,  -4.625,
        -4.375,  -4.125,  -3.875,  -3.625,  -3.375,  -3.125,  -2.875,
        -2.625,  -2.375,  -2.125,  -1.875,  -1.625,  -1.375,  -1.125,
        -0.875,  -0.625,  -0.375,  -0.125,   0.125,   0.375,   0.625,
         0.875,   1.125,   1.375,   1.625,   1.875,   2.125,   2.375,
         2.625,   2.875,   3.125,   3.375,   3.625,   3.875,   4.125,
         4.375,   4.625,   4.875,   5.125,   5.375,   5.625,   5.875,
         6.125,   6.375,   6.625,   6.875,   7.125,   7.375,   7.625,
         7.875,   8.125,   8.375,   8.625,   8.875,   9.125,   9.375,
         9.625,   9.875], dtype=float32)

In [23]:
plt.imshow(trmm[0].data)


Out[23]:
<matplotlib.image.AxesImage at 0x10ade5790>

quick and dirty mapping using the cartopy qplt function


In [24]:
proj = ccrs.PlateCarree(central_longitude=-180.0)

In [26]:
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

ax = plt.axes(projection=proj)
ax.coastlines()

qplt.contourf(trmm[0],  np.arange(5,50,5), cmap=plt.get_cmap('Blues'), extend='max')

plt.show()


if you want more control, use the matplotlib interface

The line below makes accessible some features that you may want to add to your map

see The cartopy Feature interface doc for more information


In [27]:
import cartopy.feature as cfeature

In [28]:
f = plt.figure()

ax = plt.axes(projection=proj)

im = ax.contourf(lons, lats, trmm[0].data, np.arange(5,50,5),
             transform=ccrs.PlateCarree(), cmap=plt.get_cmap('Blues'), extend='max')

ax.coastlines(linewidth=1.5)
ax.gridlines(crs=proj, draw_labels=True)
cb = plt.colorbar(im, orientation='horizontal', pad=0.05)
cb.set_label('TRMM rainfall for {}: mm/day'.format(date.strftime('%Y/%m/%d')), fontsize=14)

ax.add_feature(cfeature.LAND, alpha=0.3);


/Users/nicolasf/anaconda/envs/iris/lib/python2.7/site-packages/cartopy/io/__init__.py:261: DownloadWarning: Downloading: http://naciscdn.org/naturalearth/110m/physical/ne_110m_land.zip
  warnings.warn('Downloading: {}'.format(url), DownloadWarning)

In [19]:
f = plt.figure()

ax = plt.axes(projection=proj)

im = ax.contourf(lons, lats, trmm[0].data, np.arange(5,50,5),
             transform=ccrs.PlateCarree(), cmap=plt.get_cmap('Blues'), extend='max')

ax.coastlines(lw=0.5)
ax.gridlines(crs=proj, draw_labels=True)
cb = plt.colorbar(im, orientation='horizontal', pad=0.05)
cb.set_label('TRMM rainfall for {}: mm/day'.format(date.strftime('%Y/%m/%d')), fontsize=14)

ax.stock_img()


Out[19]:
<matplotlib.image.AxesImage at 0x109cd1790>

In [ ]: