O2 corrections applied to delayed-mode data¶

This notebook performs analysis and correction of Aandreaa data from the following C-PROOF glider deployment:

In [56]:
glider_name = 'dfo-bb046'  
deploy_name = f'{glider_name}-20220707'
deploy_prefix = f'./glider/{glider_name}/{deploy_name}/'
filepath = f'deployments/{glider_name}/{deploy_name}/' # having this is important later for functions that auto-load data
openfile = f'{filepath}/{deploy_name}_CTDadjusted.nc'
opengridfile = f'{filepath}/{deploy_name}_grid_CTDadjusted.nc'
deployfile = f'{filepath}/deployment.yml'
 
description = 'Calvert'
initials = 'LT'

# O2 specs:
sensor = 'oxygen_0054'

# For conductivity filter:
accuracy = 0.0003 #accuracy of the sensor is 0.0003 S/m, used as a cutoff on the exclusion criterion

from datetime import date
processing_date = date.today().strftime('%Y%m%d')
processing_protocol = ''
processing_report = f'CTD_{deploy_name}' 
 

# Import module for loading .md files
from IPython.display import Markdown, display
import os

os.chdir(f'/Users/Lauryn/processing/')
In [66]:
import warnings
warnings.filterwarnings('ignore')

import xarray as xr
import numpy as np
import pathlib
import pyglidersensor as pgs
import pyglider.ncprocess as ncprocess
import pyglider.utils as utils
import sensor as sns 

from datetime import datetime, date
%matplotlib ipympl

import matplotlib.gridspec as gridspec

import scipy.stats as stats
import statsmodels.api as sm

import seawater
import gsw

%matplotlib notebook
%matplotlib inline
import matplotlib.pyplot as plt 
import matplotlib.dates as mdates
from matplotlib.dates import DateFormatter
import cmocean
import cartopy.crs as ccrs
import cartopy.feature as cfeature


import pandas as pd

%load_ext autoreload
%autoreload 2

from scipy import signal
import seawater as sw


from scipy.stats import linregress

from scipy import stats
import ast
import xml.etree.ElementTree as ET
from scipy.interpolate import interp1d, griddata
from scipy.signal import find_peaks

%reload_ext autoreload
%matplotlib ipympl
import pyglider.utils as utils
import logging

The autoreload extension is already loaded. To reload it, use:
  %reload_ext autoreload
In [3]:
 %reload_ext autoreload

1.0 Preamble¶

This document describes dissolved oxygen data processing steps applied to delayed mode data. C-PROOF uses both the Aanderaa 4831 and JFE Advantech Rinko AROD-FT optodes requiring slightly different processing approaches because of their different response times. Please see the C-PROOF O2 Processing Report for more details (https://cproof.uvic.ca/gliderdata/deployments/reports/).

The metadata and sensor calibration sheets are available via the deployment page on the C-PROOF website at: https://cproof.uvic.ca/gliderdata/deployments/ for each mission.

In [4]:
# Summarize info for report:
print(f'** {description}:  glider {glider_name}**')
print(f'************')
print(f'* Deployment: {deploy_name}')
print(f'* Sensor: {sensor}')

#######
ds= xr.open_dataset(openfile)
data_optode=ast.literal_eval(ds.attrs["oxygen"])
make = data_optode.get('make')
serial = data_optode.get('serial')
print(f'make:{make}')
print(f'seriel:{serial}')
print(f'')

# print(f'* Protocols are detailed in: {processing_protocol}')
print(f'* Processing steps will be saved in: O2_{deploy_name}.html')
# print(f'* Files will be located in: {deploy_prefix}')
print(f'* Processed by {initials}, Ocean Sciences Division, Fisheries and Oceans Canada')
print(f'* Processing date: {processing_date}')

** Calvert:  glider dfo-bb046**
************
* Deployment: dfo-bb046-20220707
* Sensor: oxygen_0054
make:AROD_FT
seriel:0054

* Processing steps will be saved in: O2_dfo-bb046-20220707.html
* Processed by LT, Ocean Sciences Division, Fisheries and Oceans Canada
* Processing date: 20260715

2.0 Corrections applied to delayed mode data for this mission¶

This document covers the following steps:

  • Compensate for pressure and salinity effects
  • Manually add Oxygen QC variable
  • Convert O2 from umol/L to umol/kg
  • Determine/apply the sensor response time lag correction
  • Determine/apply the optode drift correction

2.1 Compensation for Pressure and Salinity Effects¶

According to the Aanderaa Oxygen Sensor User’s Manual, dissolved oxygen concentrations reported by the sensor (oxygen_concentration, in μmol L$^{-1}$) must be corrected for both pressure and salinity effects.

Pressure Compensation¶

The pressure compensation factor is calculated as:

$$ p_{\mathrm{factor}} = \left(\frac{|P|}{100}\right) \cdot C_p + 1 $$

where

  • $P$ is the pressure from the CTD (dbar),
  • $C_p$ is the pressure correction calibration coefficient provided on the sensor calibration certificate.

Salinity Compensation¶

The salinity compensation factor is given by:

$$ s_{\mathrm{factor}} = \exp\left[ (S - S_{\mathrm{ref}}) \left( B_0 + B_1 T_s + B_2 T_s^2 + B_3 T_s^3 \right) + C_0 \left(S^2 - S_{\mathrm{ref}}^2\right) \right] $$

where

  • $S$ is the fully corrected salinity from the CTD (PSU),
  • $S_{\mathrm{ref}}$ is the internal reference salinity of the Aanderaa sensor,
  • $B_0$, $B_1$, $B_2$, $B_3$, and $C_0$ are empirical coefficients specified in the Aanderaa manual,
  • $T_s$ is the scaled temperature defined as:

$$ T_s = \ln\left(\frac{298.15 - T}{273.15 + T}\right) $$

and

  • $T$ is the in situ temperature from the CTD (°C).

Note: Aanderaa sensors assume an internal reference salinity ($S_{\mathrm{ref}}$). It is critical to confirm the value of $S_{\mathrm{ref}}$ configured for each sensor before applying the salinity correction.

P & S compensated oxygen concentration¶

The fully compensated dissolved oxygen concentration is calculated as:

$$ \mathrm{oxygen}_{\mathrm{comp}} = \mathrm{oxygen}_{\mathrm{uncomp}} \cdot s_{\mathrm{factor}} \cdot p_{\mathrm{factor}} $$

The resulting compensated oxygen concentration (called ox_P_S_corr after applying the function ), $\mathrm{oxygen}_{\mathrm{comp}}$, has units of μmol L$^{-1}$.

Apply function to make corrections¶

In [68]:
ds= xr.open_dataset(openfile)
ds2 , ds2['s_factor'], ds2['p_factor']= sns.O2_PS_compensation(ds)

## grid to plot easily 
ds2_grid = sns.make_gridfiles_adj(ds2, deployfile)
In [69]:
fig = plt.figure(figsize=(20, 8), constrained_layout=True)

gs = gridspec.GridSpec(
    2, 6,
    figure=fig,
    hspace=0.35,
    wspace=0.4
)

# Top row (3 panels)
ax00 = fig.add_subplot(gs[0, 0:2])
ax01 = fig.add_subplot(gs[0, 2:4])
ax02 = fig.add_subplot(gs[0, 4:6])

# Bottom row (2 centered panels)
ax10 = fig.add_subplot(gs[1, 1:3])
ax11 = fig.add_subplot(gs[1, 3:5])

N = np.arange(len(ds2_grid.time))
cmap = 'RdBu_r'

mean_raw = ds2_grid.oxygen_concentration.mean('time')

# ----------------------------
# Raw O2 anomaly
# ----------------------------
# Raw
sc = ax00.pcolormesh(
    N,
    ds2_grid.depth,
    ds2_grid.oxygen_concentration - mean_raw,
    cmap=cmap,
    vmin=-20,
    vmax=20,
    rasterized=True,
)
ax01.set_title('Raw [O2] anomaly',fontsize=14)

ax00.invert_yaxis()

# Salinity correction
sc = ax01.pcolormesh(
    N,
    ds2_grid.depth,
    ds2_grid.ox_P_S_corr / ds2_grid.p_factor - mean_raw,
    cmap=cmap,
    vmin=-20,
    vmax=20,
    rasterized=True,
)
ax01.invert_yaxis()
ax01.set_title('Correcting salinity bias',fontsize=14)



# Salinity + pressure correction
sc = ax02.pcolormesh(
    N,
    ds2_grid.depth,
    ds2_grid.ox_P_S_corr - mean_raw,
    cmap=cmap,
    vmin=-20,
    vmax=20,
    rasterized=True,
)

# Shared colorbar for top row
cbar = fig.colorbar(
    sc,
    ax=[ax00, ax01, ax02],
    location='right',
    shrink=0.95,
    pad=0.02
)
cbar.set_label(r'O$_2$ anomaly ($\mu$mol L$^{-1}$)',fontsize=14)

ax02.set_title('Correcting salinity & pressure bias',fontsize=14)
ax01.set_xlabel('Profile index',fontsize=14)
ax00.set_ylabel('Depth (m)', fontsize=14)
ax02.invert_yaxis()

# ----------------------------
# Oxygen profiles
# ----------------------------
ax10.grid()

ax10.scatter(
    ds2_grid.oxygen_concentration,
    ds2_grid.pressure,
    s=1,
    c='0.6',
    label='Raw',
    rasterized=True,
)

ax10.scatter(
    ds2_grid.ox_P_S_corr,
    ds2_grid.pressure,
    s=1,
    c='tab:red',
    label='Corrected',
    rasterized=True,
)

ax10.legend(markerscale=6)
ax10.set_title('Oxygen profiles',fontsize=14)
ax10.set_xlabel('O$_2$ ($\\mu$mol L$^{-1}$)',fontsize=14)
ax10.set_ylabel('Depth (m)',fontsize=14)
ax10.invert_yaxis()

# ----------------------------
# Salinity–oxygen
# ----------------------------
ax11.grid()

ax11.scatter(
    ds2_grid.oxygen_concentration,
    ds2_grid.salinity,
    s=1,
    c='0.6',
    label='Raw',
    rasterized=True,
)

ax11.scatter(
    ds2_grid.ox_P_S_corr,
    ds2_grid.salinity,
    s=1,
    c='tab:red',
    label='Corrected',
    rasterized=True,
)

ax11.legend(markerscale=6)
ax11.set_title('Salinity–oxygen diagram',fontsize=14)
ax11.set_xlabel('O$_2$ ($\\mu$mol L$^{-1}$)',fontsize=14)
ax11.set_ylabel('Salinity',fontsize=14)

# ----------------------------
# Panel labels
# ----------------------------
axes = [ax00, ax01, ax02, ax10, ax11]
labels = ['a', 'b', 'c', 'd', 'e']

for axi, label in zip(axes, labels):
    axi.text(
        0.02,
        0.98,
        label,
        transform=axi.transAxes,
        fontsize=15,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(facecolor='white', alpha=0.7, edgecolor='none', pad=1),
    )

fig.suptitle(f'{deploy_name}: Pressure and Salinity Compensation', fontsize=16)
for axi in [ax00, ax01, ax02]:
    axi.set_facecolor('0.9')
Figure
No description has been provided for this image
In [70]:
##Quantify the differnce: percent difference between PS compensated data and raw [O2] in umol/L 
difference = ((ds.ox_P_S_corr -ds2.oxygen_concentration)/ds2.oxygen_concentration)*100

print(f'min O2: {np.nanmin(ds.ox_P_S_corr.values)}')
print(f'max O2: {np.nanmax(ds.ox_P_S_corr.values)}')
print(f'mean O2: {np.nanmean(ds.ox_P_S_corr.values)}')
print(f'std O2: {np.nanstd(ds.ox_P_S_corr.values)}')

min O2: 4.209563794508424
max O2: 365.7460348397695
mean O2: 136.89931075996816
std O2: 88.72771874542293

There is a lot of up/down asymmetry that is in part associated with a sensor time lag

Add quality flags - will need to avoid making spikes larger in the next step¶

We flag good data as QC 1 and bad data as QC 4, following Argo quality-control notation (see https://argo.ucsd.edu/data/how-to-use-argo-files/). Oxygen optodes are not pumped like the Glider Pumped CTD (GPCTD) used on some gliders, where biological material can build up in the tubing, especially in highly productive regions (see the C-PROOF GPCTD processing report). Although less common, material can stick to an optode and produce bad oxygen data.

A more robust flagging method could be developed to identify small-scale spikes that are not removed by the time-response correction or by gridding, where data are averaged into one-metre depth bins.

¶

An oxygen_QC variable was created and initialized with a quality-control flag of 1 (good data) for all observations. Oxygen anomalies were then visually inspected in isopycnal space. Profiles containing anomalous spikes that extended through a substantial portion of the water column were manually assigned a quality-control flag of 4 (bad data).

In [71]:
 ### Create oxygen_QC variable. Set all values to QC1

ds2["oxygen_concentration_QC"] = xr.ones_like(
    ds2.oxygen_concentration,
    dtype=np.int8,
)
In [72]:
outfile = f'{filepath}/{deploy_name}_O2_2.nc'
ds2.to_netcdf(outfile)
ncprocess.make_gridfiles(
        outfile,
        f'{filepath}',
        deployfile,
        fnamesuffix='_O2_2',
        maskfunction=utils.maskQC4,
        max_gap=150
    )
Out[72]:
'deployments/dfo-bb046/dfo-bb046-20220707//dfo-bb046-20220707_grid_O2_2.nc'
In [73]:
ds = xr.open_dataset(f'{filepath}/{deploy_name}_grid_O2_2.nc')

#####put in isopycnal space 

pot_density = gsw.pot_rho_t_exact(ds.salinity, ds.potential_temperature, ds.pressure,0)
# Compute the mission-mean potential density at each depth.
# This mean density profile defines the isopycnal coordinate onto which
# oxygen concentrations are interpolated.rho_avg = (pot_density.mean(dim='time'))
rho_avg = (pot_density.mean(dim='time'))
iso_density= xr.DataArray(rho_avg.values, dims=['depth']).expand_dims(dim={"time": len(ds.time)}, axis=1) 

### Loop through every profile to map onto the isopycnal coordinate
interpolated_oxy = []

for i in range(len(ds.time)):
    x = iso_density.isel(time=i).values
    xp = pot_density.isel(time=i).values
    yp = ds.ox_P_S_corr.isel(time=i).values

    # Keep only finite values
    mask = np.isfinite(xp) & np.isfinite(yp)

    if mask.sum() < 2:
        interpolated_oxy.append(np.full_like(x, np.nan, dtype=float))
    else:
        interpolated_oxy.append(
            np.interp(
                x,
                xp[mask],
                yp[mask],
                left=np.nan,
                right=np.nan,
            )
        )

# Subtract the mean oxygen concentration on each isopycnal and divide by
# the corresponding standard deviation to obtain standardized anomalies.
ds["iso_oxy"] = (("depth", "time"), np.transpose(interpolated_oxy))
mn_oxy = ds["iso_oxy"].mean("time", skipna=True)
std_oxy = ds["iso_oxy"].std("time", skipna=True)

ds['oxy_anomalies'] = (ds.iso_oxy  -mn_oxy )/std_oxy
In [74]:
NUM_PROFILES = np.nanmax(ds.profile_index)

xlim_1 = [0, int(NUM_PROFILES/4)]
xlim_2 = [int(NUM_PROFILES/4), int(NUM_PROFILES/4*2)]
xlim_3 = [int(NUM_PROFILES/4*2), int(NUM_PROFILES/4*3)]
xlim_4 = [int(NUM_PROFILES/4*3), NUM_PROFILES]



Y_LIMS = [400, 0]  
cmap='RdBu_r'
vmin=-1
vmax=1


fig, axs = plt.subplots(4,  #height_ratios=[1, 4], 
                        figsize = [10,9],
                        layout='constrained', sharex=False,sharey=True)

profile_lims = xlim_1

ax = axs[0]
ds_sub = ds.where((ds.profile_index >=profile_lims[0]) & (ds.profile_index <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, -ds_sub.depth, ds_sub['oxy_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'oxygen anomaly from $[O2_{raw}]$')


########
profile_lims = xlim_2

ax = axs[1]
ds_sub = ds.where((ds.profile_index >=profile_lims[0]) & (ds.profile_index <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, -ds_sub.depth, ds_sub['oxy_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'oxygen anomaly from $[O2_{raw}]$')

######
profile_lims = xlim_3

ax = axs[2]
ds_sub = ds.where((ds.profile_index >=profile_lims[0]) & (ds.profile_index <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, -ds_sub.depth, ds_sub['oxy_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'oxygen anomaly from $[O2_{raw}]$')

#########
profile_lims = xlim_4

ax = axs[3]
ds_sub = ds.where((ds.profile_index >=profile_lims[0]) & (ds.profile_index <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, -ds_sub.depth, ds_sub['oxy_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'oxygen anomaly from $[O2_{raw}]$')
 
print('Zooming in along the glider deployment to visualize oxygen spikes.')
 

Zooming in along the glider deployment to visualize oxygen spikes.
Figure
No description has been provided for this image

There aren't any spikes in the data. We wont flag any profiles as QC4.

It does look very stripy in the deep section. This is known to be a very productive area, so there may have been something stuck to the optode. We will not optiimize the time response in this area just in case

2.3 Convert O₂ from μmol L⁻¹ to μmol kg⁻¹¶

μmol kg⁻¹ is the more common unit. We will use these units to apply the following corrections.

In [75]:
### convert glider oxygen to umol/kg
ts3 = xr.open_dataset(f'{filepath}/{deploy_name}_O2_2.nc')
ds3 = xr.open_dataset(f'{filepath}/{deploy_name}_grid_O2_2.nc')

ds3['oxygen_umol_kg'] = ds3['ox_P_S_corr']/(ds3['potential_density']/1000) #(umol/L to umol/kg)
ds3['oxygen_raw_umol_kg'] = ds3['oxygen_concentration']/(ds3['potential_density']/1000) #(umol/L to umol/kg)

ts3['oxygen_umol_kg'] = ts3['ox_P_S_corr']/(ts3['potential_density']/1000) #(umol/L to umol/kg)
ts3['oxygen_raw_umol_kg'] = ts3['oxygen_concentration']/(ts3['potential_density']/1000) #(umol/L to umol/kg)

ts3.to_netcdf(f'{filepath}/{deploy_name}_O2_3.nc')
ds3.to_netcdf(f'{filepath}/{deploy_name}_grid_O2_3.nc')

2.4 Determine the sensor response time lag correction¶

The time response correction accounts for the finite response time of the optode relative to the glider's vertical velocity. Following Bittig et al. (2014), the sensor response can be approximated as a first-order recursive low-pass filter with time constant $\tau$. The inverse of this filter, referred to here as the sharpening filter, is applied to recover the true oxygen signal.

$y_i = a\,y_{i-1} + b\,x_i , \qquad $

with coefficients $ a = \exp\!\left(-\frac{\Delta t_i}{\tau}\right), \qquad b = 1-a ,$ where $\Delta t_i = t_i - t_{i-1}$.

A mission-specific response time, $\tau$, is used when applying the sharpening filter to the pressure- and salinity-compensated oxygen concentration. The optimal value of $\tau$ varies between deployments because it depends primarily on temperature and the thickness of the stagnant boundary layer surrounding the sensor. Bittig et al. (2014) reported that $\tau$ is typically expected to range from 25–80 s for Aanderaa optodes deployed on gliders. In contrast, the RINKO AROD-FT optode has a substantially faster response time (approximately 0–4 s), and therefore a time response correction is generally unnecessary for these sensors.

To limit the amplification of isolated spikes, the pressure- and salinity-compensated oxygen record is first smoothed with a 1.5 s (three-sample) running mean before applying the sharpening filter. After the sharpening filter has been applied to oxygen measurements collected by Aanderaa sensors, the corrected oxygen record is smoothed using a running mean with a window length of approximately $\tau/2$ s. This second smoothing step suppresses the high-frequency variability introduced by the inverse filter while preserving the lower-frequency oxygen variability.

Grid data onto isopycnals and fit response time using an objective error function¶

A mission-specific response time, $\tau$, is estimated by minimizing the following objective function. For each candidate value of $\tau$, the sharpening filter is applied to the pressure- and salinity-compensated oxygen concentration, and the corrected oxygen profiles are interpolated onto a mission-mean isopycnal grid. The objective function is then evaluated over a subset of 40 consecutive profiles between the average oxycline depth (approximately 70–150 m) and 1000 m:

\begin{equation} \frac{(\Delta O_2)^2} {\sigma\!\left(O_{2,\mathrm{iso}}\right)} \frac{N_0}{N} \end{equation}

where $\Delta O_2$ is the difference in oxygen concentration between consecutive profiles on each isopycnal, $\sigma(O_{2,\mathrm{iso}})$ is the standard deviation of oxygen concentration on that isopycnal, $N_0$ is the number of finite oxygen values contributing to the calculation, and $N$ is the total number of oxygen values considered. The factor $N_0/N$ reduces the influence of sparsely sampled isopycnals.

The oxycline is identified as the depth of the maximum vertical oxygen gradient in the pressure- and salinity-compensated oxygen concentration, $O_{2,\mathrm{PS}}$. Only data below the average oxycline depth are included in the optimization because upper-ocean variability is dominated by natural spatial and temporal variability rather than sensor response.

¶

There is only 1 deep section, so we will optimize tau using those data. This section is likely bad so we may not be able to constrain the response well.

In [76]:
fig,axs= plt.subplots(1,1, figsize=(10,5),sharex=True, sharey=True)
from tqdm import tqdm

### Look at 1 subsets of data 
ranges =[450,490]

ts3 = xr.open_dataset(f'{filepath}/{deploy_name}_O2_3.nc')

### Make a 40 profile subset 
ts_sub = ts3.where((ts3.profile_index > ranges[0]) & (ts3.profile_index < ranges[1]) & (ts3.depth > 150),drop=True)
    
### Loop through tau from 0-4s following Bittig 2014 
taus = np.arange(0, 6, 1)
errors = np.zeros(len(taus))

# apply to subset with progress bar
for nx, tau in enumerate(tqdm(taus, desc="Scanning tau", unit="tau")):
    ts3 = sns.apply_tau_correction(ts_sub, tau) # apply tau correction
    outfile = sns.make_gridfiles_adj(ts3, deployfile)  # grid the data
    oxy_iso, err2, weights,err = sns.get_error_iso(outfile,outfile.oxygen_concentration_corrected) # isopycnal error in umol/kg 
    weighted_error= err2*weights.values[:, np.newaxis] #weight the error^2
    errors[nx] = np.nanmean(weighted_error)
        
axs.set_title(f"{ranges} profile subset",fontsize=18)
axs.grid(True, alpha=0.3)
axs.plot(taus, errors, marker="o")
axs.grid(True, alpha=0.3)
axs.tick_params(axis='both', labelsize=18)


## Plot 
axs.set_xlabel("tau [s]", fontsize=18)
axs.set_ylabel(
    r"$\frac{(\Delta [O_{2,\mathrm{PS}}])^2}{\sigma([O_{2,\mathrm{PS}}])}"
    r"\,\frac{N_0}{N}$",
    fontsize=22)
fig.suptitle(f"tau scan: dfo-{deploy_name} ",fontsize=18 )
plt.show()

Scanning tau: 100%|██████████| 6/6 [00:02<00:00,  2.87tau/s]
Figure
No description has been provided for this image

The tau with the lowest objective error function is 4 s, the upper limit of the physically realistic range (Bittig et al, 2014).

Investigate how the asymmetry between upcasts and downcasts changes after applying the time response correction using the optimized value of $\tau$ and values of $\tau \pm 1$ s. Data are plotted as anomalies relative to the raw oxygen concentration (µmol kg$^{-1}$).

In [77]:
def tau_correction_plot(ds, title, ax):
    ts = ds.where(
        (ds.profile_index >= 450) &
        (ds.profile_index <= 490),
        drop=True
    )

    outfile3 = sns.make_gridfiles_adj(ts, deployfile)
    oxy_lag = outfile3.oxygen_concentration_corrected

    ss0, err02, weights0 ,err0= sns.get_error_iso(
        outfile3,
        outfile3.oxygen_raw_umol_kg # mumol/L --> mumol/kg
    )
    ss, err2, weights ,err= sns.get_error_iso(outfile3, oxy_lag)

    # Raw isopycnal oxygen anomaly
    sp0 = np.nanmean(ss0, axis=1)

    pc0 = ax[0].pcolormesh(
        ss0.profile,
        ss0.depth,
        ss0 - sp0[:, np.newaxis],
        cmap="RdBu_r",
        vmin=-20,
        vmax=20,
        rasterized=True
    )

    ax[0].set_ylabel("Isopycnal depth [m]", fontsize=16)
    ax[0].set_xlabel("Profile", fontsize=16)
    ax[0].set_title(r"Raw $[O_{2,\mathrm{iso}}]$", fontsize=16)

    # Corrected isopycnal oxygen anomaly
    sp = np.nanmean(ss, axis=1)

    pc1 = ax[1].pcolormesh(
        ss.profile,
        ss.depth,
        ss - sp[:, np.newaxis],
        cmap="RdBu_r",
        vmin=-20,
        vmax=20,
        rasterized=True
    )

    ax[1].set_title(title, fontsize=16)
    ax[1].set_xlabel("Profile", fontsize=16)

    # Weighted error
    weighted_error02 = err02 #* weights0.values[:, np.newaxis]
    weighted_error2 = err2 #* weights.values[:, np.newaxis]

    ax[2].plot(
        np.nanmean(weighted_error2, axis=1),#[weights != 0],
        outfile3.depth,#[weights != 0],
        label="Corr"
    )
    ax[2].plot(
        np.nanmean(weighted_error02, axis=1),#[weights0 != 0],
        outfile3.depth,#[weights0 != 0],
        label="Orig"
    )

    ax[2].grid(axis="x")
    ax[2].set_ylim([1000, 0])
    ax[2].legend(fontsize=14)
    ax[2].set_xlabel("Mean weighted error", fontsize=16)
    ax[2].set_title(
        r"$\frac{(\Delta [O_{2,\mathrm{PS}}])^2}{\sigma([O_{2,\mathrm{PS}}])}$",
        #r"\frac{N_0}{N}$",
        fontsize=16
    )
    ax[2].axvline(x=0, c="k", ls="--")

    # Unweighted error
    ax[3].plot(
        np.nanmean(err, axis=1),
        outfile3.depth,
        label="Corr"
    )
    ax[3].plot(
        np.nanmean(err0, axis=1),
        outfile3.depth,
        label="Orig"
    )

    ax[3].grid(axis="x")
    ax[3].set_ylim([1000, 0])
    ax[3].legend(fontsize=14, loc="lower left")
    ax[3].set_xlabel("Mean error", fontsize=16)
    ax[3].set_title(
        r"$\frac{(\Delta [O_{2,\mathrm{PS}}])}{\sigma([O_{2,\mathrm{PS}}])}$",
        fontsize=16
    )
    ax[3].axvline(x=0, c="k", ls="--")

    return pc1


taus = (2, 3, 4)

fig, axs = plt.subplots(
    len(taus),
    4,
    figsize=(14, 9),
    sharey=True,
    layout="constrained"
)

ts3 = xr.open_dataset(f"{filepath}/{deploy_name}_O2_3.nc")

for i, tau in enumerate(taus):
    ds_test = sns.apply_tau_correction(ts3, tau)

    pc = tau_correction_plot(
        ds_test,
        title=rf"$\tau$ = {tau} s",
        ax=axs[i, :]
    )

cb = fig.colorbar(
    pc,
    ax=axs[:, :2],
    shrink=0.6
)

cb.ax.tick_params(labelsize=16)
cb.set_label(
    r"$[O_{2,\mathrm{PS}}] - \overline{[O_{2,\mathrm{PS}}]}$",
    fontsize=16
)

for ax in axs.flat:
    ax.tick_params(axis="both", labelsize=14)
    ax.title.set_fontsize(14)
    ax.xaxis.label.set_size(14)
    ax.yaxis.label.set_size(14)

labels = [chr(97 + i) for i in range(12)]

for ax, label in zip(axs.flat, labels):
    ax.text(
        0.02,
        0.98,
        label,
        transform=ax.transAxes,
        fontsize=18,
        fontweight="bold",
        va="top",
        ha="left",
        bbox=dict(
            facecolor="white",
            edgecolor="none",
            alpha=0.6,
            pad=0.2
        )
    )
Figure
No description has been provided for this image

We are going to apply tau=4 since that optimized the objective error function for a good section of data (including part of the productive area as well).

We are not going to smooth the data since a short time response was corrected for.

In [78]:
#####We are going to ensure we are looking at the correct upcast/downcast pairs 

fig,axs=plt.subplots(5)

ranges = np.arange(451,461)
for i, val in enumerate(ranges): 
    
    ts_sub = ds_test.where(ds_test.profile_index == val,drop=True
    )
    
    if ts_sub.profile_direction.mean() == -1:
        c = 'red'
        label='up: '
    else: 
        c = 'blue'
        label='down: '
    
    ax_idx = i // 2  
    
    axs[ax_idx].plot(ts_sub.time, ts_sub.depth,c=c, ls='dashdot',label =f' {label} PS corrected' if i<2  else None  )
    axs[ax_idx].set_ylim(1000,0)

axs[0].legend(loc= 'lower left')
axs[4].set_xlabel('Time (s)')
axs[0].set_ylabel('Depth (m)')
 

fig.suptitle('Find up/downcast pairs')
Out[78]:
Text(0.5, 0.98, 'Find up/downcast pairs')
Figure
No description has been provided for this image
In [79]:
tau =4 
ds_test = sns.apply_tau_correction(ts3, tau)


fig, axs = plt.subplots(1, 5, figsize=(11,8), sharex=True, sharey=True)
ranges = np.arange(880, 890)
for i, val in enumerate(ranges): 
    
    ts_sub = ds_test.where(ds_test.profile_index == val).where(
        (ds_test.depth <= 200) & (ds_test.depth >= 50)
    )
    
    if ts_sub.profile_direction.mean() == -1:
        c = 'red'
        label='up: '
        cc = 'orange'
    else: 
        c = 'blue'
        label='down: '
        cc = 'lightgreen'
    
    ax_idx = i // 2  
    
    axs[ax_idx].plot(ts_sub.ox_P_S_corr_smoothed, -ts_sub.depth, c=c, ls='dashdot',label =f' {label} PS correct' if i<2  else None  )
    axs[ax_idx].plot(ts_sub.oxygen_concentration_corrected, -ts_sub.depth, c=c, label =f' {label} tau={tau}' if i<2 else None )
    #axs[ax_idx].plot(ts_sub.oxygen_concentration_corrected_smoothed_tau_2_box,- ts_sub.depth, c=cc, label =f' smoothed tau/2' if i<2 else None )


# formatting
prof = []

for j in range(5):
    axs[j].grid(True)
    axs[j].set_title(f'Profiles {str(ranges[2*j:2*j+2])}' )

axs[0].legend(loc= 'lower left')
axs[1].set_xlabel('$[O\_{2,PS}]$ ($\mu$mol/kg)')
axs[0].set_ylabel('Depth')

fig.suptitle(f'dfo-{deploy_name}; tau = {tau}')

labels = ['a', 'b', 'c', 'd','e']

for axi, label in zip(axs.flat, labels):
    axi.text(
        0.02, 0.98,
        label,
        transform=axi.transAxes,
        fontsize=14,
        fontweight='bold',
        va='top',
        ha='left'
    )
Figure
No description has been provided for this image

Save data and grid

In [80]:
outfile = f'{filepath}/{deploy_name}_O2_4.nc'
ds_test.to_netcdf(outfile)
ncprocess.make_gridfiles(
        outfile,
        f'{filepath}',
        deployfile,
        fnamesuffix='_O2_4',
        maskfunction=utils.maskQC4,
        max_gap=150
    )
Out[80]:
'deployments/dfo-bb046/dfo-bb046-20220707//dfo-bb046-20220707_grid_O2_4.nc'

Plot oxygen anomalies in depth space to illustrate the observed upcast–downcast asymmetry.

In [81]:
def depth_space(ds, min_ind, max_ind, tau):
    sub = ds.where(
        (ds.profile_index >= min_ind) &
        (ds.profile_index <= max_ind),
        drop=True
    )

    vmin, vmax = -50, 50
    cmap = "RdBu_r"

    fig, ax = plt.subplots(
        1, 2,
        layout="constrained",
        figsize=(6, 3),
        sharey=True,
        sharex=True
    )

    mean = ds.oxygen_raw_umol_kg.mean("time")

    x = sub.profile_index[5,:]

    im0 = ax[0].pcolormesh(
        x,
        sub.depth,
        sub.oxygen_umol_kg - mean,
        rasterized=True,
        vmin=vmin,
        vmax=vmax,
        cmap=cmap
    )
    ax[0].set_title(r"$[O_{PS}]$")
    ax[0].set_ylabel("Depth (m)")
    ax[0].set_ylim(1000, 0)

    im1 = ax[1].pcolormesh(
        x,
        sub.depth,
        sub.oxygen_concentration_corrected - mean,
        rasterized=True,
        vmin=vmin,
        vmax=vmax,
        cmap=cmap
    )
    ax[1].set_title(rf"[$O_2,\tau], \tau = {tau}$ s")

    cbar = fig.colorbar(
        im1,
        ax=ax,
        shrink=0.9,
        label=r"$[O_2] - \overline{[O_2]}$ ($\mu$mol kg$^{-1}$)"
    )

    labels = ["a", "b", "c"]

    for axi, label in zip(ax.flat, labels):
        axi.text(
            0.02, 0.98,
            label,
            transform=axi.transAxes,
            fontsize=14,
            fontweight="bold",
            va="top",
            ha="left",
            bbox=dict(facecolor="white", edgecolor="none", alpha=0.6, pad=0.2)
        )
        axi.set_xlabel("Profile index")

    return fig
In [82]:
ds_test_grid = xr.open_dataset(f'{filepath}/{deploy_name}_grid_O2_4.nc')
fig = depth_space(ds_test_grid,450,510,tau )
Figure
No description has been provided for this image

2.6 Determine the sensor drift correction¶

The sensor drift correction accounts for changes in the optode response that occur as the sensor ages. Although sensor drift generally increases with the total number of measurements collected, optodes have been shown to drift more rapidly during storage than during deployment. The same general drift-correction procedure is applied to both Aanderaa and RINKO AROD-FT optodes.

The drift correction is applied after the pressure, salinity, and time response corrections. A linear calibration is used to relate the corrected glider oxygen concentration to a trusted comparison profile,

[ O_{2,\mathrm{adj}} = m,O_{2} + b, ]

where (O_{2}) is the pressure-, salinity-, and time-response corrected oxygen concentration, (m) is the gain, and (b) is the offset.

The correction coefficients are determined by comparing the glider oxygen observations with a trusted reference profile collected nearby in space and time. Ideally, the comparison profile consists of Winkler-calibrated CTD oxygen measurements collected during the deployment. When a suitable shipboard profile is unavailable, comparison profiles are selected using a hierarchy of reference datasets. First, a nearby CTD profile collected with the same optode is used when available. If no suitable CTD profile exists, a previously calibrated deployment from the same optode is considered. Finally, for missions extending sufficiently offshore where deep oxygen is relatively stable, the World Ocean Atlas (WOA) climatology may be used. Comparisons are rejected when no reference profile is considered sufficiently representative.

To reduce the effects of vertical heave and internal waves, comparisons are performed in potential density space rather than depth space. Both the glider and comparison profile are interpolated onto common potential density bins before fitting. Data within the pycnocline are excluded from the comparison because strong vertical oxygen gradients make the comparison sensitive to small density mismatches. Only near-surface waters and observations below the pycnocline are retained.

A linear least-squares regression is then performed between the binned glider and comparison oxygen concentrations to determine the gain and offset. The resulting correction coefficients are evaluated by examining the agreement between the corrected glider profile and the comparison profile throughout the retained portions of the water column. The gain is expected to remain close to unity, while the offset is typically small. If the fitted correction is physically unreasonable or does not improve agreement with the comparison profile, the comparison is re-evaluated or rejected.

Although the correction coefficients are estimated using oxygen observations interpolated into potential density space, the resulting linear correction is applied to the original, unbinned oxygen measurements.

In [83]:
### open trustworthy comparison CTD 
ctd = xr.open_dataset(f'{filepath}/auxdata/O2_CTD_comparison.nc')

### open WOA P16 data to compare (since glider travels sufficiently far offshore, approximately west of 134W.)
woaP16 = xr.open_dataset(f'~/processing/deployments/otherdata/WOA_climatology49_5N_134_5W.nc',decode_times=False)

ts4= xr.open_dataset(f'{filepath}/{deploy_name}_O2_4.nc')
ds4 = xr.open_dataset(f'{filepath}/{deploy_name}_grid_O2_4.nc')

Find closest glider profiles to the comparison casts

In [84]:
def to_datetime64ns(x):
    """
    Convert scalar or array-like to numpy datetime64[ns].
    Returns np.datetime64('NaT') for unparseable values.
    """
    try:
        if np.isscalar(x):
            return pd.to_datetime(x, errors="coerce").to_datetime64()

        return pd.to_datetime(x, errors="coerce").to_numpy(dtype="datetime64[ns]")

    except Exception:
        return np.array([], dtype="datetime64[ns]")


def haversine_km(lon1, lat1, lon2, lat2):
    """
    Great-circle distance in km.
    """
    R = 6371.0

    lon1 = np.deg2rad(lon1)
    lat1 = np.deg2rad(lat1)
    lon2 = np.deg2rad(lon2)
    lat2 = np.deg2rad(lat2)

    dlon = lon2 - lon1
    dlat = lat2 - lat1

    a = (
        np.sin(dlat / 2) ** 2
        + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
    )

    return 2 * R * np.arcsin(np.sqrt(a))


def _get_lat_lon(ds):
    """
    Get latitude/longitude from dataset using either lat/lon or
    latitude/longitude names.
    """
    if "lat" in ds:
        lat = ds["lat"]
    elif "latitude" in ds:
        lat = ds["latitude"]
    else:
        raise ValueError("Dataset must contain 'lat' or 'latitude'.")

    if "lon" in ds:
        lon = ds["lon"]
    elif "longitude" in ds:
        lon = ds["longitude"]
    else:
        raise ValueError("Dataset must contain 'lon' or 'longitude'.")

    return lat, lon


def _representative_time(time_da):
    """
    Return one representative datetime64[ns] from a time DataArray.
    """
    tvals = to_datetime64ns(time_da.values)

    if np.isscalar(tvals):
        if np.isnat(tvals):
            return np.datetime64("NaT")
        return np.datetime64(tvals, "ns")

    tvals = tvals[~np.isnat(tvals)]

    if len(tvals) == 0:
        return np.datetime64("NaT")

    return np.sort(tvals)[len(tvals) // 2]


def closest_glider_comparison_profile(
    glider_ds,
    comparison_ds,
    use_time=True,
    d0_km=20,
    t0_days=3,
):
    """
    Find closest glider profile(s) to comparison profile location(s).

    If comparison_ds has valid time, latitude, and longitude, the closest
    profile is selected using both space and time:

        score = (distance / d0_km)^2 + (time_difference / t0_days)^2

    If time is unavailable, invalid, or use_time=False, only distance is used.

    Parameters
    ----------
    glider_ds : xarray.Dataset
        Glider dataset containing ``latitude``, ``longitude``, ``profile_index``,
        and ``time``.

    comparison_ds : xarray.Dataset
        Comparison dataset containing either ``lat``/``lon`` or
        ``latitude``/``longitude``. May also contain ``time``.

    use_time : bool, optional
        If True, use both distance and time when valid comparison time is
        available. If False, use distance only.

    d0_km : float, optional
        Distance scale used in the combined space-time score.

    t0_days : float, optional
        Time scale used in the combined space-time score.

    Returns
    -------
    closest_profile : int or numpy.ndarray
        Closest glider profile index. Returns one value for one comparison
        profile/location, or one value per CTD station/cast.
    """

    glider_ds = glider_ds.set_coords("profile_index")

    # One lat/lon/time per glider profile
    lat_prof = glider_ds.latitude.groupby(glider_ds.profile_index).mean()
    lon_prof = glider_ds.longitude.groupby(glider_ds.profile_index).mean()
    time_prof = glider_ds.time.groupby(glider_ds.profile_index).min()

    comp_lat, comp_lon = _get_lat_lon(comparison_ds)

    # Keep these dimensions if they exist, so CTD stations/casts stay separate
    keep_dims = ["station_ind", "cast", "station", "mission_ind"]

    # Collapse point dimensions such as depth or row
    point_dims = [d for d in comp_lat.dims if d not in keep_dims]

    if len(point_dims) > 0:
        comp_lat = comp_lat.mean(dim=point_dims, skipna=True)
        comp_lon = comp_lon.mean(dim=point_dims, skipna=True)

    # Space-only part of score
    dist_km = haversine_km(
        lon_prof,
        lat_prof,
        comp_lon,
        comp_lat,
    )

    score = (dist_km / d0_km) ** 2

    # Add time part of score only when valid comparison time exists
    if use_time and "time" in comparison_ds:
        comp_time_da = comparison_ds["time"]

        point_dims_time = [d for d in comp_time_da.dims if d not in keep_dims]

        if len(point_dims_time) > 0:
            comp_time_da = comp_time_da.min(dim=point_dims_time, skipna=True)

        # If comparison has multiple stations/casts, loop over them
        if len(comp_time_da.dims) > 0:
            dt_score = xr.full_like(score, np.nan, dtype=float)

            for index in np.ndindex(comp_time_da.shape):
                selector = {
                    dim: comp_time_da[dim].values[i]
                    for dim, i in zip(comp_time_da.dims, index)
                }

                ctime = _representative_time(comp_time_da.sel(selector))

                if np.isnat(ctime):
                    continue

                dt_days = np.abs(
                    (time_prof - ctime) / np.timedelta64(1, "D")
                )

                dt_score.loc[selector] = (dt_days / t0_days) ** 2

            if np.isfinite(dt_score).any():
                score = score + dt_score

        else:
            ctime = _representative_time(comp_time_da)

            if not np.isnat(ctime):
                dt_days = np.abs(
                    (time_prof - ctime) / np.timedelta64(1, "D")
                )

                score = score + (dt_days / t0_days) ** 2

    closest_profile = score.argmin("profile_index")

    # Return plain numbers, not xarray objects
    if closest_profile.size == 1:
        return int(closest_profile.values)

    return closest_profile.values.astype(int)
In [85]:
woa = woaP16.interp(depth=np.arange(0, 1000))
closest_profile_woa = closest_glider_comparison_profile(ts4, woa) 
print("Closest profile_index to WOA P16:", closest_profile_woa)

ctd = ctd.interp(depth=np.arange(0, 1000))
closest_profile_ctd = closest_glider_comparison_profile(ts4, ctd) 
print("Closest profile_index to comparison CTD:", closest_profile_ctd)

Closest profile_index to WOA P16: 370
Closest profile_index to comparison CTD: 896
In [86]:
ind0 = closest_profile_ctd
inds0 =np.arange(ind0-5,ind0)

ind_woa = closest_profile_woa
inds_woa = np.arange(ind_woa-2,ind_woa+3)

###
best_ctd1 = ctd

best_ctd1['time'] =np.mean(ctd['time'] )

Plot comparisons¶

In [87]:
fig = plt.figure(
    figsize=(10, 8),
    constrained_layout=True
)

gs = fig.add_gridspec(
    2,
    2,
    height_ratios=[1, 1]
)

# Top row (full width map)
ax0 = fig.add_subplot(gs[0, :])

# Bottom row
ax1 = fig.add_subplot(gs[1, 0])
ax2 = fig.add_subplot(gs[1, 1])

# So the rest of your code still works
axs = [ax0, ax1, ax2]

# ==========================================================
# PANEL 1 — MAP
# ==========================================================

axs[0].plot(
    ds4.longitude,
    ds4.latitude,
    '.',
    color='lightgray',
    label='All glider profiles'
)

axs[0].plot(
    ds4.longitude[inds0],
    ds4.latitude[inds0],
    'ok',
    label='Glider profiles to compare to CTD1'
)

axs[0].plot(
    woa.lon,
    woa.lat,
    color='orange',
    marker='*',
    markersize=12,
    linestyle='None',
    label='WOA decadal climatology'
)


diff1 = np.round(
    np.abs(
        (best_ctd1.time.values - ds4.time.values[inds0])
        / np.timedelta64(1, "D")
    ).mean(),
    2
)

axs[0].plot(
    best_ctd1.longitude,
    best_ctd1.latitude,
    color='pink',
    marker='*',
    markersize=12,
    linestyle='None',
    label=f'CTD1: {diff1} day diff'
)



axs[0].set_ylabel('Latitude [$^o$N]')
axs[0].set_xlabel('Longitude [$^o$W]')
axs[0].set_title('Profile locations')
axs[0].grid()
axs[0].legend(loc='lower left')
#axs[0].legend(loc='center left', bbox_to_anchor=(0.2, -2))

mean_lat = float(ts4.latitude.mean())

axs[0].set_aspect(
    1 / np.cos(np.deg2rad(mean_lat))
)

# ==========================================================
# PANEL 2 — DEPTH
# ==========================================================

axs[1].plot(
    ds4.oxygen_concentration_corrected[:, inds0],
    ds4.depth,
    'k',
    linewidth=0.5
)
axs[1].plot([],[],'k', linewidth=2, label = '5 closest glider profiles to CTD1')


axs[1].plot(
    best_ctd1.oxygen_umol_kg,
    best_ctd1.depth,
    color='pink',
    linewidth=3,
    label='CTD1'
)




axs[1].set_ylim(1000, 0)
axs[1].set_ylabel('Depth [m]')
axs[1].set_xlabel('Oxygen [$\\mu$mol/kg]')
axs[1].set_title('Oxygen vs depth')
axs[1].grid()
axs[1].legend(fontsize=10, loc='lower right')

# ==========================================================
# PANEL 3 — ISOPYCNAL
# ==========================================================

axs[2].plot(
    ds4.oxygen_concentration_corrected[:, inds0],
    ds4.potential_density_adjusted[:, inds0],
    'k',
    linewidth=0.5,
    label='Glider profiles to compare to CTD1'
)


axs[2].plot(
    best_ctd1.oxygen_umol_kg,
    best_ctd1.potential_density,
    color='pink',
    linewidth=3,
    label='CTD1'
)



axs[2].invert_yaxis()

axs[2].set_ylabel('Potential density [kg/m$^3$]')
axs[2].set_xlabel('Oxygen [$\\mu$mol/kg]')
axs[2].set_title('Oxygen vs potential density')
axs[2].grid()


fig.suptitle(
    'Glider oxygen comparison profiles',
    fontsize=16
)


labels = ['a', 'b', 'c', 'd']

for axi, label in zip(axs, labels):
    axi.text(
        0.02,
        0.98,
        label,
        transform=axi.transAxes,
        fontsize=18,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(
            facecolor='white',
            edgecolor='none',
            alpha=0.6,
            pad=0.2
        )
    )
Figure
No description has been provided for this image

The mission does not go far enough offshore to compare to the WOA

Select glider profiles and assess the depth range to compare¶

Glider profiles selected for the drift correction should first be evaluated to ensure they sampled the same water mass as the comparison CTD. This assessment is performed by comparing temperature and temperature anomalies in density space. Profiles exhibiting substantially different temperature structure or anomalies are excluded from the drift estimation because they likely sampled different water masses.

In [91]:
def glider_comparison(ax1, ax2,ax3, glider_indices, comparison, comparison_title):
    bins =  np.arange(1024.5, 1028.5, 0.1)
    comparison = comparison.interp(depth=ds4.depth)
    if comparison_title == 'WOA':
        oxy = comparison.mean_oxygen
        temp = comparison.mean_temperature
    else: 
        oxy = comparison.oxygen_umol_kg
        temp = comparison.temperature

    mask = (
        np.isfinite(oxy) &
        np.isfinite(comparison.potential_density)
    )

    comparison_iso = np.interp(
        bins,
        comparison.potential_density[mask],
        oxy[mask],
        left=np.nan,
        right=np.nan
    )
    comparison_iso_temp = np.interp(
        bins,
        comparison.potential_density[mask],
        temp[mask],
        left=np.nan,
        right=np.nan
    )

    gr1_oxy = np.full((len(bins), len(glider_indices)), np.nan)
    gr1_temp = np.full((len(bins), len(glider_indices)), np.nan)

    c= ['red','green','orange','blue','purple','magenta']
    for i in range(len(glider_indices)):

        rho = ds4.potential_density_adjusted.values[:, glider_indices[i]]
        oxy = ds4.oxygen_concentration_corrected[:, glider_indices[i]]
        temp = ds4.temperature.values[:, glider_indices[i]]
        mask = np.isfinite(rho) & np.isfinite(oxy)

        if np.sum(mask) < 2:
            continue

        rho_valid = rho[mask]
        oxy_valid = oxy[mask]
        temp_valid = temp[mask]

        order = np.argsort(rho_valid)

        gr1_oxy[:, i] = np.interp(
            bins,
            rho_valid[order],
            oxy_valid[order],
            left=np.nan,
            right=np.nan
        )
        gr1_temp[:, i] = np.interp(
            bins,
            rho_valid[order],
            temp_valid[order],
            left=np.nan,
            right=np.nan
        )
    temp_std = np.nanstd(gr1_temp, axis=1)

    ########### 
    ax1.plot(comparison_iso_temp,bins, 'k',lw=4, label=comparison_title)
    for i in range(len(glider_indices)):
        ax1.plot(gr1_temp[:, i]  , bins, color=c[i], label=f'index: {glider_indices[i]}')
    ax1.set_ylabel('Potential density [kg/m$^3$]')
    ax1.set_xlabel('Temperature along isopycnals [°C]')
    ax1.legend(loc='lower left')
    ax1.set_ylim(1028,1025)
    ax1.grid()

    for i in range(len(glider_indices)):
        anom = (gr1_temp[:, i] -comparison_iso_temp)
        ax2.plot(anom , bins,c=c[i], label=f'index: {glider_indices[i]}', alpha=0.2)
        kernel = np.ones(3) / 3
        anom_smoothed = np.convolve(anom, kernel, mode='same')

        ax2.plot(anom_smoothed, bins,
             color=c[i])

    ax2.set_ylabel('Potential density [kg/m$^3$]')
    ax2.set_xlabel('Temperature anomaly from CTD along isopycnals [°C]')
    ax2.legend()
    ax2.set_ylim(1028,1025)
    ax2.set_xlim(-0.45,0.6)
    ax2.axvline(x=0,c='k',ls='--')
    ax2.axvline(x=0.2,c='gray',ls='--')
    ax2.axvline(x=-0.2,c='gray',ls='--')
    ax2.grid()


    ax3.plot(comparison_iso,bins, 'k',lw=4,  label=comparison_title)
    for i in range(len(glider_indices)):
        ax3.plot(gr1_oxy[:, i], bins, color=c[i], label=f'index: {glider_indices[i]}')
    ax3.set_ylabel('Potential density [kg/m$^3$]')
    ax3.set_xlabel('Oxygen [$\\mu$mol/kg]')
    ax3.legend(loc='lower right')
    ax3.set_ylim(1028,1025)
    ax3.grid()
    
    if comparison_title != 'WOA':
        diff = np.round(
            np.abs(
                (best_ctd1.time.values - ds4.time.values[glider_indices])
                / np.timedelta64(1, "D")
            ).mean(),
            2
        )
        ax1.set_title(f'{comparison_title}, mean {diff} day diff')
    else:
        ax1.set_title(f'{comparison_title}')
In [89]:
ds4 = ds4.where( (ds4.depth > 50))
In [94]:
fig,axs=plt.subplots(3,1, figsize=(6,10), constrained_layout=True)
glider_comparison(axs[0],axs[1], axs[2],inds0, best_ctd1, 'CTD1')

# 

labels = ['a', 'b', 'c', 'd', 'e', 'f','g','h','i']

for ax, label in zip(axs.flat, labels):
    ax.text(
        0.02,
        0.98,
        label,
        transform=ax.transAxes,
        fontsize=14,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(
            facecolor='white',
            edgecolor='none',
            alpha=0.6,
            pad=0.2
        )
    )
fig.suptitle(f'{deploy_name}', fontsize=16)
Out[94]:
Text(0.5, 0.98, 'dfo-bb046-20220707')
Figure
No description has been provided for this image

We will look at the profiles 891-893. The others are too different

Map [O2] onto potential density surfaces and estimate the gain and offset¶

Ensure that the surface-layer[O2] are in the same depth-bin

In [95]:
def glider_comparison_weighted(ax1,ax2, ax3, glider_indices, comparison, comparison_title, limit=None):
    bins =  np.arange(1025.5,1028.5,0.04)

    if comparison_title == 'WOA':
        oxy = comparison.mean_oxygen
        err = comparison.oxygen_se
    else:
        oxy = comparison.oxygen_umol_kg
        err = None

    comp_mask = np.isfinite(comparison.potential_density) & np.isfinite(oxy)

    comparison_iso = np.interp(
        bins,
        comparison.potential_density[comp_mask],
        oxy[comp_mask],
        left=np.nan,
        right=np.nan
    )
    comparison_dpth = np.interp(
        bins,
        comparison.potential_density[comp_mask],
        comparison.depth[comp_mask],
        left=np.nan,
        right=np.nan
    )

    if err is not None:
        err_iso = np.interp(
            bins,
            comparison.potential_density[comp_mask],
            err[comp_mask],
            left=np.nan,
            right=np.nan
        )
    else:
        err_iso = None

    gr1_oxy = np.full((len(bins), len(glider_indices)), np.nan)
    gr1_dpth = np.full((len(bins), len(glider_indices)), np.nan)


    for i, ind in enumerate(glider_indices):
        
        if limit is None:
            rho = ds4.potential_density_adjusted[:, ind].values
            oxy_glider = ds4.oxygen_concentration_corrected[:, ind]
            dpth = ds4.depth[:]
        else: 
            rho = ds4.potential_density_adjusted[:, ind].where(ds4.potential_density_adjusted[:, ind] > limit).values
            oxy_glider = ds4.oxygen_concentration_corrected[:, ind].where(ds4.potential_density_adjusted[:, ind] > limit)
            dpth = ds4.depth[:].where(ds4.potential_density_adjusted[:, ind] > limit)

    
        m = np.isfinite(rho) & np.isfinite(oxy_glider)

        rho_valid = rho[m]
        oxy_valid = oxy_glider[m]

        order = np.argsort(rho_valid)

        gr1_oxy[:, i] = np.interp(
            bins,
            rho_valid[order],
            oxy_valid[order],
            left=np.nan,
            right=np.nan
        )
        gr1_dpth[:, i] = np.interp(
            bins,
            rho_valid[order],
            dpth[order],
            left=np.nan,
            right=np.nan
        )

    if limit is not None:
        ax2.axvline(x=limit, c='r', ls='--', label= 'Min. density compared')
        ax1.axhline(y=limit, c='r', ls='--', label= 'Min. density compared')


    # --------------------
    # Panel 1: profiles
    # --------------------
    ax2.plot(
        ds4.potential_density_adjusted[:, glider_indices],
        ds4.oxygen_concentration_corrected[:, glider_indices],
        'k.'
    )
    ax2.plot([], [], 'k.', label='glider profiles')

    ax2.plot(bins, gr1_oxy, '.', c='blue')
    ax2.plot([], [], '.', c='blue', label='glider profiles on density surfaces')

    raw_mask = np.isfinite(comparison.potential_density) & np.isfinite(oxy)

    if comparison_title == 'WOA':
        ax2.errorbar(
            comparison.potential_density[raw_mask],
            oxy[raw_mask],
            yerr=err[raw_mask],
            fmt='.',
            color='blue',
            ecolor='lightblue',
            elinewidth=1,
            capsize=2,
            label=f'Raw {comparison_title} ± uncertainty'
        )
    else:
        ax2.plot(
            comparison.potential_density,
            oxy,
            '.',
            label=f'Raw {comparison_title}'
        )

    iso_mask = np.isfinite(comparison_iso)

    ax2.plot(
        bins[iso_mask],
        comparison_iso[iso_mask],
        '.',
        label=f'{comparison_title}'
    )

    if comparison_title == 'WOA':
        m_err = np.isfinite(comparison_iso) & np.isfinite(err_iso)

        ax2.errorbar(
            bins[m_err],
            comparison_iso[m_err],
            yerr=err_iso[m_err],
            fmt='.',
            color='orange',
            ecolor='red',
            elinewidth=1,
            capsize=2,
            alpha=0.7,
            label=f'{comparison_title}'
        )
        ax2.axvline(x=bins[m_err].min(), c='r', ls='--', label= 'Min. density compared')
        ax1.axhline(y=bins[m_err].min(), c='r', ls='--', label= 'Min. density compared')

    ax2.legend(loc='lower left')
    ax2.set_xlabel('Potential density')
    ax2.set_ylabel('O$_2$ [$\\mu$mol/kg]')
    ax2.set_ylim(0, 350)
    ax2.set_xlim(1023, 1028)
    ax2.grid()



    # --------------------
    # Panel 3: regression
    # --------------------
    x = gr1_oxy
    y = comparison_iso

    ax3.plot(x, y[:, None], '.', c='gray', markersize=4)


    Y = np.repeat(y[:, None], x.shape[1], axis=1)
    BINS = np.repeat(bins[:, None], x.shape[1], axis=1)
    fit_mask = np.isfinite(x) & np.isfinite(Y) & np.isfinite(BINS)

    x_good = x[fit_mask]
    y_good = Y[fit_mask]


    xline = np.arange(0, 350)

    def weighted_fit(c, x_good, y_good, label, weight=None):
        X = sm.add_constant(x_good)

        if weight is None:
            weight = np.ones_like(x_good)

        model = sm.WLS(y_good, X, weights=weight)
        result = model.fit()

        intercept = result.params[0]
        slope = result.params[1]

        intercept_se = result.bse[0]
        slope_se = result.bse[1]

        t_crit = stats.t.ppf(0.975, df=result.df_resid)

        slope_ci = (
            slope - t_crit * slope_se,
            slope + t_crit * slope_se
        )

        icept_ci = (
            intercept - t_crit * intercept_se,
            intercept + t_crit * intercept_se
        )

        xline = np.arange(0, 350)
        yline = intercept + slope * xline
        y_pred = intercept + slope * x_good

        rmse_original = np.sqrt(np.nanmean((y_good - x_good)**2))
        rmse_fit = np.sqrt(np.nanmean((y_good - y_pred)**2))
        ax3.plot(
            xline,
            yline,
            c=c,
            lw=2,
            label=(
            f'{label} Fit\n'
            f'RMSE={rmse_fit:.2f} µmol/kg\n'
            f'y=({slope:.2f} ± {slope_ci[1]-slope:.2f})x + '
            f'{intercept:.0f} ± {icept_ci[1]-intercept:.0f}'
        ))

        ax3.fill_between(
            xline,
            slope_ci[0] * xline + icept_ci[0],
            slope_ci[1] * xline + icept_ci[1],
            color=c,
            alpha=0.2,
            
        )
        #label='95% CI'
        return result,np.round(slope,2), np.round(intercept,0)

    fit1,slope, intercept = weighted_fit('green', x_good, y_good, 'Unweighted')
    depth_weight = np.abs(np.gradient(gr1_dpth, axis=0))
    depth_weight2 = gr1_dpth


    fit_mask = (
        np.isfinite(x) &
        np.isfinite(Y) &
        np.isfinite(depth_weight)
    )

    x_good = x[fit_mask]
    y_good = Y[fit_mask]
    w_good = depth_weight[fit_mask]
    w_good2 = depth_weight2[fit_mask]


    #fit2 = weighted_fit(
    #     'r',
    #    x_good,
    #    y_good, 'Weighted thickness',
    
    #    weight=w_good
    #)

    fit3, slope_weight, intercept_weight = weighted_fit(
        'r',
        x_good,
        y_good, 'Weighted distance from sfc',
        weight=w_good2
    )

    ### force through ML point

    if comparison_title != 'WOA': 
        top_i = np.nanargmin(BINS[fit_mask])
        x0 = x_good[top_i]
        y0 = y_good[top_i]

        dx = x_good - x0
        dy = y_good - y0

        slope_force = np.sum(dx * dy) / np.sum(dx**2)
        intercept_force = y0 - slope_force * x0

        yhat = intercept_force + slope_force * x_good
        resid = y_good - yhat

        df = len(x_good) - 1   # only slope estimated; intercept fixed by anchor
        sigma2 = np.sum(resid**2) / df

        slope_se = np.sqrt(sigma2 / np.sum(dx**2))

        t_crit = stats.t.ppf(0.975, df=df)

        slope_ci = (
            slope_force - t_crit * slope_se,
            slope_force + t_crit * slope_se
        )

        xline = np.arange(0, 350)

        yline = intercept_force + slope_force * xline
        yline_low = y0 + slope_ci[0] * (xline - x0)
        yline_high = y0 + slope_ci[1] * (xline - x0)



    ax3.set_xlabel('Glider O$_2$ on density surfaces [$\\mu$mol/kg]')
    ax3.set_ylabel(f'{comparison_title} O$_2$ on density surfaces [$\\mu$mol/kg]')

    ax3.set_xlim(0, 350)
    ax3.set_ylim(0, 350)

    ax3.grid()

    rmse_11 = np.sqrt(np.nanmean((y_good - x_good)**2))
    ax3.plot(
        [0, 350],
        [0, 350],
        'k--',
        lw=1,
        label=f'1:1 line\nRMSE={rmse_11:.2f}'
    )
    ax3.legend(
        loc='lower right'
    )
    #########
    #Water Mass 
    ##########
    bins2 =  np.arange(1023, 1028.5, 0.1)
    comparison = comparison.interp(depth=ds4.depth)

    if comparison_title == 'WOA': 
        temp_comp = comparison.mean_temperature
    else:
        temp_comp = comparison.temperature  

    mask = (
        np.isfinite(temp_comp) &
        np.isfinite(comparison.potential_density)
    )
    comparison_iso_temp = np.interp(
        bins2,
        comparison.potential_density[mask],
        temp_comp[mask],
        left=np.nan,
        right=np.nan
    )
    gr1_temp = np.full((len(bins2), len(glider_indices)), np.nan)

    c= ['red','green','orange','blue','purple','magenta']
    for i in range(len(glider_indices)):

        rho = ds4.potential_density_adjusted.values[:, glider_indices[i]]
        temp = ds4.temperature.values[:, glider_indices[i]]
        mask = np.isfinite(rho) & np.isfinite(temp)

        rho_valid = rho[mask]
        temp_valid = temp[mask]

        order = np.argsort(rho_valid)

        gr1_temp[:, i] = np.interp(
            bins2,
            rho_valid[order],
            temp_valid[order],
            left=np.nan,
            right=np.nan
        )
    temp_std = np.nanstd(gr1_temp, axis=1)

    ax1.plot(comparison_iso_temp,bins2, 'k',lw=4, label=comparison_title)
    for i in range(len(glider_indices)):
        ax1.plot(gr1_temp[:, i]  , bins2, color=c[i], label=f'index: {glider_indices[i]}')
    ax1.set_ylabel('Potential density [kg/m$^3$]')
    ax1.set_xlabel('Temperature along isopycnals [°C]')
    ax1.legend(loc='upper left')
    ax1.set_ylim(1028,1023)
    ax1.grid()

    if comparison_title != 'WOA':
        diff = np.round(
        np.abs(
                (comparison.time.values - ds4.time.values[glider_indices])
                / np.timedelta64(1, "D")
            ).min(),
            2
        )
        ax1.set_title(f'{comparison_title}, mean {diff} day diff')
    else:
        ax1.set_title(f'{comparison_title}')

    return slope, intercept,  slope_weight ,intercept_weight
In [96]:
fig,axs=plt.subplots(3,1, figsize=(6,10), constrained_layout=True)
bins =  np.arange(1024,1028,0.04)


b1, intercept1, b1_w,intercept1_w = glider_comparison_weighted(axs[0],axs[1],axs[2],np.arange(891,894), best_ctd1, 'CTD1')
# 

labels = ['a', 'b', 'c', 'd', 'e', 'f','g','h','i']

for ax, label in zip(axs.flat, labels):
    ax.text(
        0.02,
        0.98,
        label,
        transform=ax.transAxes,
        fontsize=14,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(
            facecolor='white',
            edgecolor='none',
            alpha=0.6,
            pad=0.2
        )
    )
fig.suptitle(f'{deploy_name}', fontsize=16)
Out[96]:
Text(0.5, 0.98, 'dfo-bb046-20220707')
Figure
No description has been provided for this image

The unweighted fit is a reasonable estimate. The weighted distance from the surface fit is reasonable within error (i.e. gain can't be < 1).

In [97]:
inds = [np.arange(891,894) ]
gains= [b1]
offsets= [intercept1]
gains_w= [1]
offsets_w= [intercept1_w]

label=['CTD1']
comparison = [best_ctd1]

fig,axs=plt.subplots(2,1, figsize=(13,10),sharey=True)
i=0 
oxy_before = ds4.oxygen_concentration_corrected[:,inds[i]]
axs[0].plot(oxy_before,ds4.depth,'k.',markersize=4)
axs[0].plot([],[],c='k',label='glider [O2]')


drift_w = ds4.oxygen_concentration_corrected[:,inds[i]]*gains_w[i] + offsets_w[i]
axs[0].plot( drift_w,ds4.depth,'r.', alpha=0.1)
axs[0].plot([],[],c='r',label=f'Glider [O2] weighted fit: \n [02]*{gains_w[i]} + {offsets_w[i]}')

drift = ds4.oxygen_concentration_corrected[:,inds[i]]*gains[i] + offsets[i]
axs[0].plot( drift,ds4.depth,'g.', alpha=0.1)
axs[0].plot([],[],c='g',label=f'Glider [O2] fit: \n [02]*{gains[i]} + {offsets[i]}')


axs[0].set_xlim(0,350)

axs[0].grid()
axs[0].invert_yaxis()
            
axs[0].plot(comparison[i].oxygen_umol_kg,comparison[i].depth, c='orange')
axs[0].plot([],[],c='orange',label=label[i])
axs[0].invert_yaxis()
### Differences ###
ctd_o2_on_glider_depth = np.interp(
            ds4.depth.values,
            comparison[i].depth.values,
            comparison[i].oxygen_umol_kg,
            left=np.nan,
            right=np.nan
)

axs[0].set_title(f'Drift from {label[i]}')
axs[0].legend(loc='lower right')


diff_before = (oxy_before - ctd_o2_on_glider_depth[:, None])/ctd_o2_on_glider_depth[:, None] *100
diff_after = (drift - ctd_o2_on_glider_depth[:, None])/ctd_o2_on_glider_depth[:, None]*100
diff_after_w = (drift_w - ctd_o2_on_glider_depth[:, None])/ctd_o2_on_glider_depth[:, None] *100


axs[1].plot(diff_before,ds4.depth,'k.',markersize=4)
axs[1].plot([],[],c='k',label=f'glider [O2]')
axs[1].plot( diff_after_w,ds4.depth,'r.', alpha=0.1)
axs[1].plot([],[],c='r',label=f'glider with weighted drift')
axs[1].plot( diff_after,ds4.depth,'g.', alpha=0.1)
axs[1].plot([],[],c='g',label=f'glider with drift')
axs[1].set_xlabel('Percent difference',fontsize=12)   
axs[1].invert_yaxis()
axs[1].legend(loc='lower right')



#ax.set_xlim(1025,1028)
axs[0].set_ylim(1000,0)
axs[1].set_ylim(1000,0)

axs[0].set_ylabel('Depth (m)',fontsize=12)
axs[1].set_ylabel('Depth (m)',fontsize=12)
axs[0].set_xlabel('oxygen [$\mu$mol/kg]',fontsize=12)
axs[1].set_xlabel('Percent difference',fontsize=12)


labels = ['a', 'b', 'c', 'd', 'e', 'f','g','h','i']

for ax, label in zip(axs.flat, labels):
    ax.text(
        0.02,
        0.98,
        label,
        transform=ax.transAxes,
        fontsize=14,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(
            facecolor='white',
            edgecolor='none',
            alpha=0.6,
            pad=0.2
        )
    )
Figure
No description has been provided for this image

We are going to apply the results from the fit from CTD1 since it provides a good fit near the surface

In [98]:
ds4['oxygen_concentration_adjusted'] =  ds4.oxygen_concentration_corrected* b1 + intercept1


    ############
    
fig,axs=plt.subplots(2,2, figsize=(12,6), sharex=True, sharey=True)
N = len(ds4.time)
out = ds4.isel(time=slice(N//2-20,N//2+20))
    
im=axs[0,0].pcolormesh(out.profile,-out.depth,out.oxygen_raw_umol_kg,vmin=0,vmax=300)
min = np.nanmin(out.oxygen_raw_umol_kg)
max= np.nanmax(out.oxygen_raw_umol_kg)
axs[0,0].set_title(f'Raw O2: [{min:.2f}-{max:.2f} µmol/kg]')
    
    
im=axs[0,1].pcolormesh(out.profile,-out.depth,out.oxygen_umol_kg,vmin=0,vmax=300)
min = np.nanmin(out.oxygen_umol_kg)
max= np.nanmax(out.oxygen_umol_kg)
axs[0,1].set_title(f'P & S compensated O2: [{min:.2f}-{max:.2f} µmol/kg]')
    
im=axs[1,0].pcolormesh(out.profile,-out.depth,out.oxygen_concentration_corrected,vmin=0,vmax=300)
min = np.nanmin(out.oxygen_concentration_corrected)
max= np.nanmax(out.oxygen_concentration_corrected)
axs[1,0].set_title(f'Sensor response time corr: [{min:.2f}-{max:.2f} µmol/kg]')
    
im=axs[1,1].pcolormesh(out.profile,-out.depth,out.oxygen_concentration_adjusted,vmin=0,vmax=300)
min = np.nanmin(out.oxygen_concentration_adjusted)
max= np.nanmax(out.oxygen_concentration_adjusted)
axs[1,1].set_title(f'[O2]*{b1} + {intercept1}: [{min:.2f}-{max:.2f} µmol/kg]') 
    
        
axs[0,0].set_ylabel('Depth (m)') 
axs[1,0].set_xlabel('Profile index') 

fig.suptitle(f'{deploy_name}: O2 correction steps')

cbar= fig.colorbar(im,ax=axs.ravel().tolist())


labels = ['a', 'b', 'c', 'd', 'e', 'f','g','h','i']

for ax, label in zip(axs.flat, labels):
    ax.text(
        0.02,
        0.98,
        label,
        transform=ax.transAxes,
        fontsize=14,
        fontweight='bold',
        va='top',
        ha='left',
        bbox=dict(
            facecolor='white',
            edgecolor='none',
            alpha=0.6,
            pad=0.2
        )
    )
Figure
No description has been provided for this image
In [99]:
fig = plt.figure(figsize=(16, 6))

gs = fig.add_gridspec(
    2, 4,
    width_ratios=[0.3, 1, 1, 0.05],
    wspace=0.25,
    hspace=0.28
)

ax_mean = fig.add_subplot(gs[:, 0])

axs = np.empty((2, 2), dtype=object)
axs[0, 0] = fig.add_subplot(gs[0, 1])
axs[0, 1] = fig.add_subplot(gs[0, 2], sharex=axs[0, 0], sharey=axs[0, 0])
axs[1, 0] = fig.add_subplot(gs[1, 1], sharex=axs[0, 0], sharey=axs[0, 0])
axs[1, 1] = fig.add_subplot(gs[1, 2], sharex=axs[0, 0], sharey=axs[0, 0])

cax = fig.add_subplot(gs[:, 3])

raw_mean = (out.oxygen_raw_umol_kg).mean("time")

# Mean profile on left
ax_mean.plot(raw_mean, -out.depth, "k", lw=2)
ax_mean.set_title("Mean raw O2",fontsize=14)
ax_mean.set_xlabel("O2 [µmol/kg]",fontsize=14)
ax_mean.set_ylabel("Depth (m)",fontsize=14)
ax_mean.grid(alpha=0.3)

# Anomaly panels
im = axs[0, 0].pcolormesh(
    out.profile, -out.depth,out.oxygen_raw_umol_kg - raw_mean,
    vmin=-20, vmax=20, cmap=cmap, rasterized=True
)
im.set_rasterized(True)
axs[0, 0].set_title(
    f"Raw O2[µmol/kg]",fontsize=14
)

im = axs[0, 1].pcolormesh(
    out.profile, -out.depth,
    (out.oxygen_umol_kg) - raw_mean,
    vmin=-20, vmax=20, cmap=cmap, rasterized=True
)
im.set_rasterized(True)
axs[0, 1].set_title(
    f"P & S compensated O2 [µmol/kg]",fontsize=14
)

im = axs[1, 0].pcolormesh(
    out.profile, -out.depth,
    (out.oxygen_concentration_corrected) - raw_mean,
    vmin=-20, vmax=20, cmap=cmap, rasterized=True
)
im.set_rasterized(True)
axs[1, 0].set_title(
    f"Sensor response time corr: [µmol/kg]",fontsize=14
)

im = axs[1, 1].pcolormesh(
    out.profile, -out.depth,
    (out.oxygen_concentration_adjusted) - raw_mean,
    vmin=-20, vmax=20, cmap=cmap, rasterized=True
)
im.set_rasterized(True)
axs[1, 1].set_title(
    f"[O2]*{b1:.2f} + {intercept1:.1f} [µmol/kg]",fontsize=14
)

# Axis labels: only left column gets y labels, only bottom row gets x labels
axs[0, 0].set_ylabel("Depth (m)",fontsize=14)
axs[1, 0].set_ylabel("Depth (m)",fontsize=14)
axs[1, 0].yaxis.set_label_coords(-0.1, 0.5)
axs[0, 0].yaxis.set_label_coords(-0.1, 0.5)

axs[1, 0].set_xlabel("Profile index", fontsize=14)
axs[1, 1].set_xlabel("Profile index",fontsize=14)

# Remove top-row x tick labels
for ax in axs[0, :]:
    ax.tick_params(labelbottom=False)

# Remove right-column y labels to prevent overlap/clutter
for ax in axs[:, 1]:
    ax.tick_params(labelleft=False)

# Colorbar in its own axis, so it will not be covered
cbar = fig.colorbar(im, cax=cax)
cbar.set_label("O2 anomaly [µmol/kg]", fontsize=14)
cbar.ax.tick_params(labelsize=13)
fig.suptitle(f"{deploy_name}: O2 anomalies from mean raw O2", y=0.98,fontsize=14)

labels = ["a", "b", "c", "d", "e"]

for ax, label in zip([ax_mean] + list(axs.flat), labels):
    ax.text(
        0.02, 0.98, label,
        transform=ax.transAxes,
        fontsize=14,
        fontweight="bold",
        va="top",
        ha="left",
        bbox=dict(facecolor="white", edgecolor="none", alpha=0.6, pad=0.2)
    )
Figure
No description has been provided for this image

Save and grid data¶

In [100]:
ts4['oxygen_concentration_adjusted'] =  ts4.oxygen_concentration_corrected* b1 + intercept1
ds_final = ts4.drop_vars(['ox_P_S_corr','ox_P_S_corr_smoothed','oxygen_concentration_corrected',
                          'oxygen_umol_kg', 's_factor','p_factor' ])
In [101]:
import json

tau = 3

ds_final["oxygen_concentration_QC"].attrs = ds_final["density_QC"].attrs.copy()

corrections = { "O2 corrections document" : ('https://cproof.uvic.ca/gliderdata/deployments/reports/'),
    "pressure and salinity effects" : ('Applied following manufacturer instructions (see report)'),
    "optode_time_response": (
        f"Applied tau = {tau} s following Bittig et al. (2014)."
    ),
    "boxcar_smoothing": (
        f"Applied {tau/2:g} s boxcar smoothing prior to the optode time response correction."
    ),
    "drift_correction": (
        f"Linear regression using Water Properties CTD cast collected on "
        f"{best_ctd1.time.values} at "
        f"{float(best_ctd1.latitude[0].values):.4f}, "
        f"{float(best_ctd1.longitude[0].values):.4f}; "
        f"slope = {b1:.4f}, intercept = {intercept1:.2f}."
    ),
}

ds_final["oxygen_concentration_adjusted"].attrs = {
    "long_name": (
        "Dissolved oxygen concentration corrected for pressure and salinity effects, "
        "optode time response, converted from μmol L⁻¹ to μmol kg⁻¹, "
        "and corrected for sensor drift"
    ),
    "units": "micromoles_per_kilogram",
    "corrections": json.dumps(corrections),
}
In [102]:
outfile = (f'{filepath}/{deploy_name}_adjusted.nc')
ds_final.to_netcdf(outfile)
ncprocess.make_gridfiles(
        outfile,
        f'{filepath}',
        deployfile,
        fnamesuffix='_adjusted',
        maskfunction=utils.maskQC4,
        max_gap=150
    )
Out[102]:
'deployments/dfo-bb046/dfo-bb046-20220707//dfo-bb046-20220707_grid_adjusted.nc'
In [ ]: