In [1]:
import iris
import iris.plot as iplt
import numpy
import matplotlib.pyplot as plt
import seaborn
from decimal import Decimal
import os, sys
cwd = os.getcwd()
repo_dir = '/'
for directory in cwd.split('/')[1:]:
repo_dir = os.path.join(repo_dir, directory)
if directory == 'ocean-analysis':
break
modules_dir = os.path.join(repo_dir, 'modules')
sys.path.append(modules_dir)
import convenient_universal as uconv
In [2]:
%matplotlib inline
In [3]:
def calc_anomaly(cube):
"""Calculate the anomaly."""
anomaly = cube.copy()
anomaly.data = anomaly.data - anomaly.data[0]
anomaly = anomaly[-1, ::]
anomaly.remove_coord('time')
return anomaly
In [4]:
def inferred_wfo(s_orig, s_new, volume):
"""Calculate the inferred cumulative global total wfo for a given change in soga.
wfo = net water flux into sea water
soga = global mean sea water salinity
v = volume of ocean (m3)
"""
p = 1027 # kg/m3; average density of global ocean - could calculate from rhopoto data
m_globe = volume * p
delta_m = -1 * m_globe * (1 - (s_orig / s_new))
return delta_m
In [7]:
def change_budget_terms(model, experiment, direction):
"""Calculate the freshwater change budget terms."""
so_file = '/g/data/r87/dbi599/DRSv2/CMIP5/%s/%s/yr/ocean/r1i1p1/so/latest/dedrifted/so-vertical-%s-mean_Oyr_%s_%s_r1i1p1_all.nc' %(model, experiment, direction, model, experiment)
wfo_file = '/g/data/r87/dbi599/DRSv2/CMIP5/%s/%s/yr/ocean/r1i1p1/wfo/latest/dedrifted/wfo-%s-sum_Oyr_%s_%s_r1i1p1_cumsum-all.nc' %(model, experiment, direction, model, experiment)
volcello_file = '/g/data/r87/dbi599/DRSv2/CMIP5/%s/historical/fx/ocean/r0i0p0/volcello/latest/volcello-vertical-sum_fx_%s_historical_r0i0p0.nc' %(model, model)
wfo_cube = iris.load_cube(wfo_file)
so_cube = iris.load_cube(so_file)
collapse_coord = 'longitude' if direction == 'zonal' else 'latitude'
volume_cube = iris.load_cube(volcello_file)
volume_zonal_sum = volume_cube.collapsed(collapse_coord, iris.analysis.SUM)
#print('global volume:', "{:.2e} m3".format(volume_zonal_sum.data.sum()))
# Anomaly
wfo_anomaly = calc_anomaly(wfo_cube)
so_anomaly = wfo_anomaly.copy()
so_anomaly.data = inferred_wfo(so_cube[0, ::].data, so_cube[-1, ::].data, volume_zonal_sum.data)
ocean_convergence = so_anomaly - wfo_anomaly
ocean_transport_anomaly = ocean_convergence.copy()
ocean_transport_anomaly.data = numpy.ma.cumsum(-1 * ocean_convergence.data)
return wfo_anomaly, so_anomaly, ocean_transport_anomaly
def clim_budget_terms(model, experiment, direction, clim_bounds=None):
"""Calculate the freshwater climatology budget terms"""
wfo_file = '/g/data/r87/dbi599/DRSv2/CMIP5/%s/%s/yr/ocean/r1i1p1/wfo/latest/wfo-%s-sum_Oyr_%s_%s_r1i1p1_all.nc' %(model, experiment, direction, model, experiment)
wfo_cube = iris.load_cube(wfo_file)
if clim_bounds:
start, end = clim_bounds
wfo_clim = wfo_cube[start:end].collapsed('time', iris.analysis.MEAN)
else:
wfo_clim = wfo_cube.collapsed('time', iris.analysis.MEAN)
transport_clim = wfo_clim.copy()
transport_clim.data = numpy.ma.cumsum(wfo_clim.data)
return wfo_clim, transport_clim
def plot_budget(model, experiment, direction):
"""Plot the freshwater budget"""
assert direction in ['zonal', 'meridional']
keep_coord = 'latitude' if direction == 'zonal' else 'longitude'
fig, axs = plt.subplots(1, 2, figsize=[16, 4])
# Change
wfo_anomaly, so_anomaly, transport_anomaly = change_budget_terms(model, experiment, direction)
print('global wfo cumulative anomaly:', "{:.2e} kg".format(wfo_anomaly.data.sum()))
print('global so cumulative anomaly:', "{:.2e} kg".format(so_anomaly.data.sum()))
plt.sca(axs[0])
plt.axhline(0, color='black', linewidth=0.5)
iplt.plot(wfo_anomaly, color='orange', label='wfo')
iplt.plot(so_anomaly, color='blue', label='wfo')
plt.title('Freshwater change (1850-2005)')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
#ax.yaxis.major.formatter._useMathText = True
plt.xlabel(keep_coord)
plt.ylabel('cumulative anomaly (kg)')
# Climatology
wfo_clim, transport_clim = clim_budget_terms(model, experiment, direction)
plt.sca(axs[1])
plt.axhline(0, color='black', linewidth=0.5)
iplt.plot(wfo_clim, color='orange', label='wfo')
plt.title('wfo climatology')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
#ax.yaxis.major.formatter._useMathText = True
plt.xlabel(keep_coord)
plt.ylabel('total annual flux (kg)')
plt.suptitle(model + ', ' + experiment)
plt.show()
def plot_transport(model, experiment, direction):
"""Plot the freshwater transport"""
assert direction in ['zonal', 'meridional']
keep_coord = 'latitude' if direction == 'zonal' else 'longitude'
positive = 'northward' if direction == 'zonal' else 'eastward'
fig, axs = plt.subplots(1, 2, figsize=[16, 4])
# Change
wfo_anomaly, so_anomaly, transport_anomaly = change_budget_terms(model, experiment, direction)
transport_so_zero = wfo_anomaly.copy()
transport_so_zero.data = numpy.ma.cumsum(wfo_anomaly.data)
plt.sca(axs[0])
plt.axhline(0, color='black', linewidth=0.5)
iplt.plot(transport_anomaly, color='green', label='real salinity data')
iplt.plot(transport_so_zero, color='green', linestyle='--', label='salinity set zero')
title = 'Anomalous %s freshwater transport' %(positive)
plt.title(title)
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
#ax.yaxis.major.formatter._useMathText = True
plt.xlabel(keep_coord)
plt.ylabel('cumulative transport (kg)')
plt.legend()
# Climatology
wfo_clim, transport_clim = clim_budget_terms(model, experiment, direction)
plt.sca(axs[1])
plt.axhline(0, color='black', linewidth=0.5)
iplt.plot(transport_clim, color='green')
plt.title(positive + ' transport climatology (so = zero)')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
#ax.yaxis.major.formatter._useMathText = True
plt.xlabel(keep_coord)
plt.ylabel('total annual %s flux (kg)' %(positive))
plt.suptitle(model + ', ' + experiment)
plt.show()
In [8]:
plot_budget('CSIRO-Mk3-6-0', 'historical', 'zonal')
In [9]:
plot_transport('CSIRO-Mk3-6-0', 'historical', 'zonal')
In [11]:
plot_budget('CSIRO-Mk3-6-0', 'historical', 'meridional')
In [12]:
plot_transport('CSIRO-Mk3-6-0', 'historical', 'meridional')
In [16]:
plot_budget('CSIRO-Mk3-6-0', 'historicalGHG')
In [17]:
plot_transport('CSIRO-Mk3-6-0', 'historicalGHG')
In [8]:
def align_yaxis(ax1, ax2):
"""Align zeros of the two axes, zooming them out by same ratio
Taken from: https://stackoverflow.com/questions/10481990/matplotlib-axis-with-two-scales-shared-origin
"""
axes = (ax1, ax2)
extrema = [ax.get_ylim() for ax in axes]
tops = [extr[1] / (extr[1] - extr[0]) for extr in extrema]
# Ensure that plots (intervals) are ordered bottom to top:
if tops[0] > tops[1]:
axes, extrema, tops = [list(reversed(l)) for l in (axes, extrema, tops)]
# How much would the plot overflow if we kept current zoom levels?
tot_span = tops[1] + 1 - tops[0]
b_new_t = extrema[0][0] + tot_span * (extrema[0][1] - extrema[0][0])
t_new_b = extrema[1][1] - tot_span * (extrema[1][1] - extrema[1][0])
axes[0].set_ylim(extrema[0][0], b_new_t)
axes[1].set_ylim(t_new_b, extrema[1][1])
In [7]:
wfo_anomaly_hist, so_anomaly_hist, transport_anomaly_hist = change_budget_terms('CSIRO-Mk3-6-0', 'historical')
wfo_anomaly_ghg, so_anomaly_ghg, transport_anomaly_ghg = change_budget_terms('CSIRO-Mk3-6-0', 'historicalGHG')
wfo_clim_hist_start, transport_clim_hist_start = clim_budget_terms('CSIRO-Mk3-6-0', 'historical', clim_bounds=[0, 10])
wfo_clim_hist_end, transport_clim_hist_end = clim_budget_terms('CSIRO-Mk3-6-0', 'historical', clim_bounds=[-10, -1])
wfo_clim_ghg_start, transport_clim_ghg_start = clim_budget_terms('CSIRO-Mk3-6-0', 'historicalGHG', clim_bounds=[0, 10])
wfo_clim_ghg_end, transport_clim_ghg_end = clim_budget_terms('CSIRO-Mk3-6-0', 'historicalGHG', clim_bounds=[-10, -1])
In [20]:
fig, ax = plt.subplots(figsize=[10, 6])
plt.axhline(y=0, linewidth=0.3)
iplt.plot(wfo_anomaly_hist, color='black', label='historical')
iplt.plot(wfo_anomaly_ghg, color='red', label='GHG-only')
ax.yaxis.major.formatter._useMathText = True
plt.legend()
ax.set_xlabel('latitude')
ax.set_ylabel('cumulative anomaly (kg)')
plt.title('wfo change (1850-2005)')
ax2 = ax.twinx()
plt.sca(ax2)
iplt.plot(wfo_clim_hist_end, color='grey', linestyle='--', label='climatology')
align_yaxis(ax, ax2)
ax2.set_ylabel('total annual flux (kg)')
ax2.yaxis.major.formatter._useMathText = True
plt.show()
In [17]:
fig = plt.figure(figsize=[10, 6])
plt.hlines(0, -85, 88, linestyles='--')
iplt.plot(so_anomaly_hist, color='black', label='historical')
iplt.plot(so_anomaly_ghg, color='red', label='GHG-only')
plt.title('so change (1850-2005)')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
#ax.yaxis.major.formatter._useMathText = True
plt.xlabel('latitude')
plt.ylabel('cumulative anomaly (kg)')
plt.show()
In [21]:
fig, ax = plt.subplots(figsize=[10, 6])
plt.axhline(y=0, linewidth=0.3)
iplt.plot(transport_anomaly_hist, color='black', label='historical')
iplt.plot(transport_anomaly_ghg, color='red', label='GHG-only')
ax.yaxis.major.formatter._useMathText = True
plt.legend()
ax.set_xlabel('latitude')
ax.set_ylabel('cumulative freshwater flux (kg)')
plt.title('anomalous northword freshwater transport (1850-2005)')
ax2 = ax.twinx()
plt.sca(ax2)
iplt.plot(transport_clim_hist_start, color='grey', linestyle='--', label='climatology')
align_yaxis(ax, ax2)
ax2.set_ylabel('total annual northward flux (kg)')
ax2.yaxis.major.formatter._useMathText = True
plt.show()
In [19]:
fig = plt.figure(figsize=[10, 6])
plt.hlines(0, -85, 88, linestyles='--')
iplt.plot(wfo_clim_hist_end, color='black', label='historical')
iplt.plot(wfo_clim_hist_start, color='black', linestyle='--')
iplt.plot(wfo_clim_ghg_end, color='red', label='historical')
iplt.plot(wfo_clim_ghg_start, color='red', linestyle='--')
plt.title('wfo climatology')
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0), useMathText=True)
plt.xlabel('latitude')
plt.ylabel('total annual flux (kg)')
plt.show()
In [152]:
raw_soga_file = '/g/data/r87/dbi599/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/yr/ocean/r1i1p1/soga/latest/soga_Oyr_CSIRO-Mk3-6-0_historical_r1i1p1_all.nc'
raw_soga_cube = iris.load_cube(raw_soga_file)
In [153]:
raw_vm_so_file = '/g/data/r87/dbi599/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/yr/ocean/r1i1p1/so/latest/so-vertical-mean_Oyr_CSIRO-Mk3-6-0_historical_r1i1p1_185001-185912.nc'
raw_vm_so_cube = iris.load_cube(raw_vm_so_file)
In [154]:
raw_vzm_so_file = '/g/data/r87/dbi599/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/yr/ocean/r1i1p1/so/latest/so-zonal-vertical-mean_Oyr_CSIRO-Mk3-6-0_historical_r1i1p1_all.nc'
raw_vzm_so_cube = iris.load_cube(raw_vzm_so_file)
In [155]:
volume_cube = iris.load_cube('/g/data/r87/dbi599/DRSv2/CMIP5/CSIRO-Mk3-6-0/historical/fx/ocean/r0i0p0/volcello/latest/volcello-vertical-sum_fx_CSIRO-Mk3-6-0_historical_r0i0p0.nc')
In [156]:
volume_vs_weights = uconv.broadcast_array(volume_cube.data, [1, 2], raw_vm_so_cube.shape)
inferred_vm_soga = raw_vm_so_cube.collapsed(['latitude', 'longitude'], iris.analysis.MEAN, weights=volume_vs_weights)
In [157]:
volume_vzs_cube = volume_cube.collapsed('longitude', iris.analysis.SUM)
volume_vzs_weights = uconv.broadcast_array(volume_vzs_cube.data, 1, raw_vzm_so_cube.shape)
inferred_vzm_soga = raw_vzm_so_cube.collapsed('latitude', iris.analysis.MEAN, weights=volume_vzs_weights)
In [167]:
fig, ax = plt.subplots()
iplt.plot(raw_soga_cube[0:30])
#iplt.plot(inferred_vm_soga)
iplt.plot(inferred_vzm_soga[0:30])
ax.yaxis.major.formatter._useOffset = False
plt.show()
In [168]:
fig, ax = plt.subplots()
iplt.plot(raw_soga_cube[0:40])
#iplt.plot(inferred_vm_soga)
iplt.plot(inferred_vzm_soga[0:40])
ax.yaxis.major.formatter._useOffset = False
plt.show()
In [169]:
fig, ax = plt.subplots()
iplt.plot(raw_soga_cube)
iplt.plot(inferred_vm_soga)
iplt.plot(inferred_vzm_soga)
ax.yaxis.major.formatter._useOffset = False
plt.show()
So there doesn't appear to be a problem with how I'm integrating the salinity data (i.e. using volume weighted means). The problem must therefore be with the dedrifting... In partcular, the drfited soga timeseries does not match the global sum of the dedrifted zonal mean salinity timeseries. It should.
Perhaps the crazy values at some point are messing things up.
In [13]:
basin_file = '/g/data/ua6/DRSv3/CMIP5/CSIRO-Mk3-6-0/historical/fx/ocean/r0i0p0/basin/latest/basin_fx_CSIRO-Mk3-6-0_historical_r0i0p0.nc'
In [14]:
basin_cube = iris.load_cube(basin_file)
In [16]:
iplt.pcolormesh(basin_cube)
Out[16]:
In [17]:
basin_cube.data.max()
Out[17]:
In [18]:
basin_cube.data.min()
Out[18]:
In [21]:
lon_basin = basin_cube.collapsed('latitude', iris.analysis.MEDIAN)
In [23]:
iplt.plot(lon_basin)
Out[23]:
In [ ]:
In [ ]: