Biomass density estimator

Generates a mask, as a bitmap, that indicates if the biomass density is low enough to obtain reliable in situ TS values.

This is accomplished by:

calculating the number of fish in the acoustic sampling volume (Nv, Sawada et al., 1993; see also Gauthier and Rose, 2001)
then calculating the ratio of multiple echoes when measuring a single target (M)
then testing the resulting values against their respective threshold values (nv_threshold and m_threshold respectively)

List of operands

Operand 1 - Single beam Sv mean (cell statistic)

Operand 2 - Single beam TS mean (cell statistic)

Operand 3 - Single beam Echo density (cell statistic)

See also: About the Cell Statistic operator, the Operator.eval method

Code

"""
Echoview Code Operator source file
========================================================================

Created by Echoview(R) 10.0.218 on Monday, 29 April 2019

See the Echoview help file for Code operator documentation and
examples.

NumPy User Guide: https://docs.scipy.org/doc/numpy/user/
NumPy Reference: https://docs.scipy.org/doc/numpy/Reference/
SciPy Reference Guide: https://docs.scipy.org/doc/scipy/Reference/

Echoview(R) is a registered trademark of Echoview Software Pty Ltd.
"""

# Authorship information
__author__ = "Echoview Software Pty Ltd. 2019."
__disclaimer__ = (
    "This example code is provided AS IS, without warranty of any "
    "kind, express or implied, including but not limited to the "
    "warranties of merchantability, fitness for a particular purpose "
    "and noninfringement. In no event shall Echoview Software Pty Ltd "
    "be liable for any claim, damages or other liability, arising "
    "from, out of or in connection with the use of this example code."
)
__version__ = "1.0"

# System Imports
from typing import List

# Libraries
import echoview as ev
import numpy as np


class Operator(ev.OperatorBase):
    """
    Biomass density estimator (Sawada index)
    ====================================================================

    Generates a mask (bitmap) that indicates if the biomass density is
    low enough to obtain reliable in situ TS values.

    Operands
    ---------------------

    * Operand 1 - Single beam Sv mean (cell statistic)
    * Operand 2 - Single beam TS mean (cell statistic)
    * Operand 3 - Single beam Echo density (cell statistic)

    Notes
    --------------------------

    This is accomplished by calculating the number of fish in the
    acoustic sampling volume (Nv, Sawada et al., 1993; see also
    Gauthier and Rose, 2001) and ratio of multiple echoes when measuring
    a single target (M) and then testing the resulting values against
    their respective threshold values (nv_threshold and m_threshold
    respectively).
    """

    def __init__(self):
        """
        Initialization: Define the M and Nv threshold values for use
        during the calculation.
        """

        self.m_threshold = 0.7
        self.nv_threshold = 0.04

    def result_type(self, input_types: List[ev.MeasurementType]):
        """
        Determines the type of variable from the list of input operand
        types.

        In this case we are generating a mask for singlebeam data,
        therefore we want SINGLE_BEAM_BOOLEAN.
        """

        return ev.MeasurementType.SINGLE_BEAM_BOOLEAN

    def eval(self, inputs: List[ev.OperandInput]):
        """
        Returns a boolean NumPy array representing which samples have a
        low enough biomass density for reliable in situ TS estimates.
        """

        # Get calibration from first operand (Sv mean), assuming this is
        # consistent across all three operands.
        ref_ping = inputs[0].measurement

        # Get PulseDuration in seconds.
        tau = ref_ping.cal('PulseDuration') / 1000.0
        # Get two way beam angle, converted to linear units.
        two_way_ba_lin = 10.0**(ref_ping.cal('TwoWayBeamAngle') / 10.0)
        c = ref_ping.cal('SoundSpeed')

        # Calculate an array of range values, one for each sample of the
        # current ping. This makes use of NumPy's linspace method, which
        # calculates an evenly spaced array between the start and stop
        # ranges of the ping.
        sample_count = ref_ping.data.shape[0]
        sample_thickness = ((ref_ping.stop_range - ref_ping.start_range)
                            / sample_count)
        ranges = np.linspace(ref_ping.start_range + (sample_thickness / 2),
                             ref_ping.stop_range-(sample_thickness / 2),
                             sample_count)

        # Get the ping samples from each of the three input operands,
        # fixing any invalid values (such as nan, inf) as we do so.
        sv_mean = np.ma.fix_invalid(inputs[0].measurement.data)
        ts_mean = np.ma.fix_invalid(inputs[1].measurement.data)
        nec = np.ma.fix_invalid(inputs[2].measurement.data)

        # Convert the Sv and TS mean values to linear units.
        sv_mean = to_linear(sv_mean)
        ts_mean = to_linear(ts_mean)

        # Calculate fish density from Sv and TS mean given echo
        # integrator.
        nei = sv_mean/ts_mean

        # Calculate number of fish per acoustic sampling volume (Nv),
        # sometimes known as the Sawada index.
        # (nv= c*t*equivBA* range^2 * nei)/2
        nv = (c*tau*two_way_ba_lin*np.square(ranges)*nei) / 2

        # Ratio for multiple echoes when measuring a single target.
        M = ((nei - nec) / nei)

        # Calculate the boolean mask to return using NumPy's logical_and
        # method to evaluate samples where nv < nv Threshold and
        # M < M Threshold.
        # Note: Echoview doesn't support nodata for boolean variables so
        # set any masked out values to False.
        return np.logical_and(nv < self.nv_threshold, M < self.m_threshold) \
                 .filled(False)


def to_linear(data: np.ma.MaskedArray, out=None):
    """Converts the specified input to linear values"""

    # Equivalent to 10 ** (data / 10) but more efficient due to reusing
    # the same output array instead of making multiple copies.
    if out is not None and out.shape != data.shape:
        raise ValueError(
            'out must be a NumPy array with the same shape as the input array')

    out = np.divide(data, 10, out=out)
    out = np.power(10, out, out=out)

    return out