In [1]:
import matplotlib.pyplot as plt
import iris
import iris.plot as iplt
import numpy
import iris.coord_categorisation
import re
In [2]:
%matplotlib inline
In [3]:
infile = '/g/data/ua6/DRSv2/CMIP5/CSIRO-Mk3-6-0/rcp85/mon/ocean/r1i1p1/tauuo/latest/tauuo_Omon_CSIRO-Mk3-6-0_rcp85_r1i1p1_200601-210012.nc'
In [4]:
cube = iris.load_cube(infile, 'surface_downward_x_stress')
In [5]:
print(cube)
In [6]:
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
In [7]:
with iris.FUTURE.context(cell_datetime_objects=True):
start_cube = cube.extract(get_time_constraint(['2006-01-01', '2025-12-31']))
start_clim = start_cube.collapsed('time', iris.analysis.MEAN)
In [8]:
fig = plt.figure(figsize=[20,8])
iplt.contourf(start_clim, cmap='RdBu_r',
levels=[-0.25, -0.2, -0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15, 0.20, 0.25],
extend='both')
plt.gca().coastlines()
plt.gca().gridlines()
cbar = plt.colorbar()
cbar.set_label(str(start_clim.units))
plt.show()
In [9]:
with iris.FUTURE.context(cell_datetime_objects=True):
end_cube = cube.extract(get_time_constraint(['2081-01-01', '2100-12-31']))
end_clim = end_cube.collapsed('time', iris.analysis.MEAN)
In [10]:
fig = plt.figure(figsize=[20,8])
iplt.contourf(end_clim, cmap='RdBu_r',
levels=[-0.25, -0.2, -0.15, -0.1, -0.05, 0, 0.05, 0.1, 0.15, 0.20, 0.25],
extend='both')
plt.gca().coastlines()
plt.gca().gridlines()
cbar = plt.colorbar()
cbar.set_label(str(end_clim.units))
plt.show()
In [11]:
iris.coord_categorisation.add_year(cube, 'time')
annual_cube = cube.aggregated_by(['year'], iris.analysis.MEAN)
In [12]:
zonal_mean = annual_cube.collapsed('longitude', iris.analysis.MEAN)
In [13]:
print(zonal_mean)
In [14]:
iplt.plot(zonal_mean[0, :])
plt.show()
In [15]:
y_data = zonal_mean[0, :].data
x_data = zonal_mean.coord('latitude').points
In [16]:
y_data
Out[16]:
In [17]:
x_data
Out[17]:
In [18]:
plt.plot(x_data[130:150], y_data[130:150], 'o-')
Out[18]:
In [19]:
from scipy.interpolate import interp1d
In [20]:
#xnew = numpy.linspace(x_data[0], x_data[-1], num=1000, endpoint=True)
f_data = interp1d(x_data, y_data, kind='cubic')
In [25]:
xnew = numpy.linspace(x_data[0], x_data[-1], num=1000, endpoint=True)
ynew = f_data(xnew)
plt.plot(x_data, y_data, 'o', xnew, ynew, '--')
plt.legend(['data', 'cubic'], loc='best')
plt.xlim(30, 55)
plt.ylim(0, 0.07)
plt.show()
Interpolation still noisy... need to fit a smooth curve instead.
In [28]:
xnew = numpy.linspace(x_data[0], x_data[-1], num=1000, endpoint=True)
ynew = f_data(xnew)
plt.plot(x_data, y_data, 'o', xnew, ynew, '--')
plt.legend(['data', 'cubic'], loc='best')
plt.axhline()
plt.show()
In [29]:
type(ynew)
Out[29]:
In [41]:
def wind_stress_metrics(xdata, ydata, hemisphere, direction):
"""Calculate location and magnitude metric for surface wind stress."""
assert hemisphere in ['nh', 'sh']
assert direction in ['westerly', 'easterly']
# assert regularly spaced y values
if (hemisphere == 'sh') and (direction == 'westerly'):
print(hemisphere, direction)
selection = (xdata < 0) & (ydata > 0)
elif (hemisphere == 'nh') and (direction == 'westerly'):
print(hemisphere, direction)
selection = (xdata > 0) & (ydata > 0)
elif (hemisphere == 'sh') and (direction == 'easterly'):
print(hemisphere, direction)
selection = (xdata < 0) & (xdata > -50) & (ydata < 0)
elif (hemisphere == 'nh') and (direction == 'easterly'):
print(hemisphere, direction)
selection = (xdata > 0) & (xdata < 40) & (ydata < 0)
x_filtered = numpy.where(selection, xdata, 0)
y_filtered = numpy.where(selection, ydata, 0)
location = numpy.sum(x_filtered * y_filtered) / numpy.sum(y_filtered) # centre of gravity
magnitude = numpy.sum(y_filtered)
return location, magnitude
The centre of mass/gravity idea came from this post.
In [46]:
sw_loc, sw_mag = wind_stress_metrics(xnew, ynew, 'sh', 'westerly')
print(sw_loc, sw_mag)
In [47]:
nw_loc, nw_mag = wind_stress_metrics(xnew, ynew, 'nh', 'westerly')
print(nw_loc, nw_mag)
In [48]:
se_loc, se_mag = wind_stress_metrics(xnew, ynew, 'sh', 'easterly')
print(se_loc, se_mag)
In [49]:
ne_loc, ne_mag = wind_stress_metrics(xnew, ynew, 'nh', 'easterly')
print(ne_loc, ne_mag)
In [51]:
xnew = numpy.linspace(x_data[0], x_data[-1], num=1000, endpoint=True)
ynew = f_data(xnew)
plt.plot(x_data, y_data, 'o', xnew, ynew, '--')
plt.legend(['data', 'cubic'], loc='best')
plt.axhline()
plt.axvline(sw_loc)
plt.axvline(nw_loc)
plt.axvline(se_loc)
plt.axvline(ne_loc)
plt.show()
In [ ]: