In [1]:
import numpy
import iris
import matplotlib.pyplot as plt
from matplotlib import gridspec
import iris.quickplot as qplt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
In [2]:
infile = '/g/data/ua6/DRSv2/CMIP5/CSIRO-Mk3-6-0/historicalMisc/mon/aerosol/r1i1p4/od550aer/latest/od550aer_aero_CSIRO-Mk3-6-0_historicalMisc_r1i1p4_185001-201212.nc'
In [4]:
cube = iris.load_cube(infile, 'atmosphere_optical_thickness_due_to_ambient_aerosol')
print(cube)
In [5]:
zonal_mean_cube = cube.collapsed(['longitude'], iris.analysis.MEAN)
print(zonal_mean_cube)
In [6]:
def calc_trend(cube, running_mean=True):
"""Calculate linear trend.
A 12-month running mean can first be applied to the data.
"""
coord_names = [coord.name() for coord in cube.dim_coords]
assert coord_names[0] == 'time'
if running_mean:
cube = cube.rolling_window('time', iris.analysis.MEAN, 12)
time_axis = cube.coord('time')
trend = numpy.ma.apply_along_axis(linear_trend, 0, cube.data, time_axis.points)
trend = numpy.ma.masked_values(trend, cube.data.fill_value)
return trend
def linear_trend(data, time_axis):
"""Calculate the linear trend.
polyfit returns [a, b] corresponding to y = a + bx
"""
if data.mask[0]:
return data.fill_value
else:
return numpy.polynomial.polynomial.polyfit(time_axis, data, 1)[-1]
In [7]:
od_3D_trend = calc_trend(cube)
od_2D_trend = calc_trend(zonal_mean_cube)
In [8]:
%matplotlib inline
In [9]:
def plot_2D_trend(trends, lats, gs):
"""Plot the zonally integrated trends (i.e. zonal heat gain)"""
ax = plt.subplot(gs[1])
plt.sca(ax)
ax.plot(trends * 10**6, lats)
ax.set_xlabel('Trend in zonal average ($10^6$)', fontsize='small')
ax.set_ylabel('Latitude', fontsize='small')
ax.axvline(0, color='0.7', linestyle='solid')
ax.set_yticks([-80, -60, -40, -20, 0, 20, 40, 60, 80])
ax.set_ylim([lats[0], lats[-1]])
In [14]:
def plot_3D_trend(trends, lons, lats, gs, pretty=True):
"""Plot the trends."""
ax = plt.subplot(gs[0], projection=ccrs.PlateCarree(central_longitude=180.0))
plt.sca(ax)
cmap = plt.cm.hot_r
ticks = numpy.arange(0, 12, 1)
cf = ax.contourf(lons, lats, trends * 10**6 , transform=ccrs.PlateCarree(),
cmap=cmap, extend='both', levels=ticks)
ax.coastlines()
if pretty:
ax.set_yticks([-80, -60, -40, -20, 0, 20, 40, 60, 80], crs=ccrs.PlateCarree())
ax.set_xticks([0, 60, 120, 180, 240, 300, 360], crs=ccrs.PlateCarree())
lon_formatter = LongitudeFormatter(zero_direction_label=True)
lat_formatter = LatitudeFormatter()
ax.xaxis.set_major_formatter(lon_formatter)
ax.yaxis.set_major_formatter(lat_formatter)
ax.set_xlabel('Longitude', fontsize='small')
ax.set_ylabel('Latitude', fontsize='small')
cbar = plt.colorbar(cf)
cbar.set_label('trend ($10^6$)')
In [15]:
fig = plt.figure(figsize=[15, 3])
gs = gridspec.GridSpec(1, 2, width_ratios=[4, 1])
lons = cube.coord('longitude').points
lats = cube.coord('latitude').points
plot_3D_trend(od_3D_trend, lons, lats, gs)
plot_2D_trend(od_2D_trend, lats, gs)
plt.show()
In [20]:
fig = plt.figure(figsize=[15, 3])
gs = gridspec.GridSpec(1, 1)
lons = cube.coord('longitude').points
lats = cube.coord('latitude').points
plot_3D_trend(od_3D_trend, lons, lats, gs, pretty=False)
plt.savefig('/g/data/r87/dbi599/figures/od550aer_aero_CSIRO-Mk3-6-0_historicalMisc_r1i1p4_1850-2012-trend.png',
bbox_inches='tight', dpi=200)
In [17]:
regions = {'globe60': [-60, 60],
'tropics': [-20, 20],
'northern_extratropics60': [20, 60],
'southern_extratropics60': [-60, -20],
}
region_dict = {'globe60': ('globe (60S - 60N)', 'black', '-'),
'tropics': ('tropics (20S to 20N)', 'purple', '-'),
'northern_extratropics60': ('northern extratropics (20N - 60N)', 'red', '-'),
'southern_extratropics60': ('southern extratropics (60S - 20S)', 'blue', '-'),
}
In [14]:
def broadcast_array(array, axis_index, shape):
"""Broadcast a one dimensional array to a target shape.
Args:
array (numpy.ndarray): One dimensional array
axis_index (int): Postion in the target shape that the array
corresponds to (e.g. if array corresponds to lat in (time, depth
lat, lon) array then index = 2
shape (tuple): shape to broadcast to
"""
dim = axis_index - 1
while dim >= 0:
array = array[numpy.newaxis, ...]
array = numpy.repeat(array, shape[dim], axis=0)
dim = dim - 1
dim = axis_index + 1
while dim < len(shape):
array = array[..., numpy.newaxis]
array = numpy.repeat(array, shape[dim], axis=-1)
dim = dim + 1
return array
def region_mask(cube, south_bound, north_bound):
"""Create mask for excluding points not in region of interest.
False corresponds to points that are not masked.
Args:
cube (iris.cube.Cube): Data cube
south_bound (float): Southern boundary of region of interest
north_bound (float): Northern boundary of region of interest
"""
data_mask = cube.data.mask
dim_coord_names = [coord.name() for coord in cube.dim_coords]
aux_coord_names = [coord.name() for coord in cube.aux_coords]
lat_coord = cube.coord('latitude')
vin_flag = numpy.vectorize(in_flag)
in_region = vin_flag(lat_coord.points, south_bound, north_bound)
if 'latitude' in dim_coord_names:
lat_index = dim_coord_names.index('latitude')
in_region = broadcast_array(in_region, lat_index, cube.shape)
elif 'latitude' in aux_coord_names:
dim_diff = len(dim_coord_names) - len(aux_coord_names)
if dim_diff == 2:
assert 'time' in dim_coord_names[0:2], "Last two axes must be spatial coordinates"
assert 'depth' in dim_coord_names[0:2], "Last two axes must be spatial coordinates"
elif dim_diff == 1:
assert dim_coord_names[0] == 'depth', "Last two axes must be spatial coordinates"
while in_region.ndim < data_mask.ndim:
in_region = in_region[numpy.newaxis, ...]
mask = data_mask + in_region
return mask
def in_flag(lat_value, south_bound, north_bound):
"""Determine if a point is in the region of interest.
Returns false for points that are included (because they don't need
to be masked)
"""
if lat_value < north_bound and lat_value > south_bound:
return False
else:
return True
def regional_average(cube_copy, mask):
"""Calculate the heat content for each timestep."""
cube_copy.data.mask = mask
coord_names = [coord.name() for coord in cube_copy.coords()]
coord_names.remove('time')
ave = cube_copy.collapsed(coord_names, iris.analysis.MEAN)
return ave
In [15]:
regional_averages = {}
for region, bounds in regions.iteritems():
data_mask = region_mask(cube, bounds[0], bounds[-1])
regional_averages[region] = regional_average(cube.copy(), data_mask)
In [22]:
for region in regions.keys():
name, color, style = region_dict[region]
cube = regional_averages[region]
cube = cube.rolling_window('time', iris.analysis.MEAN, 12)
qplt.plot(cube.coord('time'), cube, label=name, color=color, linestyle=style)
plt.legend(loc='best')
plt.title('Aerosol')
plt.ylabel('ambient aerosol optical thickness at 550 nm')
plt.xlabel('year')
In [ ]: