Playing around with the CERES EBAF-TOA Ed4.0 net TOA all-sky radiative flux data downloaded from here.
In [4]:
import re
import glob
import numpy
import matplotlib.pyplot as plt
import iris
import iris.coord_categorisation
from iris.experimental.equalise_cubes import equalise_attributes
import iris.plot as iplt
import warnings
warnings.filterwarnings('ignore')
In [5]:
%matplotlib inline
In [6]:
infile = '/g/data/r87/dbi599/data_ceres/CERES_EBAF-TOA_Ed4.0_Subset_200003-201312.nc'
In [7]:
cube = iris.load(infile)[0]
cube
Out[7]:
In [8]:
cube.coord('time').guess_bounds()
In [9]:
print(cube.coord('time'))
In [10]:
cube.coord('time')
Out[10]:
In [11]:
def broadcast_array(array, axis_index, shape):
"""Broadcast an array to a target shape.
Args:
array (numpy.ndarray)
axis_index (int or tuple): Postion in the target shape that the
axis/axes of the array corresponds to
e.g. if array corresponds to (depth, lat, lon) in (time, depth, lat, lon)
then axis_index = [1, 3]
e.g. if array corresponds to (lat) in (time, depth, lat, lon)
then axis_index = 2
shape (tuple): shape to broadcast to
For a one dimensional array, make start_axis_index = end_axis_index
"""
if type(axis_index) in [float, int]:
start_axis_index = end_axis_index = axis_index
else:
assert len(axis_index) == 2
start_axis_index, end_axis_index = axis_index
dim = start_axis_index - 1
while dim >= 0:
array = array[numpy.newaxis, ...]
array = numpy.repeat(array, shape[dim], axis=0)
dim = dim - 1
dim = end_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 convert_to_joules(cube):
"""Convert units from Watts to Joules"""
assert 'W' in str(cube.units)
time_span_days = cube.coord('time').bounds[:, 1] - cube.coord('time').bounds[:, 0]
time_span_seconds = time_span_days * 60 * 60 * 24
cube.data = cube.data * broadcast_array(time_span_seconds, 0, cube.shape)
cube.units = str(cube.units).replace('W', 'J')
return cube
def convert_to_annual(cube, full_months=False):
"""Convert data to annual timescale.
Args:
cube (iris.cube.Cube)
full_months(bool): only include years with data for all 12 months
"""
iris.coord_categorisation.add_year(cube, 'time')
iris.coord_categorisation.add_month(cube, 'time')
cube = cube.aggregated_by(['year'], iris.analysis.SUM)
if full_months:
cube = cube.extract(iris.Constraint(month='Jan|Feb|Mar|Apr|May|Jun|Jul|Aug|Sep|Oct|Nov|Dec'))
cube.remove_coord('year')
cube.remove_coord('month')
return cube
def remove_m2(cube):
"""Multipy a cube by its area."""
if not cube.coord('latitude').has_bounds():
cube.coord('latitude').guess_bounds()
if not cube.coord('longitude').has_bounds():
cube.coord('longitude').guess_bounds()
area_weights = iris.analysis.cartography.area_weights(cube)
cube.data = cube.data * area_weights
cube.units = str(cube.units).replace('m-2', '')
return cube
def extract_region(cube, lat_bounds):
"""Extract region of interest from a regular lat/lon grid.
Returns subset of the original cube that contains only the latitude range of interest.
"""
southern_lat, northern_lat = lat_bounds
lat_constraint = iris.Constraint(latitude=lambda cell: southern_lat <= cell < northern_lat)
cube = cube.extract(lat_constraint)
return cube
def spatial_sum(cube, lat_bounds=None):
"""Calculate the spatial sum"""
cube = cube.copy()
if lat_bounds:
cube = extract_region(cube, lat_bounds)
cube = cube.collapsed(['latitude', 'longitude'], iris.analysis.SUM)
cube.remove_coord('latitude')
cube.remove_coord('longitude')
return cube
def calc_anomaly(cube, base_index=0, base_method='mean'):
"""Calculate the anomaly."""
if (base_index > 0) and (base_method == 'mean'):
base_value = cube.data[0 : base_index + 1].mean()
else:
base_value = cube.data[base_index]
anomaly = cube.copy()
anomaly.data = anomaly.data - base_value
return anomaly
def calc_cumsum(cube, base_index=0, base_method='mean'):
"""Calculate the cumulative sum."""
if (base_index > 0) and (base_method == 'mean'):
base_value = cube.data[0 : base_index + 1].mean()
else:
base_value = cube.data[base_index]
anomaly = cube.data - base_value
cumsum = cube.copy()
cumsum.data = numpy.cumsum(anomaly, axis=0)
return cumsum
In [12]:
cube = remove_m2(cube)
cube
Out[12]:
In [13]:
cube = convert_to_joules(cube)
cube
Out[13]:
In [14]:
cube = convert_to_annual(cube, full_months=True)
cube
Out[14]:
In [15]:
print(cube.coord('time'))
In [16]:
global_sum = spatial_sum(cube)
nh_sum = spatial_sum(cube, lat_bounds=[0, 90])
sh_sum = spatial_sum(cube, lat_bounds=[-90, 0])
In [17]:
iplt.plot(global_sum, color='black', label='globe')
iplt.plot(nh_sum, color='blue', label='NH')
iplt.plot(sh_sum, 'red', label='SH')
plt.legend()
plt.show()
A cumulative sum is easy to calculate for the globe, because the equilibrium/stationary climate value is zero:
In [27]:
global_cumsum_0 = calc_cumsum(global_sum)
global_cumsum_1 = calc_cumsum(global_sum, base_index=1, base_method=None)
global_cumsum_2 = calc_cumsum(global_sum, base_index=2, base_method=None)
global_cumsum_mean = calc_cumsum(global_sum, base_index=2, base_method='mean')
nh_cumsum_0 = calc_cumsum(nh_sum)
nh_cumsum_1 = calc_cumsum(nh_sum, base_index=1, base_method=None)
nh_cumsum_2 = calc_cumsum(nh_sum, base_index=2, base_method=None)
nh_cumsum_mean = calc_cumsum(nh_sum, base_index=2, base_method='mean')
sh_cumsum_0 = calc_cumsum(sh_sum)
sh_cumsum_1 = calc_cumsum(sh_sum, base_index=1, base_method=None)
sh_cumsum_2 = calc_cumsum(sh_sum, base_index=2, base_method=None)
sh_cumsum_mean = calc_cumsum(sh_sum, base_index=2, base_method='mean')
In [29]:
iplt.plot(global_cumsum_0, color='black', linestyle='--')
iplt.plot(global_cumsum_1, color='black', linestyle='--')
iplt.plot(global_cumsum_2, color='black', linestyle='--')
iplt.plot(global_cumsum_mean, color='black', label='global')
iplt.plot(nh_cumsum_0, color='blue', linestyle='--')
iplt.plot(nh_cumsum_1, color='blue', linestyle='--')
iplt.plot(nh_cumsum_2, color='blue', linestyle='--')
iplt.plot(nh_cumsum_mean, color='blue', label='NH')
iplt.plot(sh_cumsum_0, color='red', linestyle='--')
iplt.plot(sh_cumsum_1, color='red', linestyle='--')
iplt.plot(sh_cumsum_2, color='red', linestyle='--')
iplt.plot(sh_cumsum_mean, color='red', label='SH')
plt.legend(loc=2)
plt.show()
The same is not true for each hemisphere - in a stationary climate one hemisphere might lose heat and the other gain.
For a long and/or not noisy hemispheric timeseries, you could just use the first data point as the stationary climate value (i.e. and calculate the anomaly timeseries relative to that point before calculating the cumulative sum), but clearly the 13-year CERES record is too noisy for that.
Instead, it's probably best to apply a linear line of best fit to compare the two hemispheres...
In [17]:
def linear_trend(cube):
"""Return the linear trend line"""
trend = cube.copy()
x = cube.coord('time').points
y = cube.data
z = numpy.polyfit(x, y, 1)
p = numpy.poly1d(z)
# the line equation:
print('y=%.6fx+(%.6f)' %(z[0],z[1]))
trend.data = p(x)
return trend
In [18]:
global_trend = linear_trend(global_sum)
nh_trend = linear_trend(nh_sum)
sh_trend = linear_trend(sh_sum)
In [20]:
iplt.plot(global_sum, color='black', label='globe')
iplt.plot(global_trend, color='black', linestyle='--')
iplt.plot(nh_sum, color='blue', label='NH')
iplt.plot(nh_trend, color='blue', linestyle='--')
iplt.plot(sh_sum, 'red', label='SH')
iplt.plot(sh_trend, 'red', linestyle='--')
plt.legend(loc=2)
plt.savefig('/g/data/r87/dbi599/figures/ceres_ebaf-toa-ed40_2001-2013.png')
In [41]:
def get_time_constraint(time_list):
"""Get the time constraint used for reading an iris cube."""
if time_list:
start_date, end_date = time_list
date_pattern = '([0-9]{4})-([0-9]{1,2})-([0-9]{1,2})'
assert re.search(date_pattern, start_date)
assert re.search(date_pattern, end_date)
if (start_date == end_date):
year, month, day = start_date.split('-')
time_constraint = iris.Constraint(time=iris.time.PartialDateTime(year=int(year), month=int(month), day=int(day)))
else:
start_year, start_month, start_day = start_date.split('-')
end_year, end_month, end_day = end_date.split('-')
time_constraint = iris.Constraint(time=lambda t: iris.time.PartialDateTime(year=int(start_year), month=int(start_month), day=int(start_day)) <= t.point <= iris.time.PartialDateTime(year=int(end_year), month=int(end_month), day=int(end_day)))
else:
time_constraint = iris.Constraint()
return time_constraint
def get_data_triplet(model, experiment, ensnum):
"""Get a data triplet of interest"""
time_constraint = get_time_constraint(['1861-01-01', '2005-12-31'])
dir_experiment = 'rcp85' if experiment == 'historical-rcp85' else experiment
mydir = '/g/data/r87/dbi599/DRSv2/CMIP5/%s/%s/yr/atmos/r1i1p1/rndt/latest/dedrifted' %(model, dir_experiment)
output = {}
for region in ['sh', 'nh', 'globe']:
file_start = 'rndt-%s-sum' %(region)
files = glob.glob('%s/%s_*_all.nc' %(mydir, file_start))
assert len(files) == 1, '%s/%s_*_all.nc' %(mydir, file_start)
var_name = 'TOA Incoming Net Radiation %s sum' %(region)
print(files[0], var_name)
cube = iris.load_cube(files[0], var_name) # & time_constraint
new_aux_coord = iris.coords.AuxCoord(ensnum, long_name='ensemble_member', units='no_unit')
cube.add_aux_coord(new_aux_coord)
cube.cell_methods = ()
output[region] = cube
return output['sh'], output['nh'], output['globe']
def equalise_time_axes(cube_list):
"""Make all the time axes the same."""
iris.util.unify_time_units(cube_list)
reference_cube = cube_list[0]
new_cube_list = iris.cube.CubeList([])
for cube in cube_list:
assert len(cube.coord('time').points) == len(reference_cube.coord('time').points)
cube.coord('time').points = reference_cube.coord('time').points
cube.coord('time').bounds = reference_cube.coord('time').bounds
cube.coord('time').units = reference_cube.coord('time').units
cube.coord('time').attributes = reference_cube.coord('time').attributes
new_cube_list.append(cube)
return new_cube_list
def ensemble_aggregation(cube_list, operator):
"""Calculate the ensemble mean."""
aggregators = {'mean': iris.analysis.MEAN, 'median': iris.analysis.MEDIAN}
if len(cube_list) > 1:
equalise_attributes(cube_list)
equalise_time_axes(cube_list)
ensemble_cube = cube_list.merge_cube()
ensemble_agg = ensemble_cube.collapsed('ensemble_member', aggregators[operator])
ensemble_spread = ensemble_cube.collapsed('ensemble_member', iris.analysis.PERCENTILE, percent=[25, 75])
else:
ensemble_agg = cube_list[0]
ensemble_spread = None
return ensemble_agg, ensemble_spread
In [64]:
def plot_data(model, start_time_index=0, end_time_index=156):
"""Plot the model data"""
all_models = ['CanESM2', 'CCSM4', 'CSIRO-Mk3-6-0', 'GISS-E2-R', 'NorESM1-M']
if model == 'ensemble':
models = all_models
else:
models = [model]
sh_cube_list = iris.cube.CubeList([])
nh_cube_list = iris.cube.CubeList([])
globe_cube_list = iris.cube.CubeList([])
for ensnum, model in enumerate(models):
sh_cube, nh_cube, globe_cube = get_data_triplet(model, 'historical-rcp85', ensnum)
sh_cube_list.append(sh_cube)
nh_cube_list.append(nh_cube)
globe_cube_list.append(globe_cube)
sh_ensmid, sh_ensvar = ensemble_aggregation(sh_cube_list, 'mean')
nh_ensmid, nh_ensvar = ensemble_aggregation(nh_cube_list, 'mean')
globe_ensmid, globe_ensvar = ensemble_aggregation(globe_cube_list, 'mean')
iplt.plot(globe_ensmid[start_time_index:end_time_index], color='black', label='globe')
iplt.plot(sh_ensmid[start_time_index:end_time_index], color='red', label='SH')
iplt.plot(nh_ensmid[start_time_index:end_time_index], color='blue', label='NH')
print(globe_ensmid[start_time_index:end_time_index].coord('time')[0])
print(globe_ensmid[start_time_index:end_time_index].coord('time')[-1])
sh_trend = linear_trend(sh_ensmid[start_time_index:end_time_index])
nh_trend = linear_trend(nh_ensmid[start_time_index:end_time_index])
globe_trend = linear_trend(globe_ensmid[start_time_index:end_time_index])
iplt.plot(globe_trend, color='black', linestyle='--')
iplt.plot(sh_trend, color='red', linestyle='--')
iplt.plot(nh_trend, color='blue', linestyle='--')
plt.title(model)
plt.legend(loc=4)
plt.show()
In [65]:
plot_data('CanESM2')
In [68]:
plot_data('CanESM2', start_time_index=151, end_time_index=164)
In [71]:
plot_data('CCSM4')
In [70]:
plot_data('CCSM4', start_time_index=151, end_time_index=164)
In [72]:
plot_data('CSIRO-Mk3-6-0')
In [73]:
plot_data('CSIRO-Mk3-6-0', start_time_index=151, end_time_index=164)
In [77]:
plot_data('GISS-E2-R')
In [75]:
plot_data('GISS-E2-R', start_time_index=151, end_time_index=164)
In [78]:
plot_data('NorESM1-M')
In [79]:
plot_data('NorESM1-M', start_time_index=151, end_time_index=164)
I should possibly do this for the entire CMIP5 ensemble, since I'm only looking at the historical experiment.
What are those dips and why do they happen in single forcing experiments?
Why do some of the models not start at 0 for the global value? (And should I correct that...?)
conda install nc-time-axis helped with some previous issues but not this one...
In [44]:
time_constraint = get_time_constraint(['1861-01-01', '2005-12-31'])
In [45]:
infile = '/g/data/r87/dbi599/DRSv2/CMIP5/CCSM4/historicalGHG/yr/ocean/r1i1p1/ohc/latest/dedrifted/ohc-zonal-sum_Oyr_CCSM4_historicalGHG_r1i1p1_all.nc'
In [46]:
test = iris.load_cube(infile, 'ocean heat content')
In [47]:
test.extract(time_constraint)
In [ ]: