###############
Post-processing
###############

Overview
========
RSKtools is RBR's open source MATLAB toolbox for reading, post-processing, visualizing, and exporting RBR logger data.
Users may plot data as a time series or as depth profiles using tailored plotting utilities. Time-depth heat maps can be
plotted easily to visualise transects or moored profiler data. A full suite of data post-processing functionality, such as
functions to match sensor time constants and bin average, are available to enhance data quality.

You may review the :doc:`getting started guide </guides/getting-started-guide>` for an introduction on how to open an RSK file and access data.

RSKtools includes a series of functions to post-process RBR logger data. Below we show how to implement some common
processing steps to obtain the highest quality data possible. All post-processing functions are customisable via
name-value pair input arguments. Documentation for each function can be accessed in MATLAB using ``help``.

Note on multi-schedule RSK files
================================
Unlike derive functions, post-processing functions are schedule-specific.
Their default behaviour is to work on the auto-detected default schedule (``'schedule', 'default'``),
which selects the schedule with the most Conductivity data, then Temperature, then schedule 1.
If your RSK has multiple schedules, you can specify a schedule label explicitly
(e.g., ``'schedule', 's.CTD'``). See the :ref:`getting started guide <multi-schedule-support>`
for more information on multi-schedule support.


Basic Usage
============

.. code-block:: matlab

    rsk = RSKopen('sample.rsk');

    % Print a list of all the channels in the RSK file
    RSKprintchannels(rsk);

    % Read profiles from the RSK file
    rsk = RSKreadprofiles(rsk);

    % Derive sea pressure from total pressure
    rsk = RSKderiveseapressure(rsk);

    % Plot a few downcasts of temperature, conductivity, and chlorophyll
    RSKplotprofiles(rsk, 'channel', {'Conductivity', 'Temperature', 'Chlorophyll'}, ...
                    'profile', [1 2 3], 'direction', 'down');

.. code-block:: text

    Model:           RBRconcerto
    Serial ID:       80231
    Sampling period: 0.167 second
        index        channel             unit
        _____    ________________    _____________
        1        Conductivity        mS/cm
        2        Temperature         °C
        3        Pressure            dbar
        4        Dissolved O2        %
        5        Turbidity           NTU
        6        PAR                 µMol/m²/s
        7        Chlorophyll         µg/l

.. figure:: /img/post_processing_basic_use.png
    :scale: 70%
    :alt: plot profiles example


Derive sea pressure from total pressure
=======================================

We suggest deriving sea pressure first, especially when an atmospheric pressure other than the nominal
value of 10.1325 dbar is desired. The ``patm`` argument is the atmospheric pressure used to calculate the sea pressure.
A custom value can be used; otherwise, the default is to retrieve the value stored in the parameters field
or to assume it is 10.1325 dbar if the parameter's field is unavailable.
In this example, we take atmospheric pressure to be 10 dbar:

.. code-block:: matlab

    rsk = RSKderiveseapressure(rsk, 'patm', 10);


Data post-processing
====================

What follows is a generic recipe for post-processing RBR CTD data. It is a guideline -- RBR CTDs are
used in many environments, with many sensor packages, and are profiled from a variety of vessels
(or no vessel at all). While the basic approach here is relevant for many cases, the processing parameters
may not apply widely.

First, keep a copy of the raw data to compare with the processed data later.

.. code-block:: matlab

    raw = rsk;

Correct for A2D zero-order hold
===============================

The analogue-to-digital (A2D) converter on RBR instruments must recalibrate periodically. In the time it
takes for the calibration to finish, one or more samples are missed. The instrument firmware fills the missed
sample with the same data measured during the previous sample, a technique called a zero-order hold.
:func:`RSKcorrecthold` identifies zero-hold points by finding where consecutive differences of each channel
are equal to zero and then replaces these samples with a NaN or an interpolated value.

.. code-block:: matlab

    rsk = RSKcorrecthold(rsk, 'action', 'interp');


Low-pass filtering
==================

Low-pass filtering is commonly used to reduce noise and to match sensor time constants, typically for
temperature and conductivity. Users may also wish to filter other channels to simply reduce noise (e.g.,
optical channels such as chlorophyll-a or turbidity).

Most RBR instruments designed for profiling are equipped with thermistors that have a time constant of 100 ms,
which is "slower" than the conductivity cell. When the time constants are different, salinity will contain spikes
at strong gradients. One solution is to "slow down" the conductivity sensor to match the thermistor. In this example
dataset, the logger sampled at 6 Hz, so a 5-sample running average provides more than sufficient smoothing to match
the time response of the conductivity sensor to the thermistor.

.. code-block:: matlab

    rsk = RSKsmooth(rsk, 'channel', 'Salinity', 'windowLength', 5);

.. figure:: /img/post_processing_LPfilter.png
    :scale: 90%
    :alt: plot low-pass filter example


Alignment of conductivity and temperature
=========================================

Conductivity and temperature often need to be aligned in time to account for the fact these sensors are not
always co-located on the logger. The implication is that, under dynamic conditions (e.g., profiling), the
sensors are measuring a slightly different parcel of water at any instant.

Furthermore, sensors with long time constants introduce a time lag to the data. For example, dissolved
oxygen sensors often have a long time constant, and this delays the measurement relative to the true value.
This can be fixed to some degree by advancing the sensor data in time.

When temperature and conductivity are misaligned, salinity will contain spikes at sharp interfaces and
a bias in continuously stratified environments. Properly aligning the sensors, together with matching the
time response, will minimise salinity spiking and bias.

A common approach to determine the optimal lag is to compute and plot salinity for a range of lags and
choose the lag (often by eye) with the smallest salinity spikes at sharp temperature interfaces.
RSKtools includes :func:`RSKcalculateCTlag` that estimates the optimal lag between conductivity and
temperature by minimizing salinity spiking. We currently suggest using both approaches to check for consistency.

As a rough guide, temperature from a CTD equipped with the red combined CT cell and a fast thermistor
typically requires only a very small-time advance (perhaps tens of milliseconds). Temperature from a CTD
equipped with a cylindrical black conductivity cell (with the thermistor on the sensor endcap) typically
requires a temperature lag correction of about 0.1 s to 0.3 s (1 or 2 samples at 6 Hz).

.. code-block:: matlab

    % Required shift of C relative to T for each profile
    lag = RSKcalculateCTlag(rsk, 'seapressureRange', [1 50], 'direction', 'down');
    % Advance temperature by the median lag
    lag = -median(lag);
    rsk = RSKalignchannel(rsk, 'channel', 'Temperature', 'lag', lag, 'direction', 'down');

Users wishing to learn more about dynamic sensor corrections and RBR CTDs are encouraged to watch a
special RBR webinar on dynamic errors from May 2020 (`Youtube`_ and `PDF`_).


Remove loops
============

Working in rough seas can cause the CTD profiling rate to vary, and even change signs (i.e., the CTD
momentarily changes direction). When this happens, the CTD effectively samples its own wake,
degrading the quality of the profile in regions of strong gradients. The measurements taken when the instrument
is profiling too slowly or during a pressure reversal should not be used for further analysis.
We recommend using :func:`RSKremoveloops` to flag and treat the data when the instrument (1) falls below a
threshold speed and (2) when the pressure reverses (the CTD "loops"). Before using :func:`RSKremoveloops`, use
:func:`RSKderiveseapressure` to calculate sea pressure from total pressure, :func:`RSKderivedepth` to calculate depth
from sea pressure, and then use :func:`RSKderivevelocity` to calculate profiling rate.

.. code-block:: matlab

    rsk = RSKderiveseapressure(rsk);
    rsk = RSKderivedepth(rsk);
    rsk = RSKderivevelocity(rsk);
    % Apply the algorithm
    rsk = RSKremoveloops(rsk, 'threshold', 0.3);

.. figure:: /img/post_processing_RemoveLoops.png
    :scale: 70%
    :alt: plot removeloops example


Derived variables
=================

RSKtools includes a few convenience functions to derive common oceanographic variables. For
example, :func:`RSKderivesalinity` computes Practical Salinity using the TEOS-10 GSW function
``gsw_SP_from_C``, and adds a new channel called ``Salinity`` to the data. The `TEOS-10 GSW`_ MATLAB
toolbox must be installed separately.

.. code-block:: matlab

    rsk = RSKderiveseapressure(rsk);
    rsk = RSKderivedepth(rsk);
    rsk = RSKderivevelocity(rsk);
    rsk = RSKderivesalinity(rsk);
    rsk = RSKderivesigma(rsk);

    % Print a list of channels
    RSKprintchannels(rsk);

.. code-block:: text

    Model:           RBRconcerto
    Serial ID:       80231
    Sampling period: 0.167 second
        index          channel              unit
        _____    ___________________    _____________
        1        Conductivity           mS/cm
        2        Temperature            °C
        3        Pressure               dbar
        4        Dissolved O2           %
        5        Turbidity              NTU
        6        PAR                    µMol/m²/s
        7        Chlorophyll            µg/l
        8        Sea Pressure           dbar
        9        Salinity               PSU
        10       Depth                  m
        11       Velocity               m/s
        12       Density Anomaly        kg/m³


Bin average all channels by sea pressure
========================================

Bin averaging reduces sensor noise and ensures that each profile is referenced to a common grid. The latter
is often an advantage for plotting data as "heatmaps." :func:`RSKbinaverage` allows users to bin channels
according to any reference, but the most common choices are time, depth, and sea pressure. It also can
handle grids with a variable bin size. In the following example, the data are averaged into 0.25 dbar bins.

.. code-block:: matlab

    rsk = RSKbinaverage(rsk, 'binBy', 'Sea Pressure', 'binSize', 0.25, ...
                        'boundary', [0.5 5.5], 'direction', 'up');

.. figure:: /img/post_processing_binaverage.png
    :scale: 70%
    :alt: plot binaverage example


Compare the raw and processed data
==================================

Use :func:`RSKplotprofiles` to compare the binned data to the raw data for a few example profiles.

.. code-block:: matlab

    figure;
    RSKplotprofiles(raw, 'channel', {'Salinity'}, 'profile', 1, 'direction', 'down');
    hold on;
    RSKplotprofiles(rsk, 'channel', {'Salinity'}, 'profile', 1, 'direction', 'down');
    hold off;

.. figure:: /img/post_processing_compare.png
    :scale: 70%
    :alt: plot comparison example


Visualise data with a 2D plot
=============================

:func:`RSKimages` generates a time/depth heat-map of a channel. The x-axis is time; the y-axis is a reference
channel (default is sea pressure). All the profiles must be evaluated on the same reference channel
grid, which is accomplished with :func:`RSKbinaverage`. For more details, see :func:`RSKbinaverage` and
:func:`RSKgenerate2D`.

.. code-block:: matlab

    RSKimages(rsk, 'channel', {'Temperature', 'Salinity'}, 'direction', 'up');

.. figure:: /img/post_processing_images.png
    :scale: 70%
    :alt: plot 2D example


Export logger data to CSV files
===============================

:func:`RSK2CSV` writes logger data and metadata to one or more CSV files. The CSV files contain a header
with important logger metadata and a record of the processing steps made by
RSKtools. The data table starts with a row of variable names and units above each column of channel data.
If the ``profile`` number is specified, then one file will be written for each profile.
Furthermore, an extra column called ``cast_direction`` will be included. Users can select
which channels and profiles to write, the output directory, and also specify additional comments to be
placed after the metadata in the file header.

.. code-block:: matlab

    RSK2CSV(rsk, 'channel', {'Depth', 'Temperature', 'Salinity'}, ...
            'profile', [1 2 3], 'comment', 'Processed with RSKtools');

RSKtools also has an export function to write Ruskin RSK files (see :func:`RSK2RSK`).


Display a summary of all the processing steps
=============================================

.. code-block:: matlab

    disp(rsk.log);


Other Resources
===============

We recommend reading the:

* :doc:`RSKtools API reference </_rsk/read>` for detailed function documentation.

* :doc:`RSKtools getting started guide </guides/getting-started-guide>` for an introduction on how to open an RSK file, make plots, and access the data.

About this document
===================

The full source code for this documentation is located in the ``docs`` directory of the repository.
All documentation is built and generated using `Sphinx`_ with the `sphinxcontrib-matlabdomain`_ extension.

.. _Youtube: https://www.youtube.com/watch?v=uJXh5wV_2NE
.. _PDF: https://rbr-global.com/wp-content/uploads/2020/05/6.-2020-05-06-RBR-Webinar-CTD-Dyn.-Corr.pdf
.. _TEOS-10 GSW: https://teos-10.org/
.. _Sphinx: https://www.sphinx-doc.org/en/master/
.. _sphinxcontrib-matlabdomain: https://github.com/sphinx-contrib/matlabdomain
