O2 corrections applied to delayed-mode data

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

glider_name = 'dfo-bb046'  
deploy_name = f'{glider_name}-20210324'
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_0022'

# 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/')

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

 %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.

# 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-20210324
* Sensor: oxygen_0022
make:AROD_FT
seriel:0022

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

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

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)   

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

##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: 12.22436195470667
max O2: 351.74483338507594
mean O2: 185.92796951699387
std O2: 99.33644648137597

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).

 ### Create oxygen_QC variable. Set all values to QC1

ds2["oxygen_concentration_QC"] = xr.ones_like(
    ds2.oxygen_concentration,
    dtype=np.int8,
)

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
    )

'deployments/dfo-bb046/dfo-bb046-20210324//dfo-bb046-20210324_grid_O2_2.nc'

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

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

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

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.

### 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.

fig,axs= plt.subplots(1,1, figsize=(10,5),sharex=True, sharey=True)
from tqdm import tqdm

### Look at 1 subsets of data 
ranges =[400,460]

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:04<00:00,  1.38tau/s]

Figure

The tau with the lowest objective error function is 0 s, beyond the lowest realistic tau for this sensor according to Bittig 2014. We will use tau=0, the lower limit

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

ts3 = xr.open_dataset(f'{filepath}/{deploy_name}_O2_3.nc')
tau =0
ds_test = sns.apply_tau_correction(ts3, tau)

#####We are going to ensure we are looking at the correct upcast/downcast pairs 

fig,axs=plt.subplots(5)

ranges = np.arange(401,411)
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')

Text(0.5, 0.98, 'Find up/downcast pairs')

Figure

There is a space in the middle when the glider was surfacing.

fig, axs = plt.subplots(1, 5, figsize=(11,8), sharex=True, sharey=True)
ranges = np.arange(401,411)
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

Save data and grid

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
    )

'deployments/dfo-bb046/dfo-bb046-20210324//dfo-bb046-20210324_grid_O2_4.nc'

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

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 = -40, 40
    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[15,:]

    im0 = ax[0].pcolor(
        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].pcolor(
        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

ds_test_grid = xr.open_dataset(f'{filepath}/{deploy_name}_grid_O2_4.nc')
fig = depth_space(ds_test_grid,275,320,tau )

Figure

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.

There is no comparison profile.

We will use the offset and gain found for mission dfo-bb046-20210212, a mission using the same optode.

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

b1 = 1.01
intercept1 = 3

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

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=-70, 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=-70, 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=-70, 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=-70, 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

Save and grid data

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' ])

import json

tau =0 

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"There were no comparison profiles to constrain the linear regression. Constants found for the last mission with a trustworthy comparion profile, dfo-bb046-20210212 were used as drift constants.See dfo-bb046-20210212 report for more information."
        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),
}

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
    )

'deployments/dfo-bb046/dfo-bb046-20210324//dfo-bb046-20210324_grid_adjusted.nc'