@ClarkGRichards, one of oce developers, pointed out two mistakes I made in the original Water Mass Mixing post. One was the fact that I did not use conservative temperature from TEOS-10 (or at least potential temperature from EOS-10) for the T$-$S plot. That resulted in an inversion at the end of the T$-$S diagram. The other mistake was the fact that I used units in the salinity label (g kg$^{-1}$), meaning Absolute Salinity, but instead I plotted Practical Salinity (no units).
I am re-posting the original notebook here fixing those mistakes.
In [2]:
import numpy as np
import numpy.ma as ma
def mixing(T, S, inds):
"""
Compute the water mass mixing percentage using Mamayev's (1975) mixing
triangle.
Parameters
----------
T : Conservative Temperature
S : Absolute Salinity
inds : 2x3 array with thermohaline indices
[T1 T2 T3
S1 S2 S3]
Returns
-------
m1, m2, m3 : Water mass percentage for masses 1, 2 e 3.
"""
a = np.r_[inds, np.ones((1, 3))]
b = np.c_[T.ravel(), S.ravel(), np.ones(T.shape).ravel()].T
m = np.linalg.solve(a, b)
m1 = m[0].reshape(T.shape)
m2 = m[1].reshape(T.shape)
m3 = m[2].reshape(T.shape)
# Mask values outside mixing triangle.
m1 = ma.masked_outside(ma.masked_invalid(m1), 0, 1)
m2 = ma.masked_outside(ma.masked_invalid(m2), 0, 1)
m3 = ma.masked_outside(ma.masked_invalid(m3), 0, 1)
m1 = 100 * m1
m2 = 100 * m2
m3 = 100 * m3
return m1, m2, m3
In [3]:
import os
from oceans.datasets import woa_subset
lon = -25.5
kw = dict(clim_type='annual', resolution='1deg', levels=slice(0, 40),
llcrnrlat=-70, urcrnrlat=-10, llcrnrlon=lon, urcrnrlon=lon)
if os.path.isfile('./data/data.npz'):
npz = np.load('./data/data.npz')
depth = npz['p']
SP, t, p = npz['SP'], npz['t'], npz['p']
lon, lat = npz['lon'], npz['lat']
else:
salinity = woa_subset(var='salinity', **kw)
salinity = salinity['OA Climatology']['annual'].squeeze()
temperature = woa_subset(var='temperature', **kw)
temperature = temperature['OA Climatology']['annual'].squeeze()
SP = salinity.values.filled(fill_value=np.nan)
t = temperature.values.filled(fill_value=np.nan)
p = temperature.columns.values.astype(float)
lat = temperature.index.values.astype(float)
lon = np.array([lon]*len(lat))
np.savez('./data/data.npz', t=t, SP=SP, lon=lon, lat=lat, p=p)
The T$-$S data we just downloaded above are in-situ temperature and practical salinity respectively.
Before we plot a $\Theta-S$ plot we need to compute conservative temperature and Absolute Salinity, if we choose to follow TEOS-10, or we should compute at least Potential Temperature ($\theta$) from EOS-80.
In [4]:
SP, p, lon, lat = np.broadcast_arrays(SP, p[None, :], lon[:, None], lat[:, None])
In [5]:
import gsw
SA = gsw.SA_from_SP(SP, p, lon, lat)
CT = gsw.CT_from_t(SA, t, p)
In [6]:
cores = np.array([[26.42, 5, 4.05], [37.14, 33.90, 35.07]])
sacw, aaiw, nadw = mixing(CT, SA, cores)
In [7]:
s = ma.masked_invalid(SA).mean(axis=0)
t = ma.masked_invalid(CT).mean(axis=0)
Te = np.linspace(t.min(), t.max(), 25)
Se = np.linspace(s.min(), s.max(), 25)
Tg, Sg = np.meshgrid(Te, Se)
sigma_theta = gsw.sigma0(Sg, Tg)
cnt = np.linspace(sigma_theta.min(), sigma_theta.max(), 10)
In [8]:
import brewer2mpl
Reds = brewer2mpl.get_map('Reds', 'Sequential', 9).mpl_colormap
Blues = brewer2mpl.get_map('Blues', 'Sequential', 9).mpl_colormap
Greens = brewer2mpl.get_map('Greens', 'Sequential', 9).mpl_colormap
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator
from mpl_toolkits.basemap import Basemap
from mpl_toolkits.axes_grid1 import make_axes_locatable
from mpl_toolkits.axes_grid1.inset_locator import inset_axes
In [9]:
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(14, 6), facecolor='w')
ax.invert_yaxis()
divider = make_axes_locatable(ax)
# Main figure.
percent = [50, 60, 70, 80, 90, 100]
m1 = ax.contourf(lat, p, sacw, percent, cmap=Reds, zorder=3)
m2 = ax.contourf(lat, p, aaiw, percent, cmap=Greens, zorder=2)
m3 = ax.contourf(lat, p, nadw, percent, cmap=Blues, zorder=1)
m1.set_clim(percent[0], percent[-1])
m2.set_clim(percent[0], percent[-1])
m3.set_clim(percent[0], percent[-1])
# Colorbars.
dy = 0.04
left, bottom, height, width = 0.14, 0.3, 0.13, 0.02
rect1 = [left, bottom, height, width] # Top.
rect2 = [left, bottom - dy, height, width] # Center.
rect3 = [left, bottom - 2*dy, height, width] # Bottom.
cax1 = plt.axes(rect1, axisbg='w')
cax2 = plt.axes(rect2, axisbg='w')
cax3 = plt.axes(rect3, axisbg='w')
kw = dict(orientation='horizontal', extend='min')
cb1 = fig.colorbar(m1, cax=cax1, **kw)
cb2 = fig.colorbar(m2, cax=cax2, **kw)
cb3 = fig.colorbar(m3, cax=cax3, **kw)
cax1.xaxis.set_ticklabels([])
cax2.xaxis.set_ticklabels([])
new_labels = ['%s%%' % l.get_text() for l in cax3.xaxis.get_ticklabels()]
cax3.xaxis.set_ticklabels(new_labels)
kw = dict(rotation=0, labelpad=20, y=1,
verticalalignment='center', horizontalalignment='left')
cb1.ax.set_ylabel('SACW', **kw)
cb1.ax.yaxis.set_label_position("right")
cb2.ax.set_ylabel('AAIW', **kw)
cb2.ax.yaxis.set_label_position("right")
cb3.ax.set_ylabel('NADW', **kw)
cb3.ax.yaxis.set_label_position("right")
# Some tweaking.
ax.xaxis.set_ticks_position('top')
ax.xaxis.set_label_position('top')
ax.yaxis.set_ticks_position('left')
ax.yaxis.set_label_position('left')
ax.xaxis.set_tick_params(tickdir='out')
ax.yaxis.set_tick_params(tickdir='out')
ax.set_xlim(right=lat.max())
ax.set_ylim(bottom=1750)
plt.draw()
new_labels = [r'%s $^\circ$S' % l.get_text()[1:] for l in
ax.xaxis.get_ticklabels()]
ax.xaxis.set_ticklabels(new_labels)
new_labels = [r'%s m' % l.get_text() for l in
ax.yaxis.get_ticklabels()]
ax.yaxis.set_ticklabels(new_labels)
ax.set_xlabel(u'WOA09 water masses at longitude %3.1f\u00B0' % lon[0, 0])
axin = inset_axes(ax, width="35%", height="35%", loc=4)
inmap = Basemap(projection='ortho', lon_0=lon[0, 0], lat_0=0,
ax=axin, anchor='NE')
inmap.fillcontinents()
inmap.plot(lon, lat, 'r', latlon=True)
# T-S Diagram.
axTS = divider.append_axes("right", 2, pad=0.1)
cs = axTS.contour(Sg, Tg, sigma_theta, colors='grey', levels=cnt, zorder=1)
kw = dict(color='r', fontsize=14, fontweight='black')
axTS.text(35.00, 11.9, 'SACW', **kw)
axTS.text(34.50, 3.12, 'AAIW', **kw)
axTS.text(34.82, 1.72, 'NADW', **kw)
axTS.plot(s, t, 'k')
axTS.yaxis.set_label_position("right")
axTS.yaxis.set_ticks_position("right")
axTS.set_xlabel("Absolute Salinity [g kg$^{-1}$]")
axTS.set_ylabel("Conservative Temperature [$^\circ$C]", rotation=-90, labelpad=20)
axTS.set_title("$\Theta$$-$S Diagram")
axTS.xaxis.set_major_locator(MaxNLocator(nbins=5))
In [10]:
HTML(html)
Out[10]:
This post was written as an IPython notebook. It is available for download or as a static html.
python4oceanographers by Filipe Fernandes is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
Based on a work at https://ocefpaf.github.io/.