In [1]:
import iris
filename = iris.sample_data_path("A1B_north_america.nc")
cube = iris.load_cube(filename)
print(cube)
2. Extract just data from the year 1980 and beyond from the loaded cube
In [2]:
tcoord = cube.coord('time')
def since_1980(cell):
return tcoord.units.num2date(cell.point).year >= 1980
tcon = iris.Constraint(time=since_1980)
cube = cube.extract(tcon)
In [4]:
tcoord = cube.coord('time')
print(tcoord.units.num2date(tcoord.points.min()))
print(tcoord.units.num2date(tcoord.points.max()))
3. Define a function which takes a coordinate and a single time point as arguments, and returns the decade. For example, your function should return 2010 for the following:
time = iris.coords.DimCoord([10], 'time', units='days since 2018-01-01')
print your_decade_function(time, time.points[0])
In [5]:
def get_decade(coord, point):
year = coord.units.num2date(point).year
return (year // 10) * 10
time = iris.coords.DimCoord([10], 'time', units='days since 2018-01-01')
print(get_decade(time, time.points[0]))
4. Add a "decade" coordinate to the loaded cube using your function and the coord categorisation module
In [6]:
import iris.coord_categorisation as coord_cat
coord_cat.add_categorised_coord(cube, 'decade', 'time', get_decade)
print(cube.coord('decade'))
5. Calculate the decadal means cube for this scenario
In [7]:
import iris.analysis
cube = cube.aggregated_by('decade', iris.analysis.MEAN)
print(cube)
6. Create a figure with 3 rows and 4 columns displaying the decadal means, with the decade displayed prominently in each axes' title (hint: the slices
or slices_over
method of the cube will be helpful, especially combined with the built-in enumerate
function in a for-loop)
In [9]:
import matplotlib.pyplot as plt
import iris.plot as iplt
plt.figure(figsize=(12, 6))
plt.suptitle('Decadal means for the A1B scenario')
for i, decade_cube in enumerate(cube.slices(['latitude', 'longitude'])):
plt.subplot(3, 4, i+1)
iplt.contourf(decade_cube, 20, cmap='viridis')
plt.title('{}'.format(decade_cube.coord('decade').points[0]))
plt.gca().coastlines()
plt.show()
In [ ]: