CTD corrections applied to delayed-mode data

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

glider_name = 'dfo-bb046'
deploy_name = f'{glider_name}-20200908'
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}/L0-timeseries/{deploy_name}_delayed.nc'
opengridfile = f'{filepath}/L0-gridfiles/{deploy_name}_grid_delayed.nc'
deployfile = f'{filepath}/deployment.yml'

description = 'Calvert'
initials = 'LT'

# CTD specs:
sensor = 'CTD_0278'

# 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 = 'C-PROOF_SBE_CTDProcessingReport_v0.2.pdf'
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/')

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

# print(f'* Protocols are detailed in: {processing_protocol}')
print(f'* Processing steps will be saved in: CTD_{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-20200908
* Sensor: CTD_0278

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

display(Markdown("./docs/CTD_1_Preamble.md"))

1.0 Preamble

This document describes conductivity, temperature, and pressure data processing steps applied to delayed mode data collected using Sea-Bird Scientific Glider Payload Conductivity Temperature Depth (GPCTD) sensors mounted on C-PROOF Slocum and SeaExplorer autonomous ocean gliders. This sensor has a nominal sampling rate of 1 Hz and was designed specifically for Slocum gliders. This document covers the application of the sensor alignment correction and the thermal lag correction, as well as removal of questionable conductivity values and salinity profiles.

1.1 Set up the processing

The processing steps below are applied to delayed mode data stored in a single netcdf timeseries file created using the Pyglider data processing package (https://github.com/c-proof/pyglider).

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.

import warnings
warnings.filterwarnings('ignore')

import xarray as xr
import numpy as np
import pathlib
import pyglidersensor as pgs
from pyglider.ncprocess import make_gridfiles

from datetime import datetime, date
%matplotlib ipympl


import scipy.stats as stats

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

%reload_ext autoreload
%matplotlib ipympl

1.2 Profile Check

Check that upcasts and downcasts are being properly identified. Negative values should be associated with upcasts.

print(f'Loading: {openfile}')

caption = ('Identifying upcasts and downcasts. The left panel shows '
           'pressure vs. time and the right panel shows profile direction vs. '
           'time for a small subset of the time series:')


fname = openfile

with xr.open_dataset(fname) as ds0:
    # SAVE SOME PARAMS FOR PLOTTING DOWN BELOW!
    N = len(ds0.time)
    MAX_DEPTH = np.nanmax(ds0.depth)
    NUM_PROFILES = np.nanmax(ds0.profile_index)
    
    print('************')
    print(f'* There are {N} data points in total, with {NUM_PROFILES} profiles')
    print(f'* Time period: {pd.to_datetime(np.nanmin(ds0.time)).strftime("%Y-%m-%d")} to {pd.to_datetime(np.nanmax(ds0.time)).strftime("%Y-%m-%d")}')
    print(f'* Depth range: {round(np.nanmin(ds0.depth))} - {round(MAX_DEPTH)} metres')
    print('************')
    if N > 50000:
        todo = slice(int(N/2)-5000, int(N/2)+5000)
    else:
        todo = slice(int(N/3), int(2*N/3))
    
    fig, axs = plt.subplots(nrows=1, ncols=2, 
                            constrained_layout=True, 
                            figsize=(9, 4))

    ds = ds0.isel(time=todo)
    axs[0].plot(ds.time, ds.pressure, '.', markersize=1)
    axs[0].set_ylim([MAX_DEPTH, 0])
    axs[0].set_ylabel('Pressure [dbar]')
    axs[0].tick_params(axis='both', labelsize=8)
    axs[0].grid(axis='x')

    axs[1].plot(ds.time, ds.profile_direction, '.', markersize=1)
    axs[1].set_ylabel('Profile Direction')
    axs[1].tick_params(axis='both', labelsize=8)
    axs[1].grid(axis='x') 
    print(caption)

Loading: deployments/dfo-bb046/dfo-bb046-20200908//L0-timeseries/dfo-bb046-20200908_delayed.nc
************
* There are 1929203 data points in total, with 913.0 profiles
* Time period: 2020-09-08 to 2020-09-30
* Depth range: 0 - 992 metres
************
Identifying upcasts and downcasts. The left panel shows pressure vs. time and the right panel shows profile direction vs. time for a small subset of the time series:

Figure

1.3 Delayed-mode data prior to corrections

Checking fields (temperature, salinity, conductivity and density) in the delayed-mode data, before any CTD corrections:

tds = opengridfile
ds = xr.open_dataset(tds)
list(ds.keys())

fig, axs = plt.subplots(4, 1, figsize=(11, 10), sharey=True, sharex=True)

xlims = [0, NUM_PROFILES]
ylims=[MAX_DEPTH,0]
# ylims=[50,0]

pc = axs[0 ].pcolormesh(ds.profile, ds.depth, ds['salinity'],rasterized=True)
axs[0].set_ylim(ylims)
axs[1].set_xlim(xlims)
fig.colorbar(pc, ax=axs[0], label = 'Salinity [psu]')
axs[0].set_title('Salinity',loc='left')

pc = axs[1].pcolormesh(ds.profile, ds.depth, ds['temperature'],rasterized=True,cmap='plasma')

fig.colorbar(pc, ax=axs[1], label = 'Temperature [$^o$C]')
axs[1].set_title('Temperature',loc='left')
pc = axs[2].pcolormesh(ds.profile, ds.depth, ds['conductivity'],rasterized=True,cmap='cividis')
fig.colorbar(pc, ax=axs[2], label = 'Conductivity [S/m]')
axs[2].set_title('Conductivity',loc='left')

# pc = axs[3].pcolormesh(ds.profile, ds.depth, ds['oxygen_concentration'],rasterized=True,cmap='inferno')
# fig.colorbar(pc, ax=axs[3])
# axs[3].set_title('Oxygen Concentration',loc='left')

pc = axs[3].pcolormesh(ds.profile, ds.depth, ds['density'],rasterized=True,cmap='inferno')
fig.colorbar(pc, ax=axs[3], label = 'Density [kg/m$^3$]')
axs[3].set_title('Density',loc='left')

axs[0].set_ylabel('Depth [m]');
axs[1].set_ylabel('Depth [m]');
axs[2].set_ylabel('Depth [m]');
axs[3].set_ylabel('Depth [m]');

Figure

display(Markdown("./docs/CTD_2_Steps.md"))

2.0 Corrections applied to delayed mode data for this mission

Processing steps:

1. Identification of anomalous conductivity values
2. Identification of questionable temperature profiles
3. Sensor alignment correction
4. Thermal lag correction

Apply QC flags to data following Argo notation

1 : Good data
3 : Bad data, potentially correctable (e.g. sensor clogged)
4 : Bad data (e.g. sensor failed)
8 : Estimated value (interpolated, extrapolated or other estimation)

2.1 Identification and flagging anomalous conductivity values

We identify and flag any conductivity values that are obviously unphysical, which is typically caused by air bubbles in the conductivity cell. We use a simple criterion applied to the raw conductivity data. The criterion temporarily flags any data points that are more than 5 standard deviations away from the overall time series mean for a given depth bin and profile bin, then recomputes the mean and standard deviation, excluding the temporarily flagged values. Conductivity values that still differ from the mean by more than 3 standard deviations are flagged as 'bad' (QC 4). If the difference between the 'bad' values and the mean is less than the accuracy of the sensor, which is 0.0003 S/m for the GPCTD, then those points are not excluded.

This criterion is applied to data binned first by profile index, in increments of 50 profiles, then binned by depth, in increments of 5 m. The use of profile index bins rather than time or temperature bins is designed to allow for the removal of unphysical values.

Adjustments to this correction are based on examining the data and making a judgment call about which conductivity values are undeniably 'bad'. In this case, we want to exclude the extremely low values occurring at the surface consistent with air bubbles in the cell. Some unphysical values are missed by this correction, and may be caught during the removal of unphysical salinity profiles in further stepsbelow.

Note that for this mission:

srate = stats.mode((np.diff(ds0.time)).astype('timedelta64[s]')).mode
fs = 1/srate.astype(float) #the sampling frequency = 1/(delta t)
print('************')
print(f'The mode of the sampling rate for the GPCTD is one sample every {srate}.')
print('************')

************
The mode of the sampling rate for the GPCTD is one sample every 1 seconds.
************

# Identify the questionable conductivity values
flag_stdev = 5 #number of standard deviations to temporarily flag bad salinity values 
clean_stdev = 3 #number of standard deviations to flag bad conductivity values, after removing the temporary bad values from the calc
dT = 50 #size of the profile bins
dz = 5 #size of the depth bins

ts0 = ds0.copy() 
ts0.conductivity[ts0.conductivity<0.1] = np.nan
ts = pgs.get_conductivity_clean(ts0, dT, dz, flag_stdev, clean_stdev, accuracy)


##### make new variable for conductivity QC 

ts = ts.assign(conductivity_QC = np.nan*ts.conductivity)
ts['conductivity_QC'] = xr.where(np.isfinite(ts.conductivityClean), 1, 4)

#####

# Figures to look at the comparison    
fig, ax = plt.subplots(1,2,figsize=(10,4), constrained_layout=True)

ax[0].plot(ts.conductivity, ts.temperature, color='r', marker='.', linestyle='none', label='QC 4')
ax[0].plot(ts.conductivityClean, ts.temperature, color='k', marker='.', linestyle='none', label='QC 1')
ax[0].set_ylabel('Temperature [$^o$C]', fontsize=16)
ax[0].set_xlabel('Conductivity [S/m]', fontsize=16)
ax[0].grid(axis='both', color='0.5')

ax[1].plot(ts.profile_index, ts.conductivity, color='r', marker='.', linestyle='none')
ax[1].plot(ts.profile_index, ts.conductivityClean, color='k', marker='.', linestyle='none')
ax[1].set_xlabel('Profile index', fontsize=16)
ax[1].set_ylabel('Conductivity [S/m]', fontsize=16)
ax[1].grid(axis='both', color='0.5')

ax[0].legend()

print('Fig 2: Temperature vs. conductivity (left), depth vs. conductivity (middle), '
      'and conductivity vs. profile index (right), '
      'with the red dots showing the unphysical values flagged as bad and are flagged as QC 4:')

Fig 2: Temperature vs. conductivity (left), depth vs. conductivity (middle), and conductivity vs. profile index (right), with the red dots showing the unphysical values flagged as bad and are flagged as QC 4:

Figure

Adjustments to this correction are based on examining the data and making a judgment call about which conductivity values are undeniably 'bad'. In this case, we want to flag the extremely low values occurring at the surface (Fig. 2) consistent with air bubbles in the cell. Some unphysical values are missed by this correction, and may be caught during the removal of unphysical profiles below.

###### grid to make finding bad profiles easier
ts.to_netcdf(f'{filepath}/{deploy_name}_QC.nc')
# Save a gridded version as well
outfile = make_gridfiles(f'{filepath}/{deploy_name}_QC.nc', f'{filepath}', deployfile, fnamesuffix='QC') 

ts = xr.open_dataset(f'{filepath}/{deploy_name}_QC.nc')
ds=xr.open_dataset(f'{filepath}/{deploy_name}_gridQC.nc') 


# Now adding zoomed plot
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 = [300, 0] #######modified 

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

##################
ds2 = ds.where(ds.conductivity_QC != 4) 
ds0 = xr.open_dataset(tds)
cmap='jet'
vmin=3.2
vmax = 3.75

profile_lims = xlim_1
#ds_sub = ds_sub.isel(depth=range(200,800))

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

ax = axs[0,1]
ds_sub = ds2.where((ds2.profile_index >=profile_lims[0]) & (ds2.profile_index <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, ds_sub.depth, ds_sub['conductivity'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
axs[0,0].set_ylim(Y_LIMS)
axs[0,0].set_title('Conductivity',loc='left')
fig.colorbar(pc, ax=ax, label = 'Conductivity ')



########
profile_lims = xlim_2

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

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

fig.colorbar(pc, ax=ax, label = 'Conductivity')

######
profile_lims = xlim_3

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

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

fig.colorbar(pc, ax=ax, label = 'Conductivity')

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

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

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

fig.colorbar(pc, ax=ax, label = 'Conductivity')

 
print('Zooming in along the glider deployment to visualize conductivity spikes.')
 

Zooming in along the glider deployment to visualize conductivity spikes.

Figure

There are stripes from profiles 480 - 487 may indicate a clogged CTD.

Identify sections where the conductivity cell is clogged

Three ranges were manually selected to see TS patterns in up and downcasts.

fig, axs = plt.subplots(1,3, figsize=(10,4))
inds = [ range(84,90),  range(480,487), range(750,755)]


for j in range(3):
    for ind in inds[j]:    
        ds2 = ds.where((ds.profile_index==ind),drop=True)
        t = ds2.temperature
        s = ds2.salinity

        ax = axs[j]
        if ds2.profile_direction.mean() < 0:
            ax.plot(t, s, 'g',marker='.', label = 'Downcasts' if ind== 84 else None)
        elif ds2.profile_direction.mean() > 0:
            ax.plot(t, s, 'r', marker='.',label = 'Upcast' if ind== 85 else None)
    ax.set_title(str(inds[j]))


axs[0].legend()
axs[0].set_xlabel('T [C]')
axs[0].set_ylabel('Salinity [psu]')

Text(0, 0.5, 'Salinity [psu]')

Figure

Between profiles 84-90 and 750–755, the conductivity sections show a regular banded pattern. Over these ranges, the T–S diagrams indicate that upcasts and downcasts follow nearly identical curves, which is consistent with a thermal-lag artifact rather than a change in water mass properties. In contrast, between profiles 480–487 the conductivity pcolor exhibits irregular stripes and the upcast and downcast T–S curves diverge substantially in shape. This behaviour suggests a sensor problem (e.g., clogging of the conductivity cell) rather than a simple thermal-lag mismatch. We will manually set values to QC 4 in this range.

Apply flag QC 4 to conductivty and salinity in this profile range.

ds = xr.open_dataset(f'{filepath}/{deploy_name}_gridQC.nc')

prof_range = np.arange(480, 487)

ts['conductivity_QC'] = xr.where(ts.profile_index.isin(prof_range),4, ts['conductivity_QC'])
ds['conductivity_QC'] = xr.where(ds.profile_index.isin(prof_range),4, ds['conductivity_QC'])

ds['conductivity_adjusted'] = ds['conductivity']

ts['conductivity_adjusted'] =ts['conductivity']


ds['conductivity_adjusted_QC'] = ds['conductivity_QC']

ts['conductivity_adjusted_QC'] = ts['conductivity_QC']

fig,axs=plt.subplots(2, sharex=True,sharey=True)

vmin=0
vmax = 4.4
d4 = ds.where(ds.conductivity_QC!=4)

pc = axs[0].pcolormesh(d4.profile, -d4.depth, d4['conductivity'],rasterized=True,cmap='cividis',vmin=vmin,vmax=vmax
)
fig.colorbar(pc, label = 'Conductivity, no QC4 [S/m]')

pc = axs[1].pcolormesh(d4.profile, -d4.depth, d4['conductivity_adjusted'],rasterized=True,cmap='cividis',vmin=vmin,vmax=vmax
)
fig.colorbar(pc, label = 'Conductivity [S/m]')
axs[0].set_title('Conductivity, no QC4',loc='left')

Text(0.0, 1.0, 'Conductivity, no QC4')

Figure

2.2 Identifying questionable temperature profiles

To identify potentially suspicious temperature profiles, we flag measurements that indicate sensor malfunction or biofouling. Temperature data are assigned a QC flag of 4 when any of the following conditions occur:

Temperature spikes: individual points where the temperature differs from the average of adjacent samples by more than 0.75 °C, typically indicative of biofouling or flow cell obstruction.

Inter-profile rate of change: differences in temperature between adjacent profiles at the same depth that exceed 1.5 °C, suggesting that the entire profile may be biased (e.g., due to sensor lag or partial clogging), rather than reflecting true environmental variability.

Density inversions: occurrences where the density of the deeper depth bin is lower than that of the shallower bin, indicating physically unrealistic structure. Density inversions are only considered below 100 m, where stratification should be stable; inversions above this depth are ignored because mixed-layer turbulence frequently produces real, short-lived overturns.

##spike algorithum - identifies isolated spikes from biofouling or clogged cells 

temp = ts.temperature

temp_fwd = temp.shift(time=-1)
temp_bwd = temp.shift(time=1)
spike = np.abs(temp-(temp_fwd+temp_bwd)/2) > 0.75

ts['temperature_QC'] = xr.where(spike,4,1)

###### grid to make finding bad profiles easier
ts.to_netcdf(f'{filepath}/{deploy_name}_QC2.nc')
# Save a gridded version as well
outfile = make_gridfiles(f'{filepath}/{deploy_name}_QC2.nc', f'{filepath}', deployfile, fnamesuffix='QC2') 

## rate of change between profiles - catches jumps where the whole profile is off 
ts = xr.open_dataset(f'{filepath}/{deploy_name}_QC2.nc')
ds=xr.open_dataset(f'{filepath}/{deploy_name}_gridQC2.nc') 


dtemp = ds.temperature.diff(dim='time')
dtemp = dtemp.reindex(time=ds.time,fill_value=0)
rate_change = np.abs(dtemp) > 1.5 #large not to include fronts


###### density inversion - misaligned T/S within a profile, something wrong with cell 

SA = gsw.SA_from_SP(ds.salinity, ds.longitude, ds.latitude, ds.depth)
CT = gsw.CT_from_t(SA, ds.temperature, ds.depth)
rho = gsw.rho(SA, CT, ds.depth)
rho_fwd = rho.shift(depth=-1)
density_inversion = rho_fwd < rho

# Mask out the upper 75 m — don't allow inversions to count there
density_inversion = xr.where(ds.depth > 100, density_inversion, False)

ds['temperature_QC'] = xr.where(rate_change,4,ds['temperature_QC'])
ds['temperature_QC'] = xr.where(density_inversion,4,ds['temperature_QC'])

### We need to make sure the flag is being applied to the timeseries as well. These corrections were easier to apply to the profiles

#Find QC4 positions in grid
qc4_positions = (ds.temperature_QC == 4)

#convert to integer for interpolation
qc4_int = qc4_positions.astype(int)

#Interpolate mask onto timeseries (time, raw depth)
qc4_ts_int = qc4_int.interp(
    time=ts.time,
    depth=ts.depth,
    method="nearest"
)
qc4_ts = (qc4_ts_int == 1)

ts['temperature_QC'] = xr.where(qc4_ts, 4, ts['temperature_QC'])

fig, ax = plt.subplots(1,2,figsize=(9,4), 
                           constrained_layout=True)
t4 = ts.where((ts.temperature_QC == 4) | (ts.conductivity_QC == 4), drop=True) 

ax[0].plot(ts.salinity, ts.temperature, color='k', marker='.', linestyle='none', label='QC 1')
ax[0].plot(t4.salinity, t4.temperature, color='r', marker='.', linestyle='none', label='QC 4')

ax[0].set_ylabel('Temperature [$^o$C]', fontsize=12)
ax[0].set_xlabel('Salinity [psu]', fontsize=12)
ax[0].grid(axis='both', color='0.5')
ax[0].legend()


#####
bad_profiles = ts.profile_index.where(ts['temperature_QC']==4, drop=True)
t4 = ts.where((ts.temperature_QC == 4) , drop=True) 

ax[1].plot(ts.profile_index,ts.temperature,marker='.', linestyle='none',  c='k') 
ax[1].plot(t4.profile_index,t4.temperature,marker='.', linestyle='none', c='r') 

ax[1].set_ylabel('Temperature  [$^o$C]', fontsize=12)
ax[1].set_xlabel('Profile index', fontsize=12)
ax[1].grid(axis='both', color='0.5')

print ('Salinity plotted as a function of temperature '  
'(left) and vs. profile index (right), \n with the salinity profiles ' 
'flagged as QC 4 shown in red and indicated by '
'the red arrows at the top of the panel on the right:')

Salinity plotted as a function of temperature (left) and vs. profile index (right), 
 with the salinity profiles flagged as QC 4 shown in red and indicated by the red arrows at the top of the panel on the right:

Figure

ts.to_netcdf(f'{filepath}/{deploy_name}_QC3.nc')
ds.to_netcdf(f'{filepath}/{deploy_name}_gridQC3.nc')

Manually identify questionable profiles

ts = xr.open_dataset(f'{filepath}/{deploy_name}_QC3.nc')
ds=xr.open_dataset(f'{filepath}/{deploy_name}_gridQC3.nc') 

ds0 = xr.open_dataset(tds)

# Now adding zoomed plot
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]
vmin=4
vmax = 13


Y_LIMS = [800, 0] #######modified 

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

##################
ds4 = ds.where(ds.temperature_QC != 4) 

profile_lims = xlim_1
#ds_sub = ds_sub.isel(depth=range(200,800))

ax = axs[0,0]
ds_sub = ds0.where((ds0.profile >=profile_lims[0]) & (ds0.profile <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, ds_sub.depth, ds_sub['temperature'],rasterized=True,vmin=vmin,vmax=vmax)

ax = axs[0,1]
ds_sub = ds4.where((ds4.profile >=profile_lims[0]) & (ds4.profile <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, ds_sub.depth, ds_sub['temperature'],rasterized=True,vmin=vmin,vmax=vmax)
axs[0,0].set_ylim(Y_LIMS)
axs[0,0].set_title('Temperature: Adjusted',loc='left')
fig.colorbar(pc, ax=ax, label = 'Temperature [c]')



########
profile_lims = xlim_2

ax = axs[1,0]
ds_sub = ds0.where((ds0.profile >=profile_lims[0]) & (ds0.profile <= profile_lims[1]), drop=True)
pc = ax.pcolormesh(ds_sub.profile, ds_sub.depth, ds_sub['temperature'],rasterized=True,vmin=vmin,vmax=vmax)

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

fig.colorbar(pc, ax=ax, label = 'Temperature [c]')

######
profile_lims = xlim_3

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

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

fig.colorbar(pc, ax=ax, label = 'Temperature [c]')

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

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

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

fig.colorbar(pc, ax=ax, label = 'Temperature [c]')

 
print('Zooming in along the glider deployment to visualize temperature spikes. Red arrows identify some spikes to be manually flagged.')
 

Zooming in along the glider deployment to visualize temperature spikes. Red arrows identify some spikes to be manually flagged.

Figure

Temperature also seems clogged during profiles 480-487, like we saw in the conductivity signal. We will manually flag these as QC 4. We do not see the same spike at profile 790, so we will keep the QC 1 flg.

ts['temperature_QC'] = xr.where(ts.profile_index.isin(prof_range),4, ts['temperature_QC'])

ds['temperature_QC'] = xr.where(ds.profile_index.isin(prof_range),4, ds['temperature_QC'])

ts['temperature_adjusted'] =ts.temperature 
ds['temperature_adjusted'] =ds.temperature 

ds['temperature_adjusted_QC'] = ds['temperature_QC']  

ts['temperature_adjusted_QC'] = ts['temperature_QC']  

fig,axs=plt.subplots(2, sharex=True,sharey=True)

vmin=3.2
vmax = 15

ds4 = ds.where(ds.temperature_QC!=4)
pc = axs[0].pcolormesh(ds4.profile, -ds4.depth, ds4['temperature'],rasterized=True,cmap='plasma',vmin=vmin,vmax=vmax
)
fig.colorbar(pc, label = 'Temperature [S/m]')

ds4 = ds.where(ds.temperature_adjusted_QC!=4)
pc = axs[1].pcolormesh(ds4.profile, -ds4.depth, ds4['temperature_adjusted'],rasterized=True,cmap='plasma',vmin=vmin,vmax=vmax
)
fig.colorbar(pc, label = 'Temperature [S/m]')
 

<matplotlib.colorbar.Colorbar at 0x3520ddd10>

Figure

Calculate a new salinity variable using the smoothed conductivty values. Like conductivity, the resolution of QC 8 data is 9 km.

ds['salinity_adjusted'] = xr.DataArray(
    sw.salt(10 * ds.conductivity_adjusted / sw.constants.c3515,
            ds.temperature_adjusted,
            ds.pressure),
    dims=('depth', 'time'),
    coords={'time': ds.time, 'depth': ds.depth}
)
ts['salinity_adjusted'] = xr.DataArray(
    sw.salt(10 * ts.conductivity_adjusted / sw.constants.c3515,
            ts.temperature_adjusted,
            ts.pressure),
    dims=('time'),
    coords={'time': ts.time}
)

Apply salinity flags

The worst QC flag assigned to either conductivity or temperature is also applied to salinity, since salinity is derived from both variables. (QC 4 > QC 3 > QC 1)

sal_QC = xr.where(
    (ds.temperature_QC == 4) | (ds.conductivity_QC == 4),
    4,
    xr.where(
        (ds.temperature_QC == 3) | (ds.conductivity_QC == 3),
        3,
        1,
    )
)

ds['salinity_QC'] = sal_QC

sal_QC = xr.where(
    (ds.temperature_adjusted_QC == 4) | (ds.conductivity_adjusted_QC == 4),
    4,
    xr.where(
        (ds.temperature_adjusted_QC == 3) | (ds.conductivity_adjusted_QC == 3),
        3,
        1,
    )
)

ds['salinity_adjusted_QC'] = sal_QC

###############
sal_QC = xr.where(
    (ts.temperature_QC == 4) | (ts.conductivity_QC == 4),
    4,
    xr.where(
        (ts.temperature_QC == 3) | (ts.conductivity_QC == 3),
        3,
        1,
    )
)

ts['salinity_QC'] = sal_QC

sal_QC = xr.where(
    (ts.temperature_adjusted_QC == 4) | (ts.conductivity_adjusted_QC == 4),
    4,
    xr.where(
        (ts.temperature_adjusted_QC == 3) | (ts.conductivity_adjusted_QC == 3),
        3,
        1,
    )
)
ts['salinity_adjusted_QC'] = sal_QC

fig,axs=plt.subplots(2, sharex=True,sharey=True)

vmin=30
vmax = 35
ds4 = ds.where(ds.salinity_adjusted_QC!=4)

pc = axs[0].pcolormesh(ds4.profile, -ds4.depth, ds4['salinity'],rasterized=True,vmin=vmin,vmax=vmax)
axs[0].set_title('Salinity, no QC4') 
fig.colorbar(pc, label = 'Salinity [S/m]')

pc = axs[1].pcolormesh(ds4.profile, -ds4.depth, ds4['salinity_adjusted'],rasterized=True,vmin=vmin,vmax=vmax)
fig.colorbar(pc, label = 'Salinity [S/m]')

<matplotlib.colorbar.Colorbar at 0x35222bed0>

Figure

# T-S diagram to select near-surface water density range to exclude

srate = stats.mode((np.diff(ts.time)).astype('timedelta64[s]')).mode
fig,ax=plt.subplots()

ax.plot(ts.salinity,ts.temperature,'k.',markersize=2, label = 'Delayed-mode')

ts1 = ts.where((ts.salinity_QC!=4 )&(ts.temperature_QC!=4))
ax.plot(ts1.salinity_adjusted, ts1.temperature_adjusted, '.', markersize = 2, label = 'Q4 salinity profiles removed')

#Create a density grid to contour plot isopycnals
S_range = np.linspace(int(np.min(ts.salinity)-0.5), 
                      int(np.max(ts.salinity)+0.5), 1000)
T_range = np.linspace(int(np.min(ts.temperature)-1), 
                      int(np.max(ts.temperature)+1), 1000)
S_grid, T_grid = np.meshgrid(S_range, T_range)
density_grid = seawater.eos80.dens0(S_grid, T_grid)

CS = ax.contour(S_range, T_range, density_grid,
                np.arange(1014,
                          np.round(np.max(density_grid)),0.5),
                colors='k', linewidths=0.5);
ax.clabel(CS, CS.levels, inline=True, fontsize=10)
ax.set_xlabel('Salinity [psu]')
ax.set_ylabel('Temperature [$^o$C]')
# plt.xlim(28,35)
ax.grid()
ax.legend(prop={'size': 10})
print('Temperature vs. salinity diagram. ',
     'Black contours give density in kg/m^3:')

Temperature vs. salinity diagram.  Black contours give density in kg/m^3:

Figure

###### grid to make finding bad profiles easier
ts.to_netcdf(f'{filepath}/{deploy_name}_QC4.nc')
# Save a gridded version as well
ds.to_netcdf(f'{filepath}/{deploy_name}_gridQC4.nc')

display(Markdown('./docs/CTD_2_Sensor_lag_LT.md'))

2.3 Sensor alignment correction

We now test application of a sensor alignment correction. In the literature this correction is often used to align the temperature and conductivity in time, relative to the pressure. This correction reduces the occurrence of salinity spikes near sharp gradients in T and S, and ensures calculations are made using the same parcel of water for all variables. The misalignment between the sensors is caused by:

1. The physical separation between sensors causing a transit time delay for water being pumped through the CTD, and,
2. Different sensor response times

We follow the SeaBird Electronics Data Processing Manual (page 80) to determine if there is any time lag between the temperature and conductivity sensors on our pumped CTD.

Sources: Sea-Bird Electronics, Inc. SEASOFT V2: SBE Data Processing (https://misclab.umeoce.maine.edu/ftp/instruments/CTD%2037SI%20June%202011%20disk/website/pdf_documents/manuals/SBEDataProcessing_7.21d.pdf)

fig, ax = plt.subplots(figsize=(10,5))

with xr.open_dataset(f'{filepath}/{deploy_name}_QC4.nc') as ts:
    N = len(ts.time)
    
    if N > 50000:
        todo = slice(int(N/2)-1000, int(N/2)+1000)
    else:
        todo = slice(int(N/4), int(2*N/4))
    
    ts=ts.isel(time=todo)
    ts1 = ts.where(ts.conductivity_QC!=4)
    
    ax.plot(ts1.time,1.25*ts1.conductivity_adjusted, c='red', label='1.25xConductivity')
    ax.set_ylabel('Conductivity/Temperature',)
    ax.set_xlabel('Time',)

    ax.plot(ts1.time,ts1.temperature_adjusted ,c='blue',label='Temperature')
    ax.legend()
    ax.grid()
    
    print('Plot timeseries of temperature and Q1 conductivity to observe any time offset.')

Plot timeseries of temperature and Q1 conductivity to observe any time offset.

Figure

display(Markdown('./docs/CTD_2_Thermal_lag.md'))

2.4 Thermal lag correction

Thermal lag arises from the finite thermal mass of the conductivity cell, which alters the temperature of water passing through it. We follow the method of Garau et al. (2011), building upon Morison et al. (1994). Because gliders sample continuously at a stable rate, thermal-lag behaviour can be considered constant throughout a mission.

A recursive filter estimates the temperature inside the conductivity cell:

\begin{equation} T_T(n) = -b\,T_T(n-1) + a\,T(n) - a\,T(n-1), \end{equation} where \begin{equation} a = \frac{4 f_n \alpha \tau}{1 + 4 f_n \tau}, \qquad b = 1 - \frac{2a}{\alpha}. \end{equation} Here,

$\tau$ is the thermal-lag time constant (s),

$\alpha$ is the lag strength, and

$f_n$ is the sampling frequency.

The estimated temperature within the conductivity cell is:

\begin{equation} T_c(n) = T(n) - T_T(n). \end{equation}

Previous C-PROOF processing used fixed parameters $\alpha = 0.06$ and $\tau = 10$ s, values determined by Janzen and Creed (2011). Thermal-lag parameters typically vary slightly among individual GPCTD sensors—more so than alignment-correction constants. Once an optimal $\alpha$ is identified for a particular sensor, subsequent missions are first tested with those parameters, and the correction is applied only if it improves the dive–climb agreement. Otherwise, the parameters are re-evaluated.

Optimal values are found by applying the recursive filter across a grid of $\alpha$ and $\tau$ values and selecting the combination that minimizes the root-mean-squared difference (RMSD) between paired dive–climb salinity profiles. RMSD is computed as the square root of the summed squared area between pairs of salinity profiles (binned by temperature), normalized by the number of profile pairs. Following Garau et al. (2011), 20 profile pairs, evenly spaced in time, are used for the optimization.

\subsubsection*{Applying the thermal-lag correction to salinity}

Corrected salinity is computed using the newly estimated cell temperature $T_c(n)$ together with the conductivity_adjusted variable. The correction is applied only to data flagged as QC 1; regions flagged as QC 8 retain their smoothed values.

# Set up our constants

density_cutoff = 1023 #this is not excluding the top of profiles (exclude everything less dense than this from the minimization)
num_profs = 100 #number of profiles to include in the subset of data # This is not used for the correction, but is used in the q/c steps
clean_profs_start = 50 #number of profiles to exclude from the start
clean_profs_end = 0 #number of profiles to exclude from the end
dn_stdev = 1 #how many standard deviations from the mean the area between downcasts can be

# Load time series
ts = xr.open_dataset(f'{filepath}/{deploy_name}_QC4.nc')
# Save a gridded version as well
ds = xr.open_dataset(f'{filepath}/{deploy_name}_gridQC4.nc')

srate = stats.mode((np.diff(ds.time)).astype('timedelta64[s]')).mode
fs = 1/srate.astype(float) 
fn = 0.5*fs #frequency for Sea-Bird GPCTD

ts = ts.assign_coords(pind=ts.profile_index) #add a profile index coordinate
tot_profs = int(np.nanmax(ts.profile_index.values))
print('Total number of profiles:', tot_profs)

####keep values without Q4 data 
ds1 = ds.where(ds.salinity_QC!=4) 
ts1 = ts.where(ts.salinity_QC!=4) 

Total number of profiles: 913

fig, ax = plt.subplots()
vmin=31
vmax=34
pc=ax.pcolormesh(ds1.profile, -ds1.depth, ds1['salinity_adjusted'],rasterized=True,vmin=vmin,vmax=vmax)
fig.colorbar(pc,label='Salinity [psu]')
ax.set_title('Q1 salinity up/down asymmetry')
ax.set_xlim(100,140)
ax.set_ylim(-300,0)


print('There is a visible asymmetry in the up and down casts of succesive profiles, causing stipes to appear in the Q1 data' )

There is a visible asymmetry in the up and down casts of succesive profiles, causing stipes to appear in the Q1 data

Figure

We exlude the first 50 profiles, which were primarily collected in the highly energetic environment on the shelf near the shelf. We use a subet of the remianing data, consisting of 20 pairs of profiles equally spaced in time to determine the corrected.

print('Calculating profile pairs')


ts_sub, profile_bins, profile_bins_all, direction = pgs.profile_pairs(
    ts1, clean_profs_start, clean_profs_end, num_profs, prof_range
) 

# Identify boolean index for application of density cutoff 
density_bool = ts_sub.density>=density_cutoff

Calculating profile pairs

Restricting profiles > 1 standard deviation greater than the mean space between profiles

#Determine the area between subsequent downcasts to restrict profiles included in correction 
print(f'Restricting profiles')
dn_area, area_bad = pgs.TS_preprocess(
    density_bool, dn_stdev, 
    profile_bins, profile_bins_all, 
    direction, ts_sub)
print('Max and min area between downcasts = ', np.nanmax(dn_area), np.nanmin(dn_area))

ts_bad = ts_sub.where(
    ts_sub.profile_index.isin(
        profile_bins_all[area_bad]), 
        drop=True)
prof_list = ts_bad.profile_index
print('List of profiles to exclude:', np.unique(prof_list.values))

Restricting profiles
Max and min area between downcasts =  2.3917896006277424 0.00023975386366415294
List of profiles to exclude: [ 59.  60.  81.  82.  83.  84.  87.  88.  89.  90.  93.  94.  95.  96.
  99. 100. 107. 108. 117. 118. 125. 126. 169. 170. 183. 184. 185. 186.
 187. 188. 189. 190. 193. 194. 199. 200. 235. 236. 241. 242. 245. 246.
 255. 256. 262. 263. 268. 269. 270. 271. 272. 273. 276. 277. 278. 279.
 282. 283. 318. 319. 322. 323. 332. 333. 354. 355. 356. 357. 358. 359.
 364. 365. 370. 371. 372. 373. 376. 377. 394. 395. 418. 419. 430. 431.
 448. 449. 488. 489. 498. 499. 532. 533. 534. 535. 750. 751.]

ts_sub['profiles_to_exclude'] = ts_sub.profile_index.isin(prof_list.values)

# Plot the profiles that were kept for the comparison!!
print('Red indicates profile pairs that were identified in this process, where the area between \nprofile pairs was considered to be too large, and so are not included in the thermal lag correction. \nWhite bands indicate conductivity profiles with QC4 values')

subdata = ts_sub.where(ts_sub.profiles_to_exclude == True) #profile_index.isin(prof_list.values)#(profile_bins)

fig, ax = plt.subplots(1,1, figsize=(9, 3), constrained_layout=True)

ax.scatter(ts_sub.profile_index, ts_sub.pressure, marker = '.', color='black', s = 2,
           rasterized=True, label='Profiles with suspicious salinity and conductivity removed')


ax.set_ylim([MAX_DEPTH, 0])
ax.set_xlim([0,np.nanmax(ts_sub.profile_index)])
ax.set_xlabel('Profile index')
ax.set_ylabel('Depth (m)') 

ax.scatter(subdata.profile_index, subdata.pressure,
           color='red', marker = '.', s=2, rasterized=True, label = 'Profiles with large SD')

ax.legend(fontsize='small', loc='lower left');

Red indicates profile pairs that were identified in this process, where the area between 
profile pairs was considered to be too large, and so are not included in the thermal lag correction. 
White bands indicate conductivity profiles with QC4 values

Figure

# SAVING INTERMEDIATE FILE TO NETCDF
ts_sub.to_netcdf(f'{filepath}/{deploy_name}_goodprofiles.nc')
# Save a gridded version as well
outfile = make_gridfiles(f'{filepath}/{deploy_name}_goodprofiles.nc', f'{filepath}', deployfile, fnamesuffix='goodprofiles')

2.4.2 Defining the range to calculate $\alpha$ and $\tau$:

From examining the asymmetry in up-down profiles, we manually choose a range to apply the thermal lag correction to. We pick a 20 profile range where data without large standard deviations between up and downcasts.

# Select a subset of profiles to calculate tau and alpha
profile_lims = [600,620]
print(f'Using profile limits {profile_lims} for tau and alpha calculation')

Using profile limits [600, 620] for tau and alpha calculation

fname = f'{filepath}/{deploy_name}_goodprofiles.nc'
gridfname= f'{filepath}/{deploy_name}_gridgoodprofiles.nc'

ds=xr.open_dataset(f'{filepath}/{deploy_name}_gridgoodprofiles.nc') 

tbins = ds.temperature.mean(dim='time', skipna=True)
tbins = np.sort(tbins[::6])
tbins = tbins[np.isfinite(tbins)]
# print(tbins)
depbins = ds.depth[::6]

# Switches back to the time series....
# switching ds to ts
with xr.load_dataset(fname) as ds0:
    # USING PROFILE_LIMS DECIDED ABOVE!
    # print(f'Loading {fname}')
    inds=np.arange(profile_lims[0], profile_lims[1])
    
    indbins = np.arange(inds[0]-0.5, inds[-1]+0.5, 1.0)
    
    ts = ds0.where((ds0.profile_index >= inds[0]) & (ds0.profile_index <= inds[-1]), drop=False)
    # Once again, make sure to use density cutoff
    ts = ts.where((ts.density > density_cutoff), drop=True)
    # Also, profiles to exclude:
    ts = ts.where(ts.profiles_to_exclude == False, drop=True)
    
    sal = ts.salinity_adjusted

## find corrected salinity using smoothed conductivity 
def correct_sal(ds, fn, alpha, tau):
    a = 4 * fn * alpha * tau / (1 + 4*fn*tau)
    b = 1 - 2 * a / alpha
    aa = [1, b]
    bb = [a, -a]
    tempcorr = ds.temperature_adjusted.values.copy()
    tempcell = ds.temperature_adjusted.values.copy()
    good = ~np.isnan(tempcell)
    tempcorr[good] = signal.lfilter(bb, aa, ds.temperature_adjusted.values[good])
    tempcell = tempcell - tempcorr
    # tempcell = ds.temperature
    sal = sw.salt(10 * ds.conductivity_adjusted / sw.constants.c3515, tempcell, ds.pressure)
    dens = sw.dens(sal, tempcell, ds.pressure)
    return sal, tempcell, dens

Finding alpha and tau values with lowerest error estimates:

We will see if the constants found using dfo-bb046-20201006 improve this data.

dt = ts.time.diff(dim='time').mean(dim='time').astype('float') / 1e9
fn = 1.0 / dt
alpha = 0.14
tau = 8.96 

print('******')
print(f'Applying alpha = {alpha} and tau = {tau}')
print('*****')

alpha2 = 000.1
tau2 =0

########## 

with xr.load_dataset(fname) as ts0:
    inds = np.arange(0, NUM_PROFILES)
    
    indbins = np.arange(inds[0]-0.5, inds[-1]+0.5, 1.0)
    
    ts = ts0.where((ts0.profile_index >= inds[0]) & (ts0.profile_index <= inds[-1]), drop=False)
    ts = ts.where((ts.density > density_cutoff), drop=True)
    # Also, profiles to exclude:
    ts = ts.where(ts.profiles_to_exclude == False, drop=True)
    
    sal = correct_sal(ts, fn, alpha, tau)
    
    ss0, err0, totalerr = pgs.get_error(ts, ts.salinity_adjusted, tbins, indbins)
    ss, err, totalerr = pgs.get_error(ts, sal[0], tbins, indbins)

    sal2 = correct_sal(ts, fn, alpha2, tau2)
    ss2, err2, totalerr2 = pgs.get_error(ts, sal2[0], tbins, indbins)


Y_LIMS = [300,0]

fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, layout='constrained', figsize = [8,3])
sp0 = np.nanmean(ss0, axis=1)
pc = ax[0].pcolormesh(indbins[:-1], depbins[:len(tbins)-1][::-1], ss0-sp0[:, np.newaxis], cmap='RdBu_r', vmin=-0.1, vmax=0.1) #, vmin=3.34, vmax=3.52)
ax[1].set_ylim([200, 0])
ax[0].set_ylabel('ISOTHERM DEPTH [m]')
ax[0].set_xlabel('PROFILE')
ax[0].set_title('No correction')
ax[0].set_ylim(Y_LIMS)

pc = ax[1].pcolormesh(indbins[:-1], depbins[:len(tbins)-1][::-1], ss-sp0[:, np.newaxis], cmap='RdBu_r', vmin=-0.1, vmax=0.1) #, vmin=3.34, vmax=3.52)
ax[1].set_title(f'$\\alpha = {alpha}, \\tau = {tau}$')
fig.colorbar(pc, ax=ax, label='S[T(z), bin] - mean(S)[T(z)]')

fig, ax = plt.subplots(1)
ax.plot(np.nanmean(err, axis=1), depbins[:len(tbins)-1][::-1], label='Corr')
ax.plot(np.nanmean(err2, axis=1), depbins[:len(tbins)-1][::-1], label='Orig')
ax.grid('x')
ax.set_ylim([200, 0])
ax.legend(fontsize='small')
ax.set_xlabel('mean(err)')
ax.set_title('Avg. mismatch between upcast and downcast salinity')
ax.set_ylabel('Isotherm depth (m)')

******
Applying alpha = 0.14 and tau = 8.96
*****

Text(0, 0.5, 'Isotherm depth (m)')

Figure

Figure

Zoom into three areas to observe change

def zoom_in_thermal_chg(min_ind,max_ind,fname,density_cutoff):
    with xr.load_dataset(fname) as ts0:
        inds = np.arange(min_ind, max_ind)
        
        indbins = np.arange(inds[0]-0.5, inds[-1]+0.5, 1.0)
        ts = ts0.where((ts0.profile_index >= inds[0]) & (ts0.profile_index <= inds[-1]), drop=False)
        # Once again, make sure to use density cutoff? maybe?
        ts = ts.where((ts.density > density_cutoff), drop=True)
        # Also, profiles to exclude:
        ts = ts.where(ts.profiles_to_exclude == False, drop=True)
        
        sal = correct_sal(ts, fn, alpha, tau )

        ss0, err0, totalerr = pgs.get_error(ts, ts.salinity_adjusted, tbins, indbins)
        ss, err, totalerr = pgs.get_error(ts, sal[0], tbins, indbins)
    
        sal2 = correct_sal(ts, fn, alpha2, tau2)
        ss2, err2, totalerr2 = pgs.get_error(ts, sal2[0], tbins, indbins)
    
    
    Y_LIMS = [200,0]
    
    fig, ax = plt.subplots(1, 3, layout='constrained', figsize = [12,3])
    sp0 = np.nanmean(ss0, axis=1)
    pc = ax[0].pcolormesh(indbins[:-1], depbins[:len(tbins)-1][::-1], ss0-sp0[:, np.newaxis], cmap='RdBu_r', vmin=-0.1, vmax=0.1) #, vmin=3.34, vmax=3.52)
    ax[1].set_ylim([200, 0])
    ax[0].set_ylabel('ISOTHERM DEPTH [m]')
    ax[0].set_xlabel('PROFILE')
    ax[0].set_title('No correction')
    ax[0].set_ylim(Y_LIMS)
    
    pc = ax[1].pcolormesh(indbins[:-1], depbins[:len(tbins)-1][::-1], ss-sp0[:, np.newaxis], cmap='RdBu_r', vmin=-0.1, vmax=0.1) #, vmin=3.34, vmax=3.52)
    ax[1].set_title(f'$\\alpha = {alpha}, \\tau = {tau}$')
    fig.colorbar(pc, ax=ax, label='S[T(z), bin] - mean(S)[T(z)]')

    ax[2].set_title('Avg. mismatch between upcast and downcast salinity')
    ax[2].plot(np.nanmean(err, axis=1), depbins[:len(tbins)-1][::-1], label='Corr')
    ax[2].plot(np.nanmean(err2, axis=1), depbins[:len(tbins)-1][::-1], label='Orig')
    ax[2].grid('x')
    ax[2].set_ylim([200, 0])
    ax[2].legend(fontsize='small')
    ax[2].set_xlabel('mean(err)')
    ax[2].axvline(x=0,c='k',ls='--')

zoom_in_thermal_chg(200,260,fname,density_cutoff)
zoom_in_thermal_chg(500,560,fname,density_cutoff) 
zoom_in_thermal_chg(700,760,fname,density_cutoff)

Figure

Figure

Figure

The applied constants reduced the up- and down-cast asymmetry at the beginning, middle, and end of the time series. Therefore, the thermal lag correction is applied to the dataset.

Save corrected data :

These fields, adjusted and re-calculated using the new alpha and tau values, are saved in the output file as fields salinity_adjusted, temperature_adjusted and density_adjusted. The original delayed-mode, uncorrected fields are saved as salinity, temperature and density.

print('*****')
print(f'Saving with alpha = {alpha}, tau = {tau}')
print('*****')

# We're going to apply the changes to the data with quality control flags data. That step happened last here:

ts = xr.open_dataset(f'{filepath}/{deploy_name}_QC4.nc')
ds = xr.open_dataset(f'{filepath}/{deploy_name}_gridQC4.nc')

##apply thermal lag correction 
s, t, d =correct_sal(ts, fn, alpha, tau)

prof_range= [] #we didn't smooth in this mission 
######### Sub in the smoothed data from the clogged section into the new adjusted data 
adjust = xr.where(
    ts.profile_index.isin(prof_range),
    ts['salinity_adjusted'],s  # keep where True
                            # replace where False
)

*****
Saving with alpha = 0.14, tau = 8.96
*****

########
ts.attrs['processing_details'] = 'Processing details are located on the C-PROOF website for this mission under the reports tab.'  
ts.attrs['processing_tech'] = 'Lauryn Talbot; ltalbot@uvic.ca'  
ts.attrs['citation'] = '"Klymak, J., & Ross, T. (2025). C-PROOF Underwater Glider Deployment Datasets [Data set]. Canadian-Pacific Robotic Ocean Observing Facility.doi:10.82534/44DS-K310"'
ts.attrs['references'] = 'https://doi.org/10.82534/44DS-K310' 
ts.attrs['correction_constants'] = f'alpha = {alpha}; tau={tau}'
ts.attrs['quality_flags']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

# Uncorrected conductivity
ts['conductivity'].attrs['comment'] = 'uncorrected conductivity'
ts['conductivity_QC'] = ts.conductivity_QC
ts['conductivity_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']


# Adjusted (aka cleaned) conductivity

ts['conductivity_adjusted'] = ts.conductivity_adjusted
ts['conductivity_adjusted'].attrs['comment'] = 'adjusted conductivity'
ts['conductivity_adjusted'].attrs['processing_report'] = processing_report
ts['conductivity_adjusted'].attrs['processing_date'] = processing_date
ts['conductivity_adjusted'].attrs['processing_date'] = processing_protocol
ts['conductivity_adjusted_QC'] = ts['conductivity_adjusted_QC']   
ts['conductivity_adjusted_QC'].attrs['comment']=['QC 8 data smoothed due to clogged cell; resolution is approx. 9 km; 1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']


#remove conductivityClean - an intermediate step 
ts = ts.drop_vars('conductivityClean')


# Uncorrected temperature
ts['temperature'].attrs['comment'] = 'uncorrected temperature [degC]'
ts['temperature_QC'] = ts.temperature_QC 
ts['temperature_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']


# Adjusted temperature
ts['temperature_adjusted'] = ts.temperature_adjusted.copy()
ts['temperature_adjusted'].attrs['comment'] = 'temperature [degC]'
ts['temperature_adjusted'].attrs['processing_report'] = processing_report
ts['temperature_adjusted'].attrs['processing_date'] = processing_date
ts['temperature_adjusted'].attrs['processing_date'] = processing_protocol
ts['temperature_adjusted_QC'] =  ts['temperature_adjusted_QC']  
ts['temperature_adjusted_QC'].attrs['comment']=['QC 8 data smoothed due to clogged cell; resolution is approx. 9 km; 1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

# Uncorrected salinity
ts['salinity'].attrs['comment'] = 'uncorrected salinity [psu]'

# Corrected salinity
ts['salinity_adjusted'] = adjust
ts['salinity_adjusted'].attrs['comment'] = f'adjusted salinity [psu] using a thermal lag correction with alpha = {alpha} and tau = {tau} '
ts['salinity_adjusted'].attrs['method'] = ' '
ts['salinity_adjusted'].attrs['processing_report'] = processing_report
ts['salinity_adjusted'].attrs['processing_date'] = processing_date
ts['salinity_adjusted'].attrs['processing_protocol'] = processing_protocol
ts['salinity_adjusted_QC'] =  ts['salinity_adjusted_QC'] 
ts['salinity_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']
ts['salinity_adjusted_QC'].attrs['comment']=['QC 8 data smoothed due to clogged cell; resolution is approx. 9 km; 1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

# Unadjusted density
ts['density'].attrs['comment'] = 'unadjusted density'

# Adjusted density
ts['density_adjusted'] = ('time', d)
ts['density_adjusted'].attrs['comment'] = 'density from adjusted salinity [psu] and temperature [degC]'
ts['density_adjusted'].attrs['method'] = ' '
ts['density_adjusted_QC'] =  ts['salinity_adjusted_QC'] #density is only as good as the salinity values used to determine it  
ts['density_QC'] = ts['salinity_QC'] 
ts['density_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']
ts['density_adjusted_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

Calculate potential_temperature_adjusted and potential_density_adjusted using adjusted_salinity

long = ts.longitude.fillna(ts.longitude.mean(skipna=True))
lat = ts.latitude.fillna(ts.latitude.mean(skipna=True))
sa_adj = gsw.SA_from_SP(ts['salinity_adjusted'],ts['pressure'],long,lat)
ct_adj = gsw.CT_from_t(sa_adj,ts['temperature_adjusted'],ts['pressure'])
ts['potential_density_adjusted'] = (('time'),1000 + gsw.density.sigma0(sa_adj,ct_adj).values)
ts['potential_density_adjusted'].attrs['comment'] = 'calculated using adjusted salinity'
ts['potential_density_adjusted_QC'] = ts['salinity_adjusted_QC']

ts['potential_temperature_adjusted'] = (('time'),
                                        gsw.conversions.pt0_from_t(ts.salinity_adjusted,ts.temperature_adjusted,ts.pressure).values)
ts['potential_temperature_adjusted'].attrs['comment'] = 'calculated using adjusted salinity'
ts['potential_temperature_adjusted_QC'] = ts['salinity_adjusted_QC']

##### Save our final datasets
ts.to_netcdf(f'{filepath}{deploy_name}_CTDadjusted.nc')
print(f'Corrected data saved to file: {filepath}/{glider_name}/{deploy_name}/{deploy_name}_CTDadjusted.nc')
make_gridfiles(f'{filepath}{deploy_name}_CTDadjusted.nc', 
                         f'{filepath}', deployfile, fnamesuffix='_CTDadjusted')

Corrected data saved to file: deployments/dfo-bb046/dfo-bb046-20200908//dfo-bb046/dfo-bb046-20200908/dfo-bb046-20200908_CTDadjusted.nc

'deployments/dfo-bb046/dfo-bb046-20200908//dfo-bb046-20200908_grid_CTDadjusted.nc'

3.0 Summary of corrections applied to delayed mode data for this mission

Identification of anomalous conductivity values:

• Anomalous conductivity values, values that still differ from the mean by more than 3 standard deviations, are flagged as 'bad' (QC 4).
• Profiles with a clogged sensor were identified as QC 3 and smoothed over, making the resolution 9 km.

Identification of questionable temperature profiles:

• Anomolous temperature profiles were flagged as Q4 if a)individual points where the temperature differs from the average of adjacent samples by more than 0.75 °C, b)differences in temperature between adjacent profiles at the same depth that exceed 1.5 °C, and c) occurrences where the density of the deeper depth bin is lower than that of the shallower bin, indicating physically unrealistic structure. Density inversions are only considered below 100 m.
• Profiles with a clogged sensor were identified as QC 3 and smoothed over, making the resolution 9 km.

Identification of questionable salinity profiles:

• The worst QC flag assigned to either conductivity or temperature is also applied to salinity.

Sensor alignment correction:

• No sensor alignment correction was applied.

Thermal lag correction:

• The directly determined values for the thermal lag correction produced an improvement that was larger than the recommended values from Janzen and Creed (2011).
• The correction overall significantly reduced the root-mean squared difference for the area between between pairs of profiles.
• The final thermal lag correction was applied using the calculated values of:

print(f'alpha = {alpha} and tau = {tau}')

alpha = 0.14 and tau = 8.96

ds=xr.open_dataset(f'{filepath}{deploy_name}_grid_CTDadjusted.nc')
ts=xr.open_dataset(f'{filepath}{deploy_name}_CTDadjusted.nc')

#Compare the uncorrected and corrected data in T-S space

print('Temperature-salinity diagrams for all profiles, '
      'showing the difference between upcasts (red) and downcasts (blue), '
      'for the data without the thermal lag correction applied (left panel) and '
      'the data with the thermal lag correction applied (right panel):')

x_lim=[30, 34.5]

#Plotting
tsno4 = ts.where((ts.salinity_adjusted_QC !=4))
tsno4 = tsno4.where((ts.temperature_adjusted_QC !=4))

fig, ax = plt.subplots(1, 2, sharex=True, sharey=True, figsize=(9,5))

ind = np.where(tsno4.profile_direction.values== 1)[0]
ax[0].plot(tsno4.salinity[ind], tsno4.temperature[ind], 'b.', markersize=2, rasterized=True, label = 'Downcast')
ax[1].plot(tsno4.salinity_adjusted[ind], tsno4.temperature_adjusted[ind], 'b.', markersize=2, rasterized=True, label = 'Downcast')

ind = np.where(tsno4.profile_direction.values == -1)[0]
ax[0].plot(tsno4.salinity[ind], tsno4.temperature[ind], 'r.', markersize=2, alpha = 0.5, rasterized=True, label = 'Upcast')
ax[1].plot(tsno4.salinity_adjusted[ind], tsno4.temperature_adjusted[ind], 'r.', markersize=2, alpha = 0.5, rasterized=True, label = 'Upcast')

S_range = np.linspace(int(np.min(tsno4.salinity)-0.5), 
                      int(np.max(tsno4.salinity)+0.5), 1000)
T_range = np.linspace(int(np.min(tsno4.temperature)-1), 
                      int(np.max(tsno4.temperature)+1), 1000)
S_grid, T_grid = np.meshgrid(S_range, T_range)
density_grid = seawater.eos80.dens0(S_grid, T_grid)


CS = ax[0].contour(S_range, T_range, density_grid,
                np.arange(1021,np.round(np.max(density_grid)),0.5),
                colors='k', linewidths=0.5);
ax[0].clabel(CS, CS.levels, inline=True, fontsize=10)
ax[0].set_ylabel('Temperature [$^o$C]')
ax[0].set_xlabel('Salinity [psu]')
ax[0].set_title('Before correction')
ax[0].set_xlim(x_lim)
ax[0].grid()

CS = ax[1].contour(S_range, T_range, density_grid,
                np.arange(1021,np.round(np.max(density_grid)),0.5),
                colors='k', linewidths=0.5);
ax[1].clabel(CS, CS.levels, inline=True, fontsize=10)

ax[1].set_ylabel('Temperature [$^o$C]')
ax[1].set_xlabel('Salinity [psu]')
ax[1].set_title(f'(After correction: tau = {tau}, alpha = {alpha})' )
ax[1].grid()

ax[0].legend(prop={'size': 10});
ax[1].legend(prop={'size': 10});

Temperature-salinity diagrams for all profiles, showing the difference between upcasts (red) and downcasts (blue), for the data without the thermal lag correction applied (left panel) and the data with the thermal lag correction applied (right panel):

Figure

# Visualize the final data
fig, ax = plt.subplots(1, 1, figsize=(6, 6))
X_LIM = [28,34.5]
# T-S diagram for fully corrected data
ax0 = ax
ax0.plot(ts.salinity,ts.temperature,'k.',markersize=2, label = "Delayed-mode data")


tsno4 = ts.where((ts.salinity_adjusted_QC !=4))
tsno4 = tsno4.where((ts.temperature_adjusted_QC !=4))

ax0.plot(tsno4.salinity_adjusted,tsno4.temperature_adjusted,'.',markersize=2, label = "Adjusted and filtered data without Q4 data")

#Create a density grid to contour plot isopycnals
S_range = np.linspace(np.nanmin(ts.salinity_adjusted)-0.5, 
                      np.nanmax(ts.salinity_adjusted)+0.5, 1000)
T_range = np.linspace(np.nanmin(ts.temperature_adjusted)-1, 
                      np.nanmax(ts.temperature_adjusted)+1, 1000)
S_grid, T_grid = np.meshgrid(S_range, T_range)
density_grid = seawater.eos80.dens0(S_grid, T_grid)

CS = ax0.contour(S_range, T_range, density_grid,
                np.arange(1014,
                          np.round(np.max(density_grid)),0.5),
                colors='k', linewidths=0.5);
ax0.clabel(CS, CS.levels, inline=True, fontsize=10)
ax0.set_xlabel('Salinity [psu]', fontsize=18)
ax0.set_ylabel('Temperature [$^o$C]', fontsize=18)
ax0.set_xlim(X_LIM)
ax0.grid()
ax0.legend()

print('The corrected temperature and salinity fields '
      'shown in a T-S diagram with density contours:')

The corrected temperature and salinity fields shown in a T-S diagram with density contours:

Figure

fig, axs = plt.subplots(4, 1, figsize=(11, 10), sharey=True, sharex=True)

xlims = [0, NUM_PROFILES]
ylims=[MAX_DEPTH,0]
# ylims=[50,0]


dsno4 = ds.where((ds.salinity_adjusted_QC !=4))

pc = axs[0 ].pcolormesh(dsno4.profile, dsno4.depth, dsno4['salinity_adjusted'],rasterized=True)
axs[0].set_ylim(ylims)
axs[1].set_xlim(xlims)
fig.colorbar(pc, ax=axs[0], label = 'Salinity [psu]')
axs[0].set_title('Salinity',loc='left')

dsno4 = ds.where((ds.temperature_adjusted_QC !=4))

pc = axs[1].pcolormesh(dsno4.profile, dsno4.depth, dsno4['temperature_adjusted'],rasterized=True,cmap='plasma')

fig.colorbar(pc, ax=axs[1], label = 'Temperature [$^o$C]')
axs[1].set_title('Temperature',loc='left')

dsno4 = ds.where((ds.conductivity_adjusted_QC !=4))
pc = axs[2].pcolormesh(dsno4.profile, dsno4.depth, dsno4['conductivity_adjusted'],rasterized=True,cmap='cividis')
fig.colorbar(pc, ax=axs[2], label = 'Conductivity [S/m]')
axs[2].set_title('Conductivity',loc='left')

# pc = axs[3].pcolormesh(ds.profile, ds.depth, ds['oxygen_concentration'],rasterized=True,cmap='inferno')
# fig.colorbar(pc, ax=axs[3])
# axs[3].set_title('Oxygen Concentration',loc='left')

pc = axs[3].pcolormesh(dsno4.profile, dsno4.depth, dsno4['density'],rasterized=True,cmap='inferno')
fig.colorbar(pc, ax=axs[3], label = 'Density [kg/m$^3$]')
axs[3].set_title('Density',loc='left')

axs[0].set_ylabel('Depth [m]');
axs[1].set_ylabel('Depth [m]');
axs[2].set_ylabel('Depth [m]');
axs[3].set_ylabel('Depth [m]');

Figure

display(Markdown('./docs/CTD_References.md'))

References

1. Ferrari, R., and Rudnick, D. L. Thermohaline variability in the upper ocean, J. Geophys. Res., 105(C7), 16857-16883, 2000.

2. Garau, B., Ruiz, S., Zhang, W. G., Pascual, A., Heslop, E., Kerfoot, J., & Tintoré, J. Thermal Lag Correction on Slocum CTD Glider Data, J. Atmos. Oceanic Technol., 28(9), 1065-1071, 2011.

3. Janzen, C. D., and Creed, E. L. Physical oceanographic data from Seaglider trials in stratified coastal waters using a new pumped payload CTD, OCEANS'11 MTS/IEEE KONA, Waikoloa, HI, USA, 1-7, 2011.

4. Morison, J., Andersen, R., Larson, N., D’Asaro, E., & Boyd, T. The correction for thermal-lag effects in Sea-Bird CTD data, J. Atmos. Oceanic Technol., 11, 1151-1164, 1994.

5. Sea-Bird Seasoft V2:SBE Data Processing - CTD Data Processing & Plotting Software for Windows, Sea-Bird Scientific, software manual revision 7.26.8, 2017.

6. Sea-Bird User Manual - GPCTD Glider Payload CTD (optional DO) - Conductivity, Temperature, and Pressure (optional DO) Sensor with RS-232 Interface, Sea-Bird Scientific, manual version 008, 2021.