In [126]:
from __future__ import (absolute_import, division, print_function)
from six.moves import (filter, input, map, range, zip) # noqa
import iris
import cf_units
import iris.plot as iplt
import iris.quickplot as qplt
import matplotlib.pyplot as plt
import numpy
import cartopy.crs as ccrs
import cf_units
from iris.analysis.cartography import cosine_latitude_weights
import math
import seaborn
In [63]:
%matplotlib inline
In [64]:
fname = '/g/data/r87/dbi599/argo_ohc_maps.nc'
cube = iris.load_cube(fname, 'ocean heat content')
In [65]:
print(cube)
In [66]:
print(cube[0])
In [67]:
#cmap = getattr(plt.cm, 'hot_r')
fig = plt.figure(figsize=(15, 10))
plt.subplot(111)
qplt.contourf(cube[3], 15)# cmap=cmap)
plt.gca().coastlines()
# Plot #2: contourf with axes longitude from 0 to 360
#proj = ccrs.PlateCarree(central_longitude=-180.0)
#ax = plt.subplot(122, projection=proj)
#qplt.contourf(temperature, 15)
#plt.gca().coastlines()
iplt.show()
In [68]:
print(cube.shape)
In [69]:
def polyfit(data, time_axis):
"""Fit polynomial to data."""
if data.mask[0]:
return data.fill_value
else:
return numpy.polynomial.polynomial.polyfit(time_axis, data, 1)[-1]
The relationship between Watts and Joules is $W = J s^{-1}$...
In [78]:
cube_J = cube * 10**12
time_axis = cube_J.coord('time')
print(time_axis.points[0:3])
new_unit = cf_units.Unit('seconds since 2004-01-01 00:00:00', calendar='gregorian')
time_axis.convert_units(new_unit)
print(time_axis.points[0:3])
trend = numpy.ma.apply_along_axis(polyfit, 0, cube_J.data, time_axis.points)
trend = numpy.ma.masked_values(trend, cube_J.data.fill_value)
In [79]:
trend.shape
Out[79]:
In [96]:
fig = plt.figure(figsize=[15, 4])
nrows = 1
ncols = 1
plotnum = 1
ax = plt.subplot(nrows, ncols, plotnum, projection=ccrs.PlateCarree(central_longitude=180.0))
x = cube.coord('longitude').points
y = cube.coord('latitude').points
cmap = plt.cm.RdBu_r
#cmap.set_bad('blue')
#cmap.set_over('blue')
#cmap.set_under('blue')
ticks = numpy.arange(-15, 17.5, 2.5)
cf = ax.contourf(x, y, trend, transform=ccrs.PlateCarree(),
cmap=cmap, levels=ticks, extend='both')
plt.gca().coastlines()
cbar = plt.colorbar(cf)
cbar.set_label('$Wm^{-2}$')
plt.show()
In [97]:
print(trend.max())
print(trend.min())
In moving from $W m^{-2}$ to $W m^{-1}$ need to account for the fact that the length of a degree of longitude is a function of latitude. The formula for the length of one degree of longitude is $(\pi / 180) a \cos{\phi}$, where $a$ is the radius of the earth and $\phi$ is the latitude.
In [105]:
weights = cosine_latitude_weights(cube_J)
In [109]:
weights.shape
Out[109]:
In [110]:
iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS
Out[110]:
In [112]:
math.pi
Out[112]:
In [141]:
#cube_J.coord('longitude').guess_bounds()
level_bounds = cube_J.coord('longitude').bounds
print(level_bounds.shape)
In [147]:
lon_diffs = numpy.apply_along_axis(lambda x: x[1] - x[0], 1, level_bounds)
lon_diffs = lon_diffs[numpy.newaxis,:]
In [148]:
print(lon_diffs.shape)
In [149]:
lon_extents = (math.pi / 180.) * iris.analysis.cartography.DEFAULT_SPHERICAL_EARTH_RADIUS * weights[0,:,:] * lon_diffs
In [150]:
lon_extents
Out[150]:
In [151]:
type(trend)
Out[151]:
In [152]:
weighted_trends = trend * lon_extents
In [153]:
zonal_heat_gain = weighted_trends.sum(axis=1) / 10**7
In [154]:
print(zonal_heat_gain.shape)
print(zonal_heat_gain.max())
print(zonal_heat_gain.min())
In [155]:
plt.plot(zonal_heat_gain, y)
plt.xlabel('Heat gain ($10^7$ $W m^{-1}$)')
plt.ylabel('Latitude')
plt.axvline(0, color='0.7', linestyle='solid')
plt.show()
In [ ]: