Iris example plots and scripts

Graph from SCM netCDF file

#!python
import iris
import iris.plot as iplt
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator

# Use the variable long_name as a constraint on the load to get just the
# required fields (standard_name also works)
# http://scitools.org.uk/iris/docs/latest/userguide/loading_iris_cubes.html#constrained-loading

t15 = iris.load_cube('out.nc', '1.5m temperature')
tsurf = iris.load_cube('out.nc', 'Surface temperature')
tsoil = iris.load_cube('out.nc', 'Deep soil temperatures')

iplt.plot(t15[:,0], label='1.5m')
iplt.plot(tsurf[:,0], label='Surface')
iplt.plot(tsoil[:,0], label='Soil')
plt.title('ACCESS 1.3 SCM example')
plt.ylabel('Temperature (K)')

ax = plt.gca() # Current axis instance
ax.xaxis_date()
ax.xaxis.set_major_formatter(DateFormatter("%d\n%b"))
ax.xaxis.set_major_locator(DayLocator())
ax.legend(loc='lower right')

plt.show()

A 2D Plot from SCM netCDF file

#!python
import iris
import iris.plot as iplt
import iris.quickplot as qplt
import matplotlib.pyplot as plt
from matplotlib.dates import DateFormatter, DayLocator
import numpy as np

acf = iris.load_cube('out.nc', 'Area cloud fraction')

qplt.contourf(acf[:,:,0,0], np.linspace(0,1,11))

plt.title('Area cloud fraction')
plt.ylabel('Height (km)')

ax = plt.gca() # Current axis instance
ax.xaxis_date()
ax.xaxis.set_major_formatter(DateFormatter("%d\n%b"))
ax.xaxis.set_major_locator(DayLocator())

plt.show()

Using alternate coastlines

Cartopy and basemap (and many other plotting packages) use the GSHHS coastline dataset (http://www.ngdc.noaa.gov/mgg/shorelines/gshhs.html). However this has errors of up to 1 km around Australia (not simply a datum error because it isn't consistent).

GA has GEODATA Coast 100k Australian coastline dataset available in a variety of formats (http://www.ga.gov.au/meta/ANZCW0703006621.html).

With iris it's reasonably simple to modify the code that uses the GSHHS shapefile to use the GA one. Note that using the GSHHS data is very slow at the moment because iris doesn't yet subset the data according to the target plot region. On raijin this example takes about 90 s to plot the GSHHS coast and 5 s for the rest.

With basemap it's necessary to fiddle with the map function and half grid box offsets to get the coordinates right when using pcolormesh. This isn't necessary with iris.

This map shows the Sydney region with the SRTM 90m resolution orography, the GA coastline (black) and the GSHHS fine coastline in red.

Variable grid files

With version 1.5.0 iris can properly read variable grid fieldsfiles. For this example the input is a Sydney region file from Xiao.

import matplotlib.pyplot as plt
import numpy as np
from numpy import ma
import iris
import iris.plot as iplt

orog = iris.load_cube('SYDNEY2v.orog', 'surface_altitude')

orog.data = ma.masked_where(orog.data==0, orog.data)

levels = np.arange(0,1501,50)

plot = iplt.contourf(orog, levels, extend='both')
plt.colorbar(plot)
plt.title("Sydney variable grid orography")

plt.gca().coastlines(resolution='10m')
plt.savefig('sydv.png')

# Save as netcdf file
iris.save(orog,'sydv.nc')

There's an issue with saving variable grid fields to netCDF in the standard 1.5.0 (https://github.com/SciTools/iris/issues/768) but this has been patched in version on raijin. Coordinates in the netCDF file show the variable spacing. E.g. longitude spacing varies from 0.036 degrees at the edge to 0.0135 in the middle (output from ncdump)

 grid_longitude = 144.9011821253, 144.9371821253, 144.9731821253, 
...
    150.058, 150.0715, 150.085, 150.0985, 150.112, 

Attachments

• sydney_coast.png
• sydv.png