6. Some simple applications of Iris

ATMOSPHERE?OCEAN INTERACTIONS :: PYTHON NOTES

6. Some simple applications of Iris

J. S. Wright jswright@tsinghua.

6.1 THE iris MODULE

The iris module is a software package for working with climate data in python. It has been developed by the SciTools initiative at the United Kingdom Met Office. SciTools has provided a detailed tutorial introducing some of the main features of the package as well as a list of instructive examples from climatology, meteorology, and oceanography. If you are using Anaconda or one of its variants, I highly recommend using the conda-forge channel to install iris rather than the scitools channel recommended on the website:

jswright$ conda install -c conda-forge iris This channel also has several other tools that are useful for climate data analysis, and provides a convenient way to manage them. Learn more at the conda-forge github site.

The logic of iris is a bit different from the netCDF4 + numpy approach to reading and working with climate data. These differences have advantages and disadvantages. Among the advantages, iris provides data structures specifically for working with climate data. These data structures and related tools are also used as building blocks in other useful climate-related modules, such as eofs and windspharm. Disadvantages include a less intuitive syntax and a less mature code infrastructure (iris is sometimes slower and occasionally unstable, but is improving on both counts). In these notes we will work through two examples step by step.

6.2 SURFACE CURRENT ANOMALIES IN THE INDONESIAN SEAS The first step in starting to use iris is learning how to read and select data for input. After importing the module, you can load netcdf files using the iris.load() function.

In [1]: import iris In [2]: iris.load('./data/plt_current.nc')

1

You can also load data from other common file types used in climatology and related fields, such as grib2. This loads the data as an iris CubeList:

In [3]: print(ncdf) 0: eastward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50) 1: eastward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50) 2: eastward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50) 3: northward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50) 4: northward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50) 5: northward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50)

A CubeList acts like a standard python list, but it is made up of iris Cubes. Cubes are the basic building blocks of the iris module (a little like arrays in numpy). We can see more detail by printing one of the elements:

In [4]: print(ncdf[1])

eastward_sea_water_velocity / (m s**-1) (depth: 42; lat: 60; lon: 50)

Dimension coordinates:

depth

x

-

-

latitude

-

x

-

longitude

-

-

x

Attributes:

NCO: "4.6.4"

_CoordinateAxes: time depth lat lon

_Fillvalue: 9.96921e+36

average_op_ncl: dim_avg_n over dimension(s): time

interval_write: monthly

offline_operation: time average and spatial interpolation

phase: All

Here we can see that a Cube contains both the data and attributes (metadata) associated with a netcdf variable. You can also access the standard attributes:

In [5]: for vv in ncdf:

...:

print vv.long_name

Zonal Current

Zonal Current Anomaly (El Nino)

Zonal Current Anomaly (La Nina)

Meridional Current

Meridional Current Anomaly (La Nina)

Meridional Current Anomaly (El Nino)

This behavior becomes very important once we need to select and subset the data. This is done in iris using the Constraint construct, as shown in the following example. Two types of constraints are used here, one to select the variables by the standard name and one to select the variables by an attribute named phase, which differentiates the climatological currents from the anomalies during El Ni?o and La Ni?a:

2

?

?

1 import iris

2

3 # load data

4 ncdf = iris.load('./data/plt_current.nc')

5

6 # constraints

7 unme = iris.Constraint(name='eastward_sea_water_velocity')

8 vnme = iris.Constraint(name='northward_sea_water_velocity')

9 clim = iris.AttributeConstraint(phase='All')

10 nino = iris . A t t r i b u t e C o n s t r a i n t ( phase = ' Warm ')

11 nina = iris . A t t r i b u t e C o n s t r a i n t ( phase = ' Cold ')

?

?

Once we have defined the appropriate constraints, we can select the data. This is done via the

CubeList.extract() function:

?

?

1 # extract variables

2 uclim = ncdf.extract(unme & clim)[0]

3 vclim = ncdf.extract(vnme & clim)[0]

4 unino = ncdf.extract(unme & nino)[0]

5 vnino = ncdf.extract(vnme & nino)[0]

6 unina = ncdf.extract(unme & nina)[0]

7 vnina = ncdf.extract(vnme & nina)[0]

?

?

A few of the behaviors of iris Constraints are evident from this example. Among other

things, we can combine constraints using the logical and (&) and or (|). Second, when we

extract from a CubeList, we get a CubeList back. In this case each selected CubeList has

only one element, so we index it with 0 to return a Cube instead. We can also apply one or

more Constraints when loading data to avoid loading large amounts of data that we don't

need:

?

?

1 # alternate loading

2 vars_clim = iris.load('./data/plt_current.nc', constraints=clim)

3 vars_nino = iris.load('./data/plt_current.nc', constraints=nino)

4 vars_nina = iris.load('./data/plt_current.nc', constraints=nina)

?

?

This approach would result in three CubeList objects, each containing one zonal and one

meridional component of the surface currents.

One convenient feature of iris is that it tracks physical units. This is occasionally annoying,

particularly if your file uses units that are not recognized by the CF Conventions (such as gpm for geopotential height). In this case it allows us to convert from units of m s-1 to cm -1 in a

single logical step:

?

?

1 # convert to cm/s

2 for cc in ncdf:

3

cc.convert_units('centimeter -second^-1')

?

?

3

In the standard numpy approach, we would multiply by a constant to do this conversion and

mark it with a comment. Note that arithmetic with iris Cubes also checks and updates units,

as we will see in the second example.

We can also use Constraints for basic data processing:

?

?

1 # extract uppermost 100m and average

2 srfce = iris.Constraint(depth=lambda d: 0 < d < 100)

3 uclim = srfce.extract(uclim).collapsed('depth', iris.analysis.MEAN)

4 vclim = srfce.extract(vclim).collapsed('depth', iris.analysis.MEAN)

5 unino = srfce.extract(unino).collapsed('depth', iris.analysis.MEAN)

6 vnino = srfce.extract(vnino).collapsed('depth', iris.analysis.MEAN)

7 unina = srfce.extract(unina).collapsed('depth', iris.analysis.MEAN)

8 vnina = srfce.extract(vnina).collapsed('depth', iris.analysis.MEAN)

?

?

In this case the constraint selects all depths less than 100 m and calculates the average over

the depth dimension. Note that extract is used in this case as an intrinsic function of the

Constraint rather than as an intrinsic function of a CubeList (or Cube). Functions for data

processing are provided in the submodule iris.analysis, which is often imported as ia.

The iris module also includes some basic plotting tools. These are not particularly robust

-- just a simple set of wrappers around matplotlib with some convenience shortcuts for

working with geospatial data. To create the plots for this example, we need to import several

modules and prepare some data and plot parameters:

?

?

1 import matplotlib.pyplot as plt , seaborn as sns 2 import numpy as np 3 import cartopy.feature as cfeat

4

5 # cartopy continents 6 land_50m = cfeat.NaturalEarthFeature('physical', 'land', '110m',

7

edgecolor =' #666666 ',

8

facecolor='#999999', zorder=10)

9

10 # plots

11 fig = plt . figure ( figsize =(12 , 6) ) 12 fig . sub plo ts _ad ju st ( wspace =0.25)

13

14 # ================================================================

15 # climatology plot

16 #

17 # current speed

18 speed = ia . maths . exponentiate ( uclim * uclim + vclim * vclim , 0.5)

19 # contour parameters

20 clv = np . linspace (0 , 60 , 13)

21 lon = uclim . coord ( ' longitude ') . points

22 lat = vclim . coord ( ' latitude ') . points

23 # streamplot parameters

24 uclim . data . data [ uclim . data . mask ] = np . NaN

25 vclim . data . data [ vclim . data . mask ] = np . NaN 26 lws = 4* speed . data / speed . collapsed (( ' longitude ' , ' latitude ') , ia .

MAX).data

27 lws [ lws < 0.1] = 0.1

?

?

4

Once these are in place, we can create a geoaxes instance, add the continents, and plot the

currents (as streamlines using plt.streamplot) and the current speed (as a filled contour

plot using iris.plot.contourf):

?

?

1 import iris.plot as iplt 2 import cartopy.crs as ccrs 3 from cartopy.mpl.gridliner import LONGITUDE_FORMATTER ,

LATITUDE_FORMATTER

4

5 # plots

6 axa = plt.subplot(131, projection=ccrs.PlateCarree())

7 axa.add_feature(land_50m)

8 cs0 = iplt.contourf(speed , levels=clv , cmap=plt.cm.Blues , extend='

max ')

9 plt.streamplot(lon , lat , uclim.data , vclim.data , density =1.5, color=

'k', linewidth=lws ,

10

transform=ccrs.PlateCarree(), zorder=2)

11 axa . set_xticks ( r a n g e (90 , 141 , 10) , crs = ccrs . PlateCarree () )

12 axa . set_yticks ( r a n g e ( -30 , 31 , 10) , crs = ccrs . PlateCarree () )

13 axa . xaxis . se t _m a jo r _ fo r ma t te r ( LO N GI T UD E _ FO R MA T TE R )

14 ?axa . yaxis . s e t _ m a j o r _ f o r m a t t e r ( L A T I T U D E _ F O R M A T T E R )

?

Note the use of the LONGITUDE_FORMATTER and LATITUDE_FORMATTER to standardize the axis tickmarks into E and S / N. The widths of the streamlines are adjusted to reflect the speed

of the currents relative to the maximum speed. The final element in the climatology plot is a

colorbar:

?

?

1 # color bar for climatology plot

2 with sns.axes_style('white'):

3

axa_left = axa.get_position().bounds[0]

4

axa_wdth = axa.get_position().bounds[2]

5

cax = fig.add_axes([axa_left , 0.12, axa_wdth , 0.02])

6

cb0 = plt.colorbar(cs0 , cax , orientation='horizontal ')

7

cb0.set_ticks ([0, 10, 20, 30, 40, 50, 60])

8 ? cb0.set_label(r'Surface current [cm s$^{-\mathregular {1}}$]')

?

This approach uses the axes.get_position().bounds construct to help with constructing the axis for the colorbar.

After adding two additional panels that show the anomalies in surface currents (and the speed of those anomalies) during El Ni?o and La Ni?a, this script gives us Figure 6.1. See the uploaded code for a full working version of the script (contact me if you need the data file).

6.3 THE ASIAN TROPOPAUSE AEROSOL LAYER

Our second example is a basic descriptive analysis of the aerosol layer that forms near the tropopause during the Asian summer monsoon. This aerosol layer has only emerged within the last 15 years. It's existence results from several special aspects of the Asian monsoon and the surrounding regions: the severe pollution that often exists in the atmospheric boundary

5

 ................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download