https://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/data/tables/onetime/1tim_atl.htm#A16
def fix_column(df):
new_column = []
for col, unit in zip(df.columns, df.iloc[0].fillna('')):
if unit:
new_column.append(f'{col.lower()} ({unit.lower()})')
else:
new_column.append(f'{col.lower()}')
df.columns = new_column
df.drop(labels=0, inplace=True)
return df
import pandas as pd
dfs = []
for section in ['c', 'n', 's']:
url = (
f'https://www.nodc.noaa.gov/woce/woce_v3/wocedata_1/whp/data/onetime/atlantic'
f'/a16/a16{section}/a16{section}_hy1.csv'
)
df = pd.read_csv(
url,
skiprows=5,
skipfooter=1,
skipinitialspace=True,
header=0,
engine='python',
)
dfs.append(fix_column(df))
merged = pd.concat(dfs, axis=0, sort=False)
import numpy as np
from datetime import datetime
# Select a continous track.
maskN = np.logical_and(merged['stnnbr'] >= 1, merged['stnnbr'] <= 124)
maskS = np.logical_and(merged['stnnbr'] >= 237, merged['stnnbr'] <= 278)
maskC = np.logical_and(merged['stnnbr'] >= 314, merged['stnnbr'] <= 373)
mask = maskN | maskS | maskC
A16 = merged[mask].copy()
# Use time as index.
dtime = A16['date'].astype(int).astype(str) + ' ' + A16['time'].astype(int).astype(str).str.zfill(4)
dtime = pd.to_datetime(dtime, format='%Y%m%d %H%M')
# Drop unused and non-numeric columns.
A16.drop(['date', 'time', 'expocode', 'sect_id'], axis=1, inplace=True)
A16.index = dtime
# Force all misinterpret object values to be numeric.
for col in A16.select_dtypes('O').columns:
A16[col] = pd.to_numeric(A16[col])
%matplotlib inline
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
projection = {'projection': ccrs.PlateCarree()}
fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=projection)
ax.plot(merged['longitude'], merged['latitude'], marker='o', linestyle='none', zorder=0, label='all points')
ax.plot(A16['longitude'], A16['latitude'], marker='.', linestyle='none', zorder=1, label='A16')
ax.set_extent([-100, 40, -80, 70])
ax.coastlines()
ax.legend(loc='upper right', bbox_to_anchor=(1.5, 1.))
<matplotlib.legend.Legend at 0x7f6feeba9550>
# Multiple NA is used. -999.00, -999.9 and etc,
# so we are getting all data that is marked as less than -999.
A16[A16 <= -999] = np.NaN
def apply_flag(flag):
prop = flag.split('_')[0]
mask = A16[flag] == 2.
A16[prop][mask] = np.NaN
flags = [
'btlnbr_flag_w',
'ctdsal_flag_w',
'salnty_flag_w',
'silcat_flag_w',
'oxygen_flag_w',
'nitrat_flag_w',
'nitrit_flag_w',
'phspht_flag_w',
'cfc-11_flag_w',
'cfc-12_flag_w',
'tcarbn_flag_w',
'pco2_flag_w'
]
for flag in flags:
# Not sure if "exchange" format already applied bad flag.
# So I'm just dropping that column.
A16.drop([flag], axis=1, inplace=True)
import seawater as sw
def o2sat(s, pt):
"""
Calculate oxygen concentration at saturation. Molar volume of oxygen
at STP obtained from NIST website on the thermophysical properties of fluid
systems (http://webbook.nist.gov/chemistry/fluid/).
Parameters
----------
s : array_like
Salinity [pss-78]
pt : array_like
Potential Temperature [degC ITS-90]
Returns
-------
osat : array_like
Oxygen concentration at saturation [umol/kg]
References
-----
.. [1] The solubility of nitrogen, oxygen and argon in water and seawater -
Weiss (1970) Deep Sea Research V17(4): 721-735.
"""
Kelvin = 273.15
t = sw.library.T68conv(pt) + Kelvin
# Eqn (4) of Weiss 1970 (the constants are used for units of ml O2/kg).
a = (-177.7888, 255.5907, 146.4813, -22.2040)
b = (-0.037362, 0.016504, -0.0020564)
lnC = (a[0] + a[1] * (100. / t) + a[2] * np.log(t / 100.) + a[3] *
(t / 100.) +
s * (b[0] + b[1] * (t / 100.) + b[2] * (t / 100.) ** 2))
osat = np.exp(lnC) * 1000. / 22.392 # Convert from ml/kg to um/kg.
# The Apparent Oxygen Utilization (AOU) value was obtained by subtracting
# the measured value from the saturation value computed at the potential
# temperature of water and 1 atm total pressure using the following
# expression based on the data of Murray and Riley (1969):
# ln(O2 in µmol/kg) = - 173.9894 + 255.5907(100/TK) + 146.4813 ln(TK/100) -
# 22.2040(TK/100) + Sal [-0.037362 + 0.016504(TK/100) - 0.0020564(TK/100)2],
# where TK is temperature in °K and Sal in the Practical Salinity (SP) scale.
return osat
# Normalize Total DIC.
A16['tcarbn_norm'] = A16['tcarbn (umol/kg)'] * (35 / A16['ctdsal (pss-78)'])
pt = sw.ptmp(A16['ctdsal (pss-78)'], A16['ctdtmp (its-90)'], A16['ctdprs (dbar)'], pr=0)
A16['pden'] = sw.pden(A16['ctdsal (pss-78)'], A16['ctdtmp (its-90)'], A16['ctdprs (dbar)'], pr=0) - 1000
osat = o2sat(A16['ctdsal (pss-78)'], A16['theta (deg c)'])
A16['aou'] = osat - A16['oxygen (umol/kg)']
# Normalize Total DIC.
A16['tcarbn_norm'] = A16['tcarbn (umol/kg)'] * (35 / A16['ctdsal (pss-78)'])
pt = sw.ptmp(A16['ctdsal (pss-78)'], A16['ctdtmp (its-90)'], A16['ctdprs (dbar)'], pr=0)
A16['pden'] = sw.pden(A16['ctdsal (pss-78)'], A16['ctdtmp (its-90)'], A16['ctdprs (dbar)'], pr=0) - 1000
A16.head()
stnnbr | castno | sampno | btlnbr | latitude | longitude | depth | ctdprs (dbar) | ctdtmp (its-90) | ctdsal (pss-78) | ... | tcarbn (umol/kg) | pco2 (uatm) | pco2tmp (deg_c) | theta (deg c) | ctdtmp (ipts-68) | ctdoxy (umol/kg) | theta (ipts-68) | tcarbn_norm | pden | aou | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1989-03-21 06:09:00 | 314.0 | 1.0 | 1.0 | 1.0 | -32.0 | -24.9933 | 4296.0 | 1.6 | 22.1208 | 35.6488 | ... | 2031.3 | 329.3 | 20.0 | 22.1205 | NaN | NaN | NaN | 1994.330805 | 24.676985 | -4.269523 |
1989-03-21 06:09:00 | 314.0 | 1.0 | 2.0 | 2.0 | -32.0 | -24.9933 | 4296.0 | 31.0 | 21.8409 | 35.6592 | ... | 2032.6 | 323.0 | 20.0 | 21.8348 | NaN | NaN | NaN | 1995.025127 | 24.765151 | -5.796456 |
1989-03-21 06:09:00 | 314.0 | 1.0 | 3.0 | 3.0 | -32.0 | -24.9933 | 4296.0 | 61.9 | 16.9177 | 35.6133 | ... | NaN | NaN | NaN | 16.9075 | NaN | NaN | NaN | NaN | 26.007188 | -12.239147 |
1989-03-21 06:09:00 | 314.0 | 1.0 | 71.0 | 71.0 | -32.0 | -24.9933 | 4296.0 | 98.4 | 15.5412 | 35.5082 | ... | 2073.5 | 397.2 | 20.0 | 15.5259 | NaN | NaN | NaN | 2043.823680 | 26.246580 | 13.102335 |
1989-03-21 06:09:00 | 314.0 | 1.0 | 5.0 | 5.0 | -32.0 | -24.9933 | 4296.0 | 153.9 | 14.6904 | 35.4514 | ... | NaN | NaN | NaN | 14.6674 | NaN | NaN | NaN | NaN | 26.392974 | 24.393390 |
5 rows × 28 columns
def lsqfity(X, Y):
"""
Calculate a "MODEL-1" least squares fit.
The line is fit by MINIMIZING the residuals in Y only.
The equation of the line is: Y = my * X + by.
Equations are from Bevington & Robinson (1992)
Data Reduction and Error Analysis for the Physical Sciences, 2nd Ed."
pp: 104, 108-109, 199.
Data are input and output as follows:
my, by, ry, smy, sby = lsqfity(X,Y)
X = x data (vector)
Y = y data (vector)
my = slope
by = y-intercept
ry = correlation coefficient
smy = standard deviation of the slope
sby = standard deviation of the y-intercept
"""
X, Y = list(map(np.asanyarray, (X, Y)))
# Determine the size of the vector.
n = len(X)
# Calculate the sums.
Sx = np.sum(X)
Sy = np.sum(Y)
Sx2 = np.sum(X ** 2)
Sxy = np.sum(X * Y)
Sy2 = np.sum(Y ** 2)
# Calculate re-used expressions.
num = n * Sxy - Sx * Sy
den = n * Sx2 - Sx ** 2
# Calculate my, by, ry, s2, smy and sby.
my = num / den
by = (Sx2 * Sy - Sx * Sxy) / den
ry = num / (np.sqrt(den) * np.sqrt(n * Sy2 - Sy ** 2))
diff = Y - by - my * X
s2 = np.sum(diff * diff) / (n - 2)
smy = np.sqrt(n * s2 / den)
sby = np.sqrt(Sx2 * s2 / den)
return my, by, ry, smy, sby
def lsqfitx(X, Y):
"""
Calculate a "MODEL-1" least squares fit.
The line is fit by MINIMIZING the residuals in X only.
The equation of the line is: Y = mx * X + bx.
Equations are modified from those in Bevington & Robinson (1992)
Data Reduction and Error Analysis for the Physical Sciences, 2nd Ed."
pp: 104, 108-109, 199.
Data are input and output as follows:
mx, bx, rx, smx, sbx = lsqfitx(X, Y)
X = x data (vector)
Y = y data (vector)
mx = slope
bx = y-intercept
rx = correlation coefficient
smx = standard deviation of the slope
sbx = standard deviation of the y-intercept
"""
X, Y = list(map(np.asanyarray, (X, Y)))
# Determine the size of the vector.
n = len(X)
# Calculate the sums.
Sx = np.sum(X)
Sy = np.sum(Y)
Sx2 = np.sum(X ** 2)
Sxy = np.sum(X * Y)
Sy2 = np.sum(Y ** 2)
# Calculate re-used expressions.
num = n * Sxy - Sy * Sx
den = n * Sy2 - Sy ** 2
# Calculate m, a, rx, s2, sm, and sb.
mxi = num / den
a = (Sy2 * Sx - Sy * Sxy) / den
rx = num / (np.sqrt(den) * np.sqrt(n * Sx2 - Sx ** 2))
diff = X - a - mxi * Y
s2 = np.sum(diff * diff) / (n - 2)
sm = np.sqrt(n * s2 / den)
sa = np.sqrt(Sy2 * s2 / den)
# Transpose coefficients
mx = 1 / mxi
bx = -a / mxi
smx = mx * sm / mxi
sbx = np.abs(sa / mxi)
return mx, bx, rx, smx, sbx
def lsqfitgm(X, Y):
"""
Calculate a "MODEL-2" least squares fit.
The SLOPE of the line is determined by calculating the GEOMETRIC MEAN
of the slopes from the regression of Y-on-X and X-on-Y.
The equation of the line is: y = mx + b.
This line is called the GEOMETRIC MEAN or the REDUCED MAJOR AXIS.
See Ricker (1973) Linear regressions in Fishery Research, J. Fish.
Res. Board Can. 30: 409-434, for the derivation of the geometric
mean regression.
Since no statistical treatment exists for the estimation of the
asymmetrical uncertainty limits for the geometric mean slope,
I have used the symmetrical limits for a model I regression
following Ricker's (1973) treatment. For ease of computation,
equations from Bevington and Robinson (1992) "Data Reduction and
Error Analysis for the Physical Sciences, 2nd Ed." pp: 104, and
108-109, were used to calculate the symmetrical limits: sm and sb.
Data are input and output as follows:
m, b, r, sm, sb = lsqfitgm(X,Y)
X = x data (vector)
Y = y data (vector)
m = slope
b = y-intercept
r = correlation coefficient
sm = standard deviation of the slope
sb = standard deviation of the y-intercept
Note that the equation passes through the centroid: (x-mean, y-mean)
"""
X, Y = list(map(np.asanyarray, (X, Y)))
# Determine slope of Y-on-X regression.
my = lsqfity(X, Y)[0]
# Determine slope of X-on-Y regression.
mx = lsqfitx(X, Y)[0]
# Calculate geometric mean slope.
m = np.sqrt(my * mx)
if (my < 0) and (mx < 0):
m = -m
# Determine the size of the vector.
n = len(X)
# Calculate sums and means.
Sx = np.sum(X)
Sy = np.sum(Y)
xbar = Sx / n
ybar = Sy / n
# Calculate geometric mean intercept.
b = ybar - m * xbar
# Calculate more sums.
# Sxy = np.sum(X * Y) # FIXME: Assigned but never used.
Sx2 = np.sum(X ** 2)
# Sy2 = np.sum(Y ** 2) # FIXME: Assigned but never used.
# Calculate re-used expressions.
# num = n * Sxy - Sx * Sy # FIXME: Assigned but never used.
den = n * Sx2 - Sx ** 2
# Calculate r, sm, sb and s2.
r = np.sqrt(my / mx)
if (my < 0) and (mx < 0):
r = -r
diff = Y - b - m * X
s2 = np.sum(diff * diff) / (n - 2)
sm = np.sqrt(n * s2 / den)
sb = np.sqrt(Sx2 * s2 / den)
return m, b, r, sm, sb
def plot_gmfit(x, y, ax):
"""Plots the geometric mean functional regression (neutral regression.)"""
m, b, r, sm, sb = lsqfitgm(x, y)
x = np.array([x.min(), x.max()])
ax.plot(x, m * x + b, 'k-', linewidth=2)
ax.text(0.75, 0.10, r'$r^2 = %2.2f$' % r, transform=ax.transAxes)
# Plot equation.
ax.text(0.05, 0.93, r'$y = %2.2fx\pm%2.2f %+2.2f\pm%2.2f$' %
(m, sm, b, sb), transform=ax.transAxes)
from matplotlib import style
style.use('ggplot')
cbarkw = {
'extend': 'both',
'orientation': 'vertical',
}
sckw = {
'alpha': 0.7,
'vmax': 350,
'vmin': 0,
}
df = A16[['phspht (umol/kg)', 'nitrat (umol/kg)', 'oxygen (umol/kg)']].dropna()
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['phspht (umol/kg)'], df['nitrat (umol/kg)'], c=df['oxygen (umol/kg)'], **sckw)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'Oxygen ($\mu$mol kg$^{-1}$)')
ax.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'Nitrate ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['nitrat (umol/kg)'], ax)
df = A16[['phspht (umol/kg)', 'tcarbn_norm', 'aou']].dropna()
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['phspht (umol/kg)'], df['tcarbn_norm'], c=df['aou'], **sckw)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
ax.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'nTDIC (Sal normalized) ($\mu$mol kq$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['tcarbn_norm'], ax)
# Missing CONVRADCARBNAGE, using CFC ratio instead.
df = A16[['phspht (umol/kg)', 'aou', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
nage = (age + abs(age.min())) / (age.ptp())
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['phspht (umol/kg)'], df['aou'], c=nage, alpha=0.7)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'CFC 11/12')
ax.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['aou'], ax)
df = A16[['nitrat (umol/kg)', 'aou', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
nage = (age + abs(age.min())) / (age.ptp())
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['nitrat (umol/kg)'], df['aou'], c=nage, alpha=0.7)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'CFC 11/12')
ax.set_xlabel(r'Nitrate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['nitrat (umol/kg)'], df['aou'], ax)
df = A16[['tcarbn_norm', 'aou', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
nage = (age + abs(age.min())) / (age.ptp())
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['tcarbn_norm'], df['aou'], c=nage, alpha=0.7)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'CFC 11/12')
ax.set_xlabel(r'nTDIC (Sal normalized) ($\mu$mol kq$^{-1}$)')
ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['tcarbn_norm'], df['aou'], ax)
df = A16[['phspht (umol/kg)', 'aou', 'pden', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
mask = np.logical_and(df['pden'] >=25, df['pden'] <= 26.7)
df = df[mask]
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
nage = (age + abs(age.min())) / (age.ptp())
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['phspht (umol/kg)'], df['aou'], c=nage, alpha=0.7)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'CFC 11/12')
ax.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['aou'], ax)
df = A16[['nitrat (umol/kg)', 'aou', 'pden', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
mask = np.logical_and(df['pden'] >=25, df['pden'] <= 26.7)
df = df[mask]
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
nage = (age + abs(age.min())) / (age.ptp())
fig, ax = plt.subplots(figsize=(9, 9))
cs = ax.scatter(df['nitrat (umol/kg)'], df['aou'], c=nage, alpha=0.7)
cbar = fig.colorbar(cs, **cbarkw)
cbar.ax.set_ylabel(r'CFC 11/12')
ax.set_xlabel(r'Nitrate ($\mu$mol kg$^{-1}$)')
ax.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['nitrat (umol/kg)'], df['aou'], ax)
from palettable.cmocean.sequential import Dense_20
dense = Dense_20.mpl_colormap
df = A16[['phspht (umol/kg)', 'nitrat (umol/kg)', 'oxygen (umol/kg)']].dropna()
fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(
figsize=(15, 13),
ncols=2,
nrows=2,
)
# C:AOU:N:P -> 110.47:66.48:15.3:1
## N:P -> 15.3:1
cs = ax0.scatter(df['phspht (umol/kg)'], df['nitrat (umol/kg)'], c=df['oxygen (umol/kg)'], **sckw)
cbar = fig.colorbar(cs, ax=ax0, **cbarkw)
cbar.ax.set_title('AOU\n' r'($\mu$mol kg$^{-1}$)', fontsize=10)
ax0.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax0.set_ylabel(r'Nitrate ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['nitrat (umol/kg)'], ax0)
## AOU:P -> 66.48:1
df = A16[['phspht (umol/kg)', 'aou', 'pden']].dropna()
cs = ax1.scatter(df['phspht (umol/kg)'], df['aou'], c=df['pden'], alpha=0.7, cmap=dense)
cbar = fig.colorbar(cs, ax=ax1, **cbarkw)
cbar.ax.set_title(r'$\sigma_{\theta}$' '\n' r'(kg m $^{-3}$)', fontsize=10)
ax1.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax1.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['aou'], ax1)
# AOU:N -> 66.48:15.3 (4.34:1)
df = A16[['nitrat (umol/kg)', 'aou', 'pden', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
mask = np.logical_and(df['pden'] >=25, df['pden'] <= 26.7)
df = df[mask]
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
age = (age - abs(age.min())) / (age.ptp())
cs = ax2.scatter(df['nitrat (umol/kg)'], df['aou'], c=age, alpha=0.7, cmap=plt.cm.viridis_r)
cbar = fig.colorbar(cs, ax=ax2, **cbarkw)
cbar.ax.set_title(r'CFC 11/12', fontsize=10)
ax2.set_xlabel(r'Nitrate ($\mu$mol kg$^{-1}$)')
ax2.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['nitrat (umol/kg)'], df['aou'], ax2)
## AOU:P -> 66.48:1
df = A16[['phspht (umol/kg)', 'aou', 'pden', 'cfc-11 (pmol/kg)', 'cfc-12 (pmol/kg)']].dropna()
mask = np.logical_and(df['pden'] >=25, df['pden'] <= 26.7)
df = df[mask]
age = df['cfc-11 (pmol/kg)'] / df['cfc-12 (pmol/kg)']
age = np.ma.masked_invalid(age)
age = (age - abs(age.min())) / (age.ptp())
cs = ax3.scatter(df['phspht (umol/kg)'], df['aou'], c=age, alpha=0.7, cmap=plt.cm.viridis_r)
cbar = fig.colorbar(cs, ax=ax3, **cbarkw)
cbar.ax.set_title(r'CFC 11/12', fontsize=10)
ax3.set_xlabel(r'Phosphate ($\mu$mol kg$^{-1}$)')
ax3.set_ylabel(r'AOU ($\mu$mol kg$^{-1}$)')
plot_gmfit(df['phspht (umol/kg)'], df['aou'], ax3)
fig.savefig('redfield_ompa.pdf', dpi=90)
These four are where created with the The WOCE Hydrographic Program Office (WHPO) data for the Atlantic section (A16, see Out4 [4] for a map). The data bottle have records of salinity, temperature, oxygen, phosphate, nitrate, and total inorganic carbon in the water column.
The goal of this figure is to teach the concept of the Redfield ratio and how we can identify water masses based on nutrients that are not being re-mineralized as predicted by Redfield.
In the top left panel we can see that the proportion of Nitrate to Phosphate is linear and almost exactly what we expect from Redfield (15.3 to 1 and we got 15.22 to 1). This means that all the nutrients are being re-mineralized and consumed at the expected rate.
However, in the top right panel, we can see a lot of scattered points that makes this deviate from the expected ratio for Apparent Oxygen Utilization and Phosphate (AOU), expected 66.48 to 1 but got 76.1 to 1. This indicates that some nutrients may have been pre-formed before the system reached a balance. This extra nutrient won't be utilized in the biogeochemical sense and can be used as a tracer for water masses.
In order to investigate the scatter we can filter the date based on the density layer. This can help us reduce the number of water masses we are looking at. The bottom panels show the Nitrate/AOU (left) and Phosphate/AOU (right) ratios for the density layer of >=25,<= 26.7 kg m$^{-3}$.
This is a layer closer to the surface, probably where the extra nutrients are coming from. The colors the CFC-11/CFC-12 ratio normalized, that gives an idea of the water mass age. We can note easily that the "outliers" are older and therefore, the pre-formed nutrients we want to use as tracer for this particular water mass.