O2 corrections applied to delayed-mode data¶

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

In [1]:
glider_name = 'dfo-eva035'  
deploy_name = f'{glider_name}-20190718'
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'
config_sheet= f'deployments/{glider_name}/AAnderaa_665_configxml.txt' 

description = 'Explorer Seamount'
initials = 'LT'

# O2 specs:
sensor = 'oxygen_0021'

# 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 [2]:
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
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}')

** Explorer Seamount:  glider dfo-eva035**
************
* Deployment: dfo-eva035-20190718
* Sensor: oxygen_0021
make:AROD_FT
seriel:0021

* Processing steps will be saved in: O2_dfo-eva035-20190718.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¶

In [5]:
ds= xr.open_dataset(openfile)

##Apply function to correct for P and S effects. 
ds2 , ds2['s_factor'], ds2['p_factor']= sns.O2_PS_compensation(ds)


### The new vairable ox_P_S_corr is pressure and salinity compensated in umol/L 
ds2_grid = sns.make_gridfiles_adj(ds2, deployfile)
In [6]:
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 [7]:
###There are lots of gaps here too!

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 [8]:
 ### Create oxygen_QC variable. Set all values to QC1

ds2["oxygen_concentration_QC"] = xr.ones_like(
    ds2.oxygen_concentration,
    dtype=np.int8,
)
In [9]:
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[9]:
'deployments/dfo-eva035/dfo-eva035-20190718//dfo-eva035-20190718_grid_O2_2.nc'

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.

In [10]:
### 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')
In [13]:
ts3['oxygen_concentration_adjusted'] =  ts3.oxygen_concentration

ds_final = ts3.drop_vars(['ox_P_S_corr','oxygen_umol_kg', 's_factor','p_factor'])
In [14]:
import json

tau = 25

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)'),
    "Note: ": "There are a lot of gaps in the data. We were not apply to apply the time repsonse correction or drift correction."
}

ds_final["oxygen_concentration_adjusted"].attrs = {
    "long_name": (
        "Dissolved oxygen concentration corrected for pressure and salinity effects, "
        " and converted from μmol L⁻¹ to μmol kg⁻¹, "
       
    ),
    "units": "micromoles_per_kilogram",
    "corrections": json.dumps(corrections),
}
In [15]:
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[15]:
'deployments/dfo-eva035/dfo-eva035-20190718//dfo-eva035-20190718_grid_adjusted.nc'
In [ ]: