Glidertest demo

The purpose of this notebook is to demostrate the functionality of some glidertests functions. This notebook can also be used to diagnose issues within your glider data. We have added suggested processing in some cases.

The demo notebook is roughly organised as follows:

Load dataset
General mission overview
Sensors data overview

The sensor section is divided by data type or computed values

CTD
Dissolved oxygen
Bio-optics (e.g., WETLabs or Tridente)
Photosynthetically active radiation (PAR)
Vertical velocity (from a flight model)

[1]:

import matplotlib.pyplot as plt
import numpy as np
from glidertest import fetchers
from glidertest import tools, utilities, plots
from ioos_qc import qartod
import logging
import sys
import gsw
logging.basicConfig(stream=sys.stdout, level=logging.INFO)

2. General mission overview

Visualize the glider track, hydrography, sampling resolution, profile monotonicity, profile duration, and maximum depth

Glider track

[5]:

plots.plot_glider_track(ds)

/home/runner/micromamba/envs/TEST/lib/python3.13/site-packages/cartopy/io/__init__.py:263: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_land.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/runner/micromamba/envs/TEST/lib/python3.13/site-packages/cartopy/io/__init__.py:263: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_ocean.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)
/home/runner/micromamba/envs/TEST/lib/python3.13/site-packages/cartopy/io/__init__.py:263: DownloadWarning: Downloading: https://naturalearth.s3.amazonaws.com/10m_physical/ne_10m_coastline.zip
  warnings.warn(f'Downloading: {url}', DownloadWarning)

[5]:

(<Figure size 1400x1000 with 2 Axes>,
 <GeoAxes: title={'center': 'Glider Track'}, xlabel='Longitude', ylabel='Latitude'>)

Hydrographic overview

[6]:

plots.plot_basic_vars(ds,v_res=1, start_prof=0, end_prof=int(ds.PROFILE_NUMBER.max()))

[6]:

(<Figure size 1400x1000 with 5 Axes>,
 array([<Axes: xlabel='Temperature \n(°C)', ylabel='Depth (m)'>,
        <Axes: xlabel='Dissolved oxygen \n(mmol m$^{-}$³)'>], dtype=object))

[7]:

plots.plot_section(ds)

[7]:

(<Figure size 1400x1450 with 6 Axes>,
 None,
 <matplotlib.colorbar.Colorbar at 0x733c19465450>,
 <SecondaryAxis: >)

[8]:

plots.plot_ts(ds)

INFO:glidertest.utilities:Filtered values between 0.5% and 99.5% percentiles.

[8]:

(<Figure size 1400x1000 with 6 Axes>,
 array([<Axes: xlabel='Frequency', ylabel='Conservative Temperature (°C)'>,
        <Axes: title={'center': '2D Histogram \n (Log Scale)'}>,
        <Axes: ylabel='Log Counts'>, <Axes: >,
        <Axes: xlabel='Absolute Salinity', ylabel='Frequency'>, <Axes: >],
       dtype=object))

[9]:

plots.plot_max_depth_per_profile(ds)

[9]:

(<Figure size 1400x1000 with 2 Axes>,
 array([<Axes: title={'center': 'Max depth per profile'}, xlabel='Profile number', ylabel='Max depth (m)'>,
        <Axes: title={'center': 'Histogram of max depth per profile'}, xlabel='Max depth (m)', ylabel='Number of profiles'>],
       dtype=object))

Sampling resolution and spacing

[10]:

plots.plot_grid_spacing(ds)

INFO:glidertest.utilities:Depth and time differences have been filtered to the middle 99% of values.
INFO:glidertest.utilities:Numeric median/mean/max/min values are based on the original data.

[10]:

(<Figure size 1400x1000 with 2 Axes>,
 array([<Axes: title={'center': 'Histogram of Depth Spacing'}, xlabel='Depth Spacing (m)', ylabel='Frequency'>,
        <Axes: title={'center': 'Histogram of Time Spacing'}, xlabel='Time Spacing (s)', ylabel='Frequency'>],
       dtype=object))

[11]:

plots.plot_sampling_period_all(ds)

INFO:glidertest.utilities:Depth and time differences have been filtered to the middle 99% of values.
INFO:glidertest.utilities:Numeric median/mean/max/min values are based on the original data.

[11]:

(<Figure size 1400x1000 with 4 Axes>,
 array([<Axes: xlabel='Time Spacing (s)', ylabel='Frequency'>,
        <Axes: xlabel='Time Spacing (s)'>,
        <Axes: xlabel='Time Spacing (s)'>,
        <Axes: xlabel='Time Spacing (s)'>], dtype=object))

Check PROFILE_NUMBER assignment

We want to check if the profile number is increasing monotonically or not. The variable PROFILE_NUMBER is used to create an index of the profiles (dive and climb separately) within a glider dataset. Since it is used in many of the plotting functions below, it is worth checking that there aren’t obvious problems with the profile number assignment.

Note: When using pyglider to assign the profile, there can be some issues if the correct pressure and time threshold is not chosen. Users, based on their dataset will chose a specific min_dp and profile_min_time for the get_profiles_new(). If not chosen carefully, the function may struggle assigning a profile number (resulting in 0s in the middle of the dataset) or could split profiles when, for example, the glider get stuck on pycnoclne

[12]:

tools.check_monotony(ds.PROFILE_NUMBER)
plots.plot_prof_monotony(ds)

INFO:glidertest.utilities:PROFILE_NUMBER is always monotonically increasing

[12]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: xlabel='Time [UTC] (2023-Jun-04)', ylabel='Profile number'>)

Profile duration We can check for inconsistent profile duration as well which can highlight possible issues with how the profile number is assigned or weird navigation patterns

[13]:

duration = tools.compute_prof_duration(ds)
rolling_mean, overtime = tools.find_outlier_duration(duration,rolling=20, std=2)

fig, ax = plots.plot_outlier_duration(ds, rolling_mean, overtime, std=2)

INFO:glidertest.utilities:There are 4 profiles where the duration differs by 2 standard deviations of the nearby 20 profiles. Further checks are recommended

3. Sensors data overview

3.1 CTD

Sensors on gliders may suffer from biases (offsets) between when a glider is diving vs climbing. These plots average the dives and the climbs separately, to visualise whether a dataset may suffer from such a bias.

Note that there are two functions within glidertest which can be used to evaluate this bias: plot_updown_bias() uses a 2d gridding of the dataset, while quant_binavg() calculates a bin average of datapoints within a specified depth range.

See OG1 description of phase for more information on PHASE.

[14]:

updown_psal = tools.quant_updown_bias(ds, variable='PSAL', v_res=1)
updown_psal

[14]:

	dc	cd	dc_se	cd_se	dc_n	cd_n	depth
0	0.000359	NaN	0.000062	NaN	3	0	0.607222
1	-0.001706	NaN	0.001859	NaN	3	0	1.607222
2	0.006196	-0.010819	0.002059	0.000935	8	3	2.607222
3	0.010847	-0.007883	0.001334	0.002595	20	20	3.607222
4	0.011774	-0.008645	0.001531	0.002754	20	20	4.607222
...	...	...	...	...	...	...	...
77	-0.014688	0.011868	0.009004	0.004311	20	20	77.607226
78	-0.013137	0.010810	0.009321	0.006434	20	20	78.607226
79	-0.002660	0.000773	0.001752	0.002433	18	19	79.607226
80	-0.010982	0.001013	0.001308	0.002568	2	5	80.607226
81	NaN	NaN	NaN	NaN	0	0	81.607226

82 rows × 7 columns

[15]:

fig, ax = plt.subplots(1, 2, figsize=(15, 5))
plots.plot_updown_bias(ds, variable='TEMP', v_res=1, ax=ax[0])
plots.plot_updown_bias(ds, variable='PSAL', v_res=1, ax=ax[1])
plt.show()

3.2 Dissolved oxygen

Check for any possible drift in the data that might look suspicious. The great temporal and spatial variability may not allow for such check to be succesful. Evaluation using reference CTD cast data or any other data available in the study area is recommended.

[16]:

# Convert practical salinity and in-situ temperature
SA = gsw.SA_from_SP(ds["PSAL"], ds["PRES"], ds["LONGITUDE"], ds["LATITUDE"])
CT = gsw.CT_from_t(SA, ds["TEMP"], ds["PRES"])

o2sol = gsw.O2sol(SA, CT, ds["PRES"], ds["LONGITUDE"], ds["LATITUDE"])
ds["DOXY_SAT"] = 100 * ds["DOXY"] / o2sol
ds['DOXY_SAT']=ds['DOXY_SAT'].assign_attrs(units='%')

[17]:

plots.plot_updown_bias(ds, variable='DOXY', v_res=1)

[17]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: xlabel='Dissolved oxygen (mmol m$^{-}$³)', ylabel='Depth (m)'>)

[18]:

plots.check_temporal_drift(ds, variable='DOXY_SAT')

[18]:

(<Figure size 1400x1000 with 3 Axes>,
 array([<Axes: xlabel='Time [UTC] (2023-Jun-04)', ylabel='DOXY_SAT (%)'>,
        <Axes: xlabel='DOXY_SAT (%)', ylabel='Depth (m)'>], dtype=object))

[19]:

# Check any noise at surface during night time
plots.plot_daynight_avg(ds, variable='DOXY_SAT')

[#################################] 100% de421.bsp

[19]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: title={'center': '2023-06-04'}, xlabel='DOXY_SAT (%)', ylabel='Depth (m)'>)

Oxygen SOP provides a great summary of the needed processing of oxygen data and the various difference based on sensor model and platform type https://oceangliderscommunity.github.io/Oxygen_SOP/sections/oxygen_rtqc.html

We will now apply the check recommended by the SOP:

Global range check
Outlier and spike check
Stuck value test
Effects of biofouling
Noise during night time at surface

We use some functions from IOOS QC qartod tools [https://github.com/ioos/ioos_qc/blob/main/ioos_qc/qartod.py]

Global range check

[20]:

tools.compute_global_range(ds, variable='DOXY', min_val=-5, max_val=600)
plots.plot_global_range(ds, variable='DOXY', min_val=-5, max_val=600)

[20]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: title={'center': 'Global range check'}, xlabel='Dissolved oxygen (mmol m$^{-}$³)', ylabel='Frequency'>)

Outliers and spike check

[21]:

spike = qartod.spike_test( ds.DOXY,
    suspect_threshold= 25,
    fail_threshold= 50,
    method = "average",)
fig, ax = plots.plot_ioosqc(spike,suspect_threshold= [25], fail_threshold= [50], title='Spike test Dissolved oxygen')
plt.show()

Flat line test

[22]:

flat =qartod.flat_line_test(ds.DOXY,ds.TIME,1,2, 0.001)

# Add text with info
stuck_data = np.where(flat>2) # When it is not good
perct = (len(stuck_data[0]) * 100) / len(flat)
if len(stuck_data[0]) != 0:
    perct_0 = (len(np.where(ds.DOXY[stuck_data] == 0)[0]) * 100) / len(stuck_data[0])
    text = f'{np.round(perct, 1)}% of data is identical to the next one. Of those values, {np.round(perct_0, 1)}% are 0s (this is relevant in case your glider surveys \nanoxic waters)'
print(text)

fig, ax = plots.plot_ioosqc(flat,suspect_threshold= [0.01], fail_threshold=[0.1], title='Flat test Dissolved oxygen')

2.3% of data is identical to the next one. Of those values, 58.6% are 0s (this is relevant in case your glider surveys
anoxic waters)

Oxygen hysteresis

[23]:

df, diff, err_mean, err_range, rms = tools.compute_hyst_stat(ds, variable='DOXY', v_res=1)
fig, ax = plots.plot_hysteresis(ds, variable='DOXY')

[24]:

if len(err_mean.where(err_mean>2).dropna())<2:
    print('No data has an percentage error above 2% when the difference between up and downcast is computed')
else:
    print(f'{int((100*len(err_mean.where(err_mean>2).dropna()))/len(err_mean))}% of the data has a percentage error above 2 %. \nThe error reaches {int(np.nanmax(err_mean))}% at {int(df.depth[np.where(err_mean==np.nanmax(err_mean))[0][0]])}m. \nThis may idicate data shows hysteresis. \nThe RMS between climb and dive data is {np.round(rms, 1)}')

26% of the data has a percentage error above 2 %.
The error reaches 200% at 79m.
This may idicate data shows hysteresis.
The RMS between climb and dive data is 9.2

3.3 Bio-optics

Chlorophyll fluorescence

Check bottom data and see if we have stable data that can be used for calibration. We also check stability of data to assess whether or not we have suspicious drift over the mission
We check for any temporal drift when analysing the entire water column and not just the bottom data
We then check if data is affected by non photochemical quenching (NPQ). NPQ is a physiological response to high light environments used by plants and algae to protect themselves from damage and causes an evident weakening in fluorescence signal during the day. With the day_night_avg function, we compute day and night averages of chlorophyll. We then plot a selected section of chlorophyll data with plot_section_with_srss to see if any NPQ effect in the top few meters is visible and then we plot a selcted day daily and night average to check again any NPQ effect with plot_daynight_avg.

(Note: The test mission had issues with FLBBCD as it stopped working few days into the mission and got flooded)

[25]:

plots.process_optics_assess(ds, variable='CHLA')

INFO:glidertest.utilities:1.2% of scaled CHLA data is negative, consider recalibrating data
INFO:glidertest.utilities:Negative data in present from 2023-06-04T12:53 to 2023-06-04T21:33
INFO:glidertest.utilities:Negative data is present between 0.2 and 79.7
INFO:glidertest.utilities:CHLA data is missing for part of the mission

INFO:glidertest.utilities:Data from the deepest 10% of data has been analysed and data does not seem perfectly stable. An alternative solution for dark counts has to be considered.
Moreover, it is recommended to check the sensor has this may suggest issues with the sensor (i.e water inside the sensor, temporal drift etc)
INFO:glidertest.utilities:Data changed (increased or decreased) by 578.1% from the beginning to the end of the mission

[25]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: xlabel='Measurements', ylabel='Chlorophyll (mg m$^{-}$³)'>)

[26]:

plots.check_temporal_drift(ds, variable='CHLA')

[26]:

(<Figure size 1400x1000 with 3 Axes>,
 array([<Axes: xlabel='Time [UTC] (2023-Jun-04)', ylabel='Chlorophyll (mg m$^{-}$³)'>,
        <Axes: xlabel='Chlorophyll (mg m$^{-}$³)', ylabel='Depth (m)'>],
       dtype=object))

[27]:

# Let's visually check a section of chlorphyll and see if we observe any NPQ
plots.plot_quench_assess(ds, 'CHLA', ax=None, ylim=35);

[28]:

# Compute day and night average for chlorophylla and temeparture
dayT, nightT = tools.compute_daynight_avg(ds, variable='TEMP')
dayT

[28]:

	batch	depth	dat	date
0	0.0	0.0	12.055072	2023-06-04
1	0.0	1.0	12.065854	2023-06-04
2	0.0	2.0	11.984293	2023-06-04
3	0.0	3.0	11.843248	2023-06-04
4	0.0	4.0	11.793981	2023-06-04
...	...	...	...	...
77	0.0	77.0	8.474310	2023-06-04
78	0.0	78.0	8.478616	2023-06-04
79	0.0	79.0	8.492355	2023-06-04
80	0.0	80.0	8.504098	2023-06-04
81	0.0	81.0	8.505873	2023-06-04

82 rows × 4 columns

[29]:

fig, ax = plt.subplots(1, 3, figsize=(15, 5))

plots.plot_daynight_avg(ds, variable='TEMP', ax=ax[0])
plots.plot_daynight_avg(ds, variable='PSAL', ax=ax[1])
plots.plot_daynight_avg(ds, variable='CHLA', ax=ax[2])
plt.show()

Do we see any difference in chl between day and night? Can this just simply be explained by changes in water mass properties (different temp and salinity)?

Non photochemical quenching corrections

NPQ occurs only during the daytime, therefore when night measurements are available that are nearby in time and space, they can be used to correct daytime profiles. Different methods exist for NPQ correction. We recommend the paper by Thomalla et al. (2018) for a good overview on all possible methods.

GliderTools (https://glidertools.readthedocs.io/en/latest/optics.html#quenching-correction) provides good tools for correction.

Thomalla, S. J., W. Moutier, T. J. Ryan-Keogh, L. Gregor, and J. Schütt. 2018. An optimized method for correcting fluorescence quenching using optical backscattering on autonomous platforms. Limnology and Oceanography: Methods, 16: 132-144. DOI: https://doi.org/10.1002/lom3.10234

Optical Backscatter

Check bottom data and see if we have stable data that can be used for calibration. We also check stability of data to assess whether or not we have suspicious drift over the mission
We check for any temporal drift when analysing the entire water column and not just the bottom data
In case computation of particle backscattering from the scaled optical data was not done, this can be done following a function from GliderTools. this functions uses uses the coefficients from Zhang et al. (2009) to convert the raw counts into total backscatter (m-1), correcting for temperature and salinity. The \(\chi\) factor and \(\theta\) in this example were taken from Sullivan et al. (2013) and Slade & Boss (2015).

Slade, W., Boss, E. 2015. Spectral attenuation and backscattering as indicators of average particle size. Applied Optics 54: 7264-7277, doi:10.1364/AO.54.00726.

Sullivan, J., Twardowski, M., Zaneveld, J.R.V., Moore, C. 2013. Measuring optical backscattering in water. Light Scattering Reviews 7. 189-224. 10.1007/978-3-642-21907-8_6.

Zhang, X., and L. Hu. 2009. Estimating scattering of pure water from density fluctuation of the refractive index. Optics Express, 17: 1671-1678. DOI: 10.1364/OE.17.001671 7

[30]:

plots.process_optics_assess(ds, variable='BBP700')

INFO:glidertest.utilities:0.7% of scaled BBP700 data is negative, consider recalibrating data
INFO:glidertest.utilities:Negative data in present from 2023-06-04T12:58 to 2023-06-04T21:39
INFO:glidertest.utilities:Negative data is present between 37.5 and 58.8
INFO:glidertest.utilities:BBP700 data is missing for part of the mission

INFO:glidertest.utilities:Data from the deepest 10% of data has been analysed and data does not seem perfectly stable. An alternative solution for dark counts has to be considered.
Moreover, it is recommended to check the sensor has this may suggest issues with the sensor (i.e water inside the sensor, temporal drift etc)
INFO:glidertest.utilities:Data changed (increased or decreased) by 42.8% from the beginning to the end of the mission

[30]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: xlabel='Measurements', ylabel='Red backscatter, b${bp}$(700) (m$^{-}$¹)'>)

3.4 Photosyntetically Active Radiation (PAR)

[31]:

fig, ax = plt.subplots(1, 1, figsize=(5, 5))
plots.plot_updown_bias(ds, variable='DPAR', v_res=1, ax=ax)
plt.show()

Do we notice any strong up down cast bias?

Likely we do as the diving angle changes. The pitch for upcast and downcast are very different while the position of the sensor remains the same. This means that the angle at which the sensor is exposed to light is very different and data will not be comparable. Furthermore, navigation patterns have to be considered too when processing PAR data. As the glider sits at surface, the pitch (therefore the sensor angle) can be very different from the rest of the dive. Moreover, as the glider starts to dive or prepares for surfacing during a climb the pitch may be very different as well.

Discarding and reconstructing algebraically the surface PAR using an exponential equation and selecting data from only up or downcast is therefore recommended. GliderTools provides great examples and functions to address this issues (https://glidertools.readthedocs.io/en/latest/optics.html#par-replacement)

3.5 Vertical velocity

Vertical seawater velocity can be calculated from glider T&S data using the flight model. Note, the variables for glider flight from the flight model (GLIDER_VERT_VELO_MODEL) must be present in the dataset already.

[32]:

# Load glider dataset with vertical velocity parameters
ds_sg014 = fetchers.load_sample_dataset(dataset_name="sg014_20040924T182454_delayed_subset.nc") # Subset of a full Seaglider mission in the Labrador Sea

# Calculate vertical seawater velocity
# First, calculate the vertical speed of the glider from the depth data
ds_sg014 = tools.calc_w_meas(ds_sg014)

# Next, calculate the vertical seawater speed by differencing the DZDT data and the modelled vertical glider speed
ds_sg014 = tools.calc_w_sw(ds_sg014)

# Plot about 20 profiles to see the behaviour of the flight model
start_prof = 400
end_prof = 420
plots.plot_vertical_speeds_with_histograms(ds_sg014, start_prof, end_prof);

Downloading file 'sg014_20040924T182454_delayed_subset.nc' from 'https://erddap.observations.voiceoftheocean.org/examples/glidertest/sg014_20040924T182454_delayed_subset.nc' to '/home/runner/.cache/glidertest'.

[33]:

# One way to characterise the flight model performance is to look at the average
# vertical seawater velocity for climbs and for dives.  This uses a bin average of the original
# data on time intervals, and PHASE to separate dives from climbs

# Dive climb bias in the vertical speed
# Separate dives and climbs using the PHASE variable, PHASE=2 is a dive, PHASE=1 is a climb
ds_dives = ds_sg014.sel(N_MEASUREMENTS=ds_sg014.PHASE == 2)
ds_climbs = ds_sg014.sel(N_MEASUREMENTS=ds_sg014.PHASE == 1)

# Calculate bin averages
ds_out_dives = tools.quant_binavg(ds_dives, variable = 'VERT_CURR_MODEL', dz=10)
ds_out_climbs = tools.quant_binavg(ds_climbs, variable = 'VERT_CURR_MODEL', dz=10)

# Plot the profiles (compare to Fig 3 and Fig 4 in Frajka-Williams et al. 2011)
plots.plot_combined_velocity_profiles(ds_out_dives, ds_out_climbs)

[33]:

(<Figure size 1400x1000 with 1 Axes>,
 <Axes: xlabel='Vertical Velocity w$_{sw}$ (cm s$^{-1}$)', ylabel='Depth (m)'>)