glider_name = 'dfo-bb046'
deploy_name = f'{glider_name}-20210413'
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-20210413
* Sensor: CTD_0278

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

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

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

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-20210413//L0-timeseries/dfo-bb046-20210413_delayed.nc
************
* There are 474809 data points in total, with 1031.0 profiles
* Time period: 2021-04-13 to 2021-05-05
* Depth range: 0 - 703 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:

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

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

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

t4 = ts.where(ts.conductivity_QC ==4)
ax[0].plot(ts.conductivity, ts.temperature, color='k', marker='.', linestyle='none', label='QC 1')
ax[0].plot(t4.conductivity, ts.temperature, color='r', marker='.', linestyle='none', label='QC 4')
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='k', marker='.', linestyle='none')
ax[1].plot(t4.profile_index, t4.conductivity, color='r', 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), '
      '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), and conductivity vs. profile index (right), with the red dots showing the unphysical values flagged as bad and are flagged as QC 4:

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

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


#put in isopycnal space 
pot_density = gsw.pot_rho_t_exact(ds.salinity, ds.potential_temperature, ds.pressure,0)
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) 


iso_cond = np.empty(np.shape(ds.conductivity))

interpolated_cond=[]
for i in range(len(ds.time)): 
    x= np.array(iso_density.isel(time=i))
    xp = np.array(pot_density.isel(time=i))
    yp = np.array(ds.conductivity.isel(time=i))
    #interpolate 
    interpolated_cond.append(np.interp(x,xp,yp, left=np.nan, right=np.nan))

ds['iso_cond']=((('depth','time')),np.transpose(interpolated_cond))
mn_cond = ds['iso_cond'].mean('time')
std_cond = ds['iso_cond'].std('time')

ds['cond_anomalies'] = (ds.iso_cond  -mn_cond )/std_cond

import matplotlib.gridspec as gridspec

fig = plt.figure(figsize=(10, 4))

# Grid with 4 columns: widths 1, 1, 2, 2.
gs = gridspec.GridSpec(1, 4, width_ratios=[1, 1, 2, 2])

ax0 = fig.add_subplot(gs[0])
ax1 = fig.add_subplot(gs[1], sharey=ax0)    
ax2 = fig.add_subplot(gs[2:], sharey=ax0)  

# --- Left panel: mean cond ---
ax0.plot(mn_cond, -ds.depth)
ax0.set_ylabel('Mean isopycnal depth (m)')
ax0.set_title('mean C')

# --- Middle panel: std cond ---
ax1.plot(std_cond, -ds.depth)
ax1.set_title('std C')

# --- Right double-width panel: pcolor ---
im = ax2.pcolor(ds.profile_index, -ds.depth, ds.cond_anomalies,
                cmap='RdBu_r', vmin=-1, vmax=1)

ax2.set_title(r'$ \frac{C_{\text{iso}} - \bar{C}_{\text{iso}}}{\text{std}(C_{\text{iso}})}$')
ax2.set_xlabel('Profile index')
ax2.set_ylim(-300, 0)

# Colorbar
fig.colorbar(im, ax=ax2)
plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax2.get_yticklabels(), visible=False)

[None, None, None, None, None, None, None]

# Determine cond values that are not temporarily flagged as bad 
dC = np.fabs(ds.cond_anomalies - ds.cond_anomalies.mean(dim='time')) 
CNoflag = ds.cond_anomalies.where((dC < flag_stdev * ds.cond_anomalies.std(dim='time'))) 
# Determine cond values that are 'clean' - not flagged as bad 
dC = np.fabs(ds.cond_anomalies - CNoflag.mean(dim='time')) 
CClean = ds.cond_anomalies.where((dC < clean_stdev * CNoflag.std(dim='time'))) 

# Find the sum of clean values in each profile relative to the total number of values 
NC = np.isfinite(ds.cond_anomalies).sum(dim='depth')
NCClean = np.isfinite(CClean).sum(dim='depth')

# Fraction bad (avoid divide-by-zero warnings)
bad_frac = xr.where(NC > 0, (NC - NCClean) / NC, 1.0)
# NC == 0  → set bad_frac = 1 (100% bad)

clean_cutoff = 0.1
bad_profiles = np.unique(
    ds.profile_index.where(bad_frac >= clean_cutoff, drop=True)
)

print(f'Bad profiles: {bad_profiles}')
bad_profiles = xr.DataArray(bad_profiles)

Bad profiles: [ 247.  249.  397.  398.  399.  400.  401.  402.  403.  404.  405.  406.
 1030. 1031.   nan]

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.5
vmax=1.5


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

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

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['cond_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
fig.colorbar(pc, ax=ax, label = 'T anomalies')
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')
fig.colorbar(pc, ax=ax, label = 'C anomalies')
########
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['cond_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')
fig.colorbar(pc, ax=ax, label = 'C anomalies')

######
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['cond_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')
fig.colorbar(pc, ax=ax, label = 'C anomalies')

#########
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['cond_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')
fig.colorbar(pc, ax=ax, label = 'C anomalies')
 
print('Zooming in along the glider deployment to visualize temperature spikes.')

Zooming in along the glider deployment to visualize temperature spikes.

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

prof_range = np.r_[bad_profiles,246,24,346,1028 ]
#########
ds['conductivity_QC'] = xr.where(ds.profile_index.isin(prof_range), 3,  ds['conductivity_QC'])

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

##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)
    & (np.abs(temp_fwd - temp_bwd) < 0.5)
)

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

ts.to_netcdf(f'{filepath}/{deploy_name}_QC2.nc')
outfile = make_gridfiles(f'{filepath}/{deploy_name}_QC2.nc', f'{filepath}', deployfile, fnamesuffix='QC2')

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

#put in isopycnal space 
pot_density = gsw.pot_rho_t_exact(ds.salinity, ds.potential_temperature, ds.pressure,0)
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) 


iso_temp = np.empty(np.shape(ds.temperature))

interpolated_temp=[]
for i in range(len(ds.time)): 
    x= np.array(iso_density.isel(time=i))
    xp = np.array(pot_density.isel(time=i))
    yp = np.array(ds.temperature.isel(time=i))
    #interpolate 
    interpolated_temp.append(np.interp(x,xp,yp, left=np.nan, right=np.nan))

ds['iso_temp']=((('depth','time')),np.transpose(interpolated_temp))
mn_temp = ds['iso_temp'].mean('time')
std_temp = ds['iso_temp'].std('time')

ds['temp_anomalies'] = (ds.iso_temp  -mn_temp )/std_temp

import matplotlib.gridspec as gridspec

fig = plt.figure(figsize=(10, 4))

# Grid with 4 columns: widths 1, 1, 2, 2.
gs = gridspec.GridSpec(1, 4, width_ratios=[1, 1, 2, 2])

ax0 = fig.add_subplot(gs[0])
ax1 = fig.add_subplot(gs[1], sharey=ax0)    
ax2 = fig.add_subplot(gs[2:], sharey=ax0)  

# --- Left panel: mean cond ---
ax0.plot(mn_temp, -ds.depth)
ax0.set_ylabel('Mean isopycnal depth (m)')
ax0.set_title('mean temp')

# --- Middle panel: std cond ---
ax1.plot(std_temp, -ds.depth)
ax1.set_title('std temp')

# --- Right double-width panel: pcolor ---
im = ax2.pcolor(ds.profile_index, -ds.depth, ds.temp_anomalies,
                cmap='RdBu_r', vmin=-1, vmax=1)

ax2.set_title(r'$ \frac{T_{\text{iso}} - \bar{T}_{\text{iso}}}{\text{std}(T_{\text{iso}})}$')
ax2.set_xlabel('Profile index')
ax2.set_ylim(-300, 0)

# Colorbar
fig.colorbar(im, ax=ax2)
plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax2.get_yticklabels(), visible=False)

[None, None, None, None, None, None, None]

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.5
vmax=1.5


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

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

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['temp_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
fig.colorbar(pc, ax=ax, label = 'T anomalies')
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')


########
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['temp_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'T anomalies')
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')

######
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['temp_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)

fig.colorbar(pc, ax=ax, label = 'T anomalies')
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')

#########
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['temp_anomalies'],rasterized=True,vmin=vmin,vmax=vmax,cmap=cmap)
x = bad_profiles.where(bad_profiles.isin(ds_sub.profile),drop=True)
y = 0 + np.zeros_like(x)
ax.scatter(x,y,60,marker='v',color='green')

fig.colorbar(pc, ax=ax, label = 'T anomalies')
 
print('Zooming in along the glider deployment to visualize temperature spikes.')

Zooming in along the glider deployment to visualize temperature spikes.

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

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

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(
    (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

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

s4 = ts.where((ts.salinity_QC == 4) | (ts.salinity_QC == 3) ) 

ax[0].plot(ts.salinity, ts.temperature, color='k', marker='.', linestyle='none', label='QC 1')
ax[0].plot(s4.salinity, s4.temperature, color='r', marker='.', linestyle='none', label='QC 3 & 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()

t4 = ts.where((ts.salinity_QC == 4) | (ts.salinity_QC == 3) ) 

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

ax[1].set_ylabel('Salinity ', 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 3 & 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 3 & QC 4 shown in red and indicated by the red arrows at the top of the panel on the right:

fig,axs=plt.subplots()

vmin=30
vmax = 35
ds4 = ds.where(ds.salinity_QC==1)


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

<matplotlib.colorbar.Colorbar at 0x37fdc1bd0>

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

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

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)-2000, int(N/2)+2000)
    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, c='red', label='1.25xConductivity')
    ax.set_ylabel('Conductivity/Temperature',)
    ax.set_xlabel('Time',)

    ax.plot(ts1.time,ts1.temperature ,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.

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

# 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==1) 
ts1 = ts.where(ts.salinity_QC==1)

Total number of profiles: 1031

fig, ax = plt.subplots()
vmin=31
vmax=34
pc=ax.pcolormesh(ds1.profile, -ds1.depth, ds1['salinity'],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

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

#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 =  0.5752466382037338 7.030160270280424e-05
List of profiles to exclude: [ 134.  135.  156.  157.  208.  209.  230.  231.  232.  233.  248.  250.
  251.  282.  283.  284.  285.  286.  287.  294.  295.  352.  353.  354.
  355.  362.  363.  368.  369.  374.  375.  376.  377.  378.  379.  380.
  381.  390.  391.  420.  421.  434.  435.  436.  437.  440.  441.  477.
  478.  479.  480.  481.  482.  487.  488.  489.  490.  491.  492.  495.
  496.  497.  498.  499.  500.  527.  528.  551.  552.  553.  554.  555.
  556.  557.  558.  563.  564.  565.  566.  567.  568.  569.  570.  571.
  572.  573.  574.  639.  640.  641.  642.  659.  660.  667.  668.  669.
  670.  679.  680.  711.  712.  733.  734.  739.  740.  757.  758.  759.
  760.  807.  808.  831.  832.  845.  846.  968.  969. 1014. 1015. 1020.
 1021. 1022. 1023. 1024. 1025.]

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 QC3 or 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.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 QC3 or QC4 values

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

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

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

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 = pgs.correct_sal(ts, fn, alpha, tau)
    
    ss0, err0, totalerr = pgs.get_error(ts, ts.salinity, tbins, indbins)
    ss, err, totalerr = pgs.get_error(ts, sal, tbins, indbins)

    sal2 = pgs.correct_sal(ts, fn, alpha2, tau2)
    ss2, err2, totalerr2 = pgs.get_error(ts, sal2, 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.2, vmax=0.2) #, 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.2, vmax=0.2) #, 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)')

******
Applying alpha = 0.48 and tau = 2.14
*****

Text(0.5, 0, 'mean(err)')

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 = pgs.correct_sal(ts, fn, alpha, tau )

        ss0, err0, totalerr = pgs.get_error(ts, ts.salinity, tbins, indbins)
        ss, err, totalerr = pgs.get_error(ts, sal, tbins, indbins)
    
        sal2 = pgs.correct_sal(ts, fn, alpha2, tau2)
        ss2, err2, totalerr2 = pgs.get_error(ts, sal2, tbins, indbins)
    
    
    Y_LIMS = [400,100]
    
    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(Y_LIMS)
    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].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(Y_LIMS)
    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(420,480,fname,density_cutoff) 
zoom_in_thermal_chg(600,660,fname,density_cutoff)

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

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


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



##apply thermal lag correction 
s =pgs.correct_sal(ts, fn, alpha, tau)

*****
Saving with alpha = 0.48, tau = 2.14
*****

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

#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.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_QC']  
ts['temperature_adjusted_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

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


# Corrected salinity
ts['salinity_adjusted'] = xr.where(ts.salinity_QC == 1, s, ts.salinity)
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_QC'] 
ts['salinity_adjusted_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

# Unadjusted density
ts['density'].attrs['comment'] = 'unadjusted density'
ts['density_QC'] =  ts['salinity_QC'] #same as salinity since derived 
ts['density_QC'].attrs['comment']=['1 = good data; 3 = bad data, potentially correctable; 4 = bad data; 8 = estimated data']

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-20210413//dfo-bb046/dfo-bb046-20210413/dfo-bb046-20210413_CTDadjusted.nc

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

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

alpha = 0.48 and tau = 2.14

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

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

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:

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 ==1))

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('Q1 Salinity adjusted',loc='left')

dsno4 = ds.where((ds.temperature_adjusted_QC ==1))
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('Q1 Temperature adjusted',loc='left')

dsno4 = ds.where((ds.conductivity_QC ==1))
pc = axs[2].pcolormesh(dsno4.profile, dsno4.depth, dsno4['conductivity'],rasterized=True,cmap='cividis')
fig.colorbar(pc, ax=axs[2], label = 'Conductivity [S/m]')
axs[2].set_title('Q1 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')

dsno4 = ds.where((ds.density_QC ==1))
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]');

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

CTD corrections applied to delayed-mode data¶

1.0 Preamble¶

1.1 Set up the processing¶

1.2 Profile Check¶

1.3 Delayed-mode data prior to corrections¶

2.0 Corrections applied to delayed mode data for this mission¶

Apply QC flags to data following Argo notation¶

2.1.1 Identification and removal of anomalous conductivity values¶

Look at anomalies relative to the mean along an isopycanl to observe any sections with clogged cells¶

2.2 Identifying questionable temperature profiles¶

Manually identify questionable profiles¶

Apply salinity flags¶

Identifying questionable salinity profiles¶

2.3 Sensor alignment correction¶

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

Finding alpha and tau values with lowerest error estimates:¶

Save corrected data :¶

Calculate potential_temperature_adjusted and potential_density_adjusted using adjusted_salinity¶

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

References¶