Skip to content

API reference

This is the API reference for the echoSMs package.

Each type of model is contained in a separate Python class (with name ending in Model), but with common calling signatures across all model classes, as defined in ScatterModelBase. There are also classes to provide ready access to the benchmark models and reference model definitions. There are also utility functions.

echosms.ScatterModelBase()

Bases: ABC

Base class for a class that provides a scattering model.

All scattering models should inherit from this class, have a name that ends with 'Model', and provide initialisation and calculate_ts_single() functions.

Initialise.

ATTRIBUTE DESCRIPTION
long_name

The long name of the model.

TYPE: str

short_name

A short version of the model's long name, typically an acronym.

TYPE: str

analytical_type

Whether the model implements an exact or an approximate model.

TYPE: str

boundary_types

The types of boundary conditions that the model provides, e.g., 'fixed rigid', 'pressure release', 'fluid filled'

TYPE: list[boundary_type]

shapes

The target shapes that the model can represent.

TYPE: list[str]

max_ka

An approximate maximum ka value that will result in accurate target strength results. Note that ka is often not the only parameter that determines the accuracy of the model (e.g., aspect ratio and incident angle can also affect the accuracy).

TYPE: float

no_expand_parameters

The model parameters that are not expanded into Pandas DataFrame columns or Xarray DataArray coordinates. They will instead end up as a dict in the DataFrame or DataArray attrs attribute.

TYPE: list[str]

Source code in src/echosms/scattermodelbase.py 
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
@abc.abstractmethod
def __init__(self):
    """Initialise.

    Attributes
    ----------
    long_name : str
        The long name of the model.
    short_name : str
        A short version of the model's long name, typically an acronym.
    analytical_type : str
        Whether the model implements an ``exact`` or an ``approximate`` model.
    boundary_types : list[boundary_type]
        The types of boundary conditions that the model provides, e.g., 'fixed rigid',
        'pressure release', 'fluid filled'
    shapes : list[str]
        The target shapes that the model can represent.
    max_ka : float
        An approximate maximum ka value that will result in accurate target strength results.
        Note that ka is often not the only parameter that determines the accuracy of the
        model (e.g., aspect ratio and incident angle can also affect the accuracy).
    no_expand_parameters : list[str]
        The model parameters that are not expanded into Pandas DataFrame columns or
        Xarray DataArray coordinates. They will instead end up as a dict in the DataFrame or
        DataArray `attrs` attribute.
    """
    self.long_name = ''
    self.short_name = ''
    self.analytical_type = ''
    self.boundary_types = []
    self.shapes = []
    self.max_ka = np.nan
    self.no_expand_parameters = []

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single() abstractmethod

Calculate the TS for one parameter set.

Source code in src/echosms/scattermodelbase.py 
248
249
250
@abc.abstractmethod
def calculate_ts_single(self):
    """Calculate the TS for one parameter set."""

validate_parameters(p) abstractmethod

Validate the model parameters.

PARAMETER DESCRIPTION
p

The model parameters.

TYPE: dict | DataFrame | DataArray

RAISES DESCRIPTION
ValueError

If any of the model parameters are invalid.
KeyError

If any required model parameters are not present.
Source code in src/echosms/scattermodelbase.py 
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
@abc.abstractmethod
def validate_parameters(self, p: dict | pd.DataFrame | xr.DataArray):
    """Validate the model parameters.

    Parameters
    ----------
    p :
        The model parameters.

    Raises
    ------
    ValueError
        If any of the model parameters are invalid.
    KeyError
        If any required model parameters are not present.
    """

echosms.DCMModel()

Bases: ScatterModelBase

Modal series deformed cylinder model (DCM).

This class contains methods to calculate acoustic scatter from finite straight cylinders with various boundary conditions.

Source code in src/echosms/dcmmodel.py 
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'deformed cylinder model'
    self.short_name = 'dcm'
    self.analytical_type = 'approximate analytical'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled]
    self.shapes = ['finite cylinder']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, b, theta, f, boundary_type, target_c=None, target_rho=None, validate_parameters=True, **kwargs)

Calculate the scatter from a finite cylinder using the modal series deformed cylinder model.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Radius of the cylinderical target [m].

TYPE: float

b

Length of the cylinderical target [m].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The model type. Supported model types are given in the boundary_types class attribute.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the sphere [m/s]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

target_rho

Density of the fluid inside the sphere [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Section B.1 of Jech et al. (2015).

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742–3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/dcmmodel.py 
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
def calculate_ts_single(self, medium_c, medium_rho, a, b, theta, f, boundary_type: bt,
                        target_c=None, target_rho=None, validate_parameters=True,
                        **kwargs):
    """
    Calculate the scatter from a finite cylinder using the modal series deformed cylinder model.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    a : float
        Radius of the cylinderical target [m].
    b : float
        Length of the cylinderical target [m].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The model type. Supported model types are given in the `boundary_types` class attribute.
    target_c : float, optional
        Sound speed in the fluid inside the sphere [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho : float, optional
        Density of the fluid inside the sphere [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Section B.1 of Jech et al. (2015).

    References
    ----------
    Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B.,
    Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015.
    Comparisons among ten models of acoustic backscattering used in aquatic ecosystem
    research. Journal of the Acoustical Society of America 138, 3742–3764.
    <https://doi.org/10.1121/1.4937607>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if theta == 0.0:
        return nan

    theta_rad = theta*pi/180.
    kL = wavenumber(medium_c, f)*b
    K = wavenumber(medium_c, f) * sin(theta_rad)
    Ka = K*a

    m = range(30)  # TODO this needs to vary with f

    match boundary_type:
        case bt.fixed_rigid:
            series = map(lambda m: (-1)**m * Neumann(m)*(jvp(m, Ka) / h1vp(m, Ka)), m)
        case bt.pressure_release:
            series = map(lambda m: (-1)**m * Neumann(m)*(jv(m, Ka) / hankel1(m, Ka)), m)
        case bt.fluid_filled:
            g = target_rho/medium_rho
            h = target_c/medium_c
            gh = g*h
            Kda = K/h*a

            def Cm(m):
                numer = (jvp(m, Kda)*yv(m, Ka)) / (jv(m, Kda)*jvp(m, Ka))\
                    - gh*(yvp(m, Ka)/jvp(m, Ka))
                denom = (jvp(m, Kda)*jv(m, Ka)) / (jv(m, Kda)*jvp(m, Ka)) - gh
                return numer/denom

            series = map(lambda m: 1j**(2*m) * Neumann(m) / (1 + 1j*Cm(m)), m)
        case _:
            raise ValueError(f'The {self.long_name} model does not support '
                             f'a model type of "{boundary_type}".')

    fbs = 1j*b/pi * (sin(kL*cos(theta_rad)) / (kL*cos(theta_rad))) * sum(series)
    return 20*log10(abs(fbs))  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/dcmmodel.py 
26
27
28
29
30
31
32
33
34
35
36
37
38
39
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'a', 'b', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        if t == bt.fluid_filled:
            super()._present_and_positive(p, ['target_c', 'target_rho'],
                                          mask=p['boundary_type'] == t)

DWBA models

There are several models that use the distorted-wave Born approximation, documented below. There are also some functions to make cylinder and spheroid shapes for use in the DWBA models.

echosms.DWBAModel()

Bases: ScatterModelBase

Distorted-wave Born approximation (DWBA) scattering models.

This class calculates acoustic scatter from piecewise cylindical shapes with weakly scattering material contrasts. It can also add a stochastic component to provide the stochastic DWBA (SDWBA) model.

Source code in src/echosms/dwbamodel.py 
21
22
23
24
25
26
27
28
29
30
31
32
33
def __init__(self):
    super().__init__()
    self.long_name = 'distorted-wave Born approximation'
    self.short_name = 'dwba'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['piecewise cylindical']
    self.max_ka = 20
    # The distorted wave Born approximation is increasingly inaccurate outside these limits:
    self.g_range = [0.95, 1.05]
    self.h_range = [0.95, 1.05]
    self.no_expand_parameters = ['a', 'rv_pos', 'rv_tan']
    self.rng = np.random.default_rng()  # for SDWBA

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, theta, phi, f, target_c, target_rho, a, rv_pos, rv_tan, phase_sd=0, num_runs=1, validate_parameters=True, **kwargs)

Distorted-wave Born approximation scattering model.

Implements the distorted-wave Born approximation and stochastic DWBA model for calculating the acoustic backscatter from weakly scattering bodies.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to run the model at [Hz]

TYPE: float

target_c

Sound speed in the fluid inside the target [m/s].

TYPE: float

target_rho

Density of the fluid inside the target [kg/m³].

TYPE: float

a

The radii of the discs that define the target shape [m].

TYPE: iterable

rv_pos

An interable of vectors of the 3D positions of the centre of each disc that defines the target shape. Each vector should have three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: iterable[ndarray]

rv_tan

An interable of unit vectors of the tangent to the target body axis at the points given in rv_pos. Each vector should have three values corresponding to the x, y, and z components of the tangent vector.

TYPE: iterable[ndarray]

phase_sd

If non-zero, this model becomes the SDWBA (stochastic DWBA). A random phase is applied to each term in the DWBA integral, obtained from a Gaussian distribution centered on 0 with standard deviation of phase_sd [°].

TYPE: float DEFAULT: 0

num_runs

The number of times to run the SDWBA model. The mean TS (calculated in the linear domain) is returned. Intended to be used in conjunction with phase_sd.

TYPE: int DEFAULT: 1

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) [dB] of the target.
Notes

This class implements the method presented in Eqn (5) of Stanton et al. (1998). This caters to objects that are discretised into a set of adjacent discs with potentially offset centres.

The DWBA model density and sound speed values are often specified as ratios of the medium and target values (g & h). However, to maintain compatibility with other echoSMs models this implementation requires actual densities and sound speeds.

The SDWBA component is as per Demer & Conti (2003).

References

Demer, D. A., & Conti, S. G. (2003). Reconciling theoretical versus empirical target strengths of krill: Effects of phase variability on the distorted-wave Born approximation. ICES Journal of Marine Science, 60, 429-434. https://doi.org/10.1016/S1054-3139(03)00002-X

Stanton, T. K., Chu, D., & Wiebe, P. H. (1998). Sound scattering by several zooplankton groups. II. Scattering models. The Journal of the Acoustical Society of America, 103(1), 236-253. https://doi.org/10.1121/1.421110

Source code in src/echosms/dwbamodel.py 
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def calculate_ts_single(self, medium_c, medium_rho, theta, phi, f, target_c, target_rho,
                        a, rv_pos, rv_tan, phase_sd=0, num_runs=1, validate_parameters=True,
                        **kwargs) -> float:
    """Distorted-wave Born approximation scattering model.

    Implements the distorted-wave Born approximation and stochastic DWBA
    model for calculating the acoustic backscatter from weakly scattering bodies.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to run the model at [Hz]
    target_c : float
        Sound speed in the fluid inside the target [m/s].
    target_rho : float
        Density of the fluid inside the target [kg/m³].
    a : iterable
        The radii of the discs that define the target shape [m].
    rv_pos : iterable[np.ndarray]
        An interable of vectors of the 3D positions of the centre of each disc that
        defines the target shape. Each vector should have three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan : iterable[np.ndarray]
        An interable of unit vectors of the tangent to the target body axis at
        the points given in `rv_pos`. Each vector should have three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    phase_sd : float
        If non-zero, this model becomes the SDWBA (stochastic DWBA). A random phase is
        applied to each term in the DWBA integral, obtained from a Gaussian distribution
        centered on 0 with standard deviation of `phase_sd` [°].
    num_runs : int
        The number of times to run the SDWBA model. The mean TS (calculated in the
        linear domain) is returned. Intended to be used in conjunction with `phase_sd`.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) [dB] of the target.

    Notes
    -----
    This class implements the method presented in Eqn (5) of Stanton et al. (1998). This
    caters to objects that are discretised into a set of adjacent discs with potentially
    offset centres.

    The DWBA model density and sound speed values are often specified as ratios of the medium
    and target values (g & h). However, to maintain compatibility with other echoSMs models
    this implementation requires actual densities and sound speeds.

    The SDWBA component is as per Demer & Conti (2003).

    References
    ----------
    Demer, D. A., & Conti, S. G. (2003). Reconciling theoretical versus empirical target
    strengths of krill: Effects of phase variability on the distorted-wave Born approximation.
    ICES Journal of Marine Science, 60, 429-434.
    <https://doi.org/10.1016/S1054-3139(03)00002-X>

    Stanton, T. K., Chu, D., & Wiebe, P. H. (1998). Sound scattering by several zooplankton
    groups. II. Scattering models. The Journal of the Acoustical Society of America, 103(1),
    236-253. <https://doi.org/10.1121/1.421110>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    do_sdwba = False if phase_sd == 0.0 else True

    # The structure of this code follows closely the formulae in Stanton et al (1998). Where
    # relevant, the equation numbers from that paper are given.

    k1 = wavenumber(medium_c, f)
    k2 = wavenumber(target_c, f)
    kappa_1 = 1.0 / (medium_rho * medium_c**2)  # Eqn (4) for medium
    kappa_2 = 1.0 / (target_rho * target_c**2)  # Eqn (4) for target
    gamma_k = (kappa_2 - kappa_1) / kappa_1  # Eqn (2)
    gamma_rho = (target_rho - medium_rho) / target_rho  # Eqn (3)
    # Note: the (gamma_k - gamma_rho) term in Eqn (5) can be simplified using g & h to:
    # 1/gh^2 + 1/g - 2.

    # Calculate the distance between each disc using the disc position vectors. The Euclidean
    # distance is the L2 norm so use np.norm(). Could also use
    # scipy.spatial.distance.euclidean(), but that is apparently a lot slower.
    dist = np.array([np.linalg.norm(r1-r0) for r0, r1 in zip(rv_pos[0:-1], rv_pos[1:])])  # [m]

    # Thickness of each disc based on the distance between discs. The first and last
    # discs are treated differently.
    d_rv_pos = np.hstack((dist[0], dist[0:-1]/2 + dist[1:]/2, dist[-1]))  # [m]

    # Calculate direction of incident wave given theta and phi. The echoSMs convention
    # has the target rotating and the incident vector always being (0,0,1), but for the DWBA
    # we keep the body stationary and change the incident vector.
    rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
    k_hat_i = rot.apply([0, 0, 1])  # needs to be a unit vector
    k_i2 = k_hat_i * k2  # incident vector with magnitude equal to wavenumber inside the target

    # This code is a little complex because it does both the DWBA and SDWBA
    phase_factors = np.ones(len(a))  # for DWBA
    runs = []

    for run in range(int(num_runs)):  # is only ever 1 run for the DWBA
        if do_sdwba:
            phase_factors = np.exp(1j*self.rng.normal(scale=radians(phase_sd), size=len(a)))

        integral = 0.0
        for a_, r_pos, d_r_pos, r_tan, p in zip(a, rv_pos, d_rv_pos, rv_tan, phase_factors):
            # The round() is here because sometimes the dot product gets values slightly outside
            # the [-1, 1] range (due to floating point inaccuracies) and cos will complain.
            cbeta_tilt = cos(pi/2 - acos(round(k_hat_i @ r_tan, 8)))

            # This is the integral part of Eqn (5) with addition of
            # Eqn (4) from Demer & Conti (2003) for the SDWBA part
            integral += (gamma_k-gamma_rho) * exp(2j*(k_i2@r_pos))\
                * a_*j1(2*k2*a_*cbeta_tilt) / cbeta_tilt * abs(d_r_pos) * p
        runs.append(integral)

    return 20*log10(np.mean(abs(k1/4.0*np.array(runs))))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/dwbamodel.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'target_rho', 'target_c', 'f'])

    g = p['target_rho'] / p['medium_rho']
    h = p['target_c'] / p['medium_c']

    if (np.any(g < self.g_range[0])) or np.any((g > self.g_range[1])):
        warnings.warn('Ratio of target and medium densities (g) is outside the DWBA limits.')

    if (np.any(h < self.h_range[0])) or np.any((h > self.h_range[1])):
        warnings.warn('Ratio of target and medium sound speeds (h) are '
                      'outside the DWBA limits.')

    if not np.all([isclose(1.0, np.linalg.norm(v)) for v in p['rv_tan']]):
        raise ValueError('All vectors in rv_tan must be of unit length.')

echosms.PTDWBAModel()

Bases: ScatterModelBase

Phase-tracking distorted-wave Born approximation scattering model.

Source code in src/echosms/ptdwbamodel.py 
13
14
15
16
17
18
19
20
21
def __init__(self):
    super().__init__()
    self.long_name = 'phase-tracking distorted-wave Born approximation'
    self.short_name = 'pt-dwba'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['unrestricted voxel-based']
    self.max_ka = 20
    self.no_expand_parameters = ['volume', 'voxel_size', 'rho', 'c']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(volume, theta, phi, f, voxel_size, rho, c, validate_parameters=True, **kwargs)

Phase-tracking distorted-wave Born approximation scattering model.

Implements the phase-tracking distorted-wave Born approximation model for calculating the acoustic backscatter from weakly scattering bodies.

PARAMETER DESCRIPTION
volume

The object to be modelled as a 3D volume of voxels. Array contents should be 0 for the surrounding medium, then increasing by 1 for each additional material type (i.e., 1, 2, 3, etc). volume should be arranged as per the echoSMs coordinate system, where

  • axis 0 (rows) is the x-axis
  • axis 1 (columns) is the y-axis
  • axis 2: (slices) is the z-axis

Increasing axes indices correspond to increasing x, y, and z values.

TYPE: Numpy ndarray[int]

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to run the model at [Hz]

TYPE: float

voxel_size

The size of the voxels in volume [m], ordered (x, y, z). This code assumes that the voxels are cubes so y and z are currently irrelevant.

TYPE: iterable[float]

rho

Densities of each material. Must have at least the same number of entries as unique integers in volume [kg/m³].

TYPE: iterable[float]

c

Sound speed of each material. Must have at least the same number of entries as unique integers in volume [m/s].

TYPE: iterable[float]

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) [dB] of the target.
Notes

This class implements the method presented in Jones et. al. (2009). The code is based closely on the Matlab code in Jones (2006).

References

Jones, B. A. (2006). Acoustic scattering of broadband echolocation signals from prey of Blainville's beaked whales: Modeling and analysis. Master of Science, Massachusetts Institute of Technology. https://doi.org/10.1575/1912/1283

Jones, B. A., Lavery, A. C., & Stanton, T. K. (2009). Use of the distorted wave Born approximation to predict scattering by inhomogeneous objects: Application to squid. The Journal of the Acoustical Society of America, 125(1), 73-88. https://doi.org/10.1121/1.3021298

Source code in src/echosms/ptdwbamodel.py 
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
def calculate_ts_single(self, volume, theta, phi, f, voxel_size, rho, c,
                        validate_parameters=True, **kwargs):
    """Phase-tracking distorted-wave Born approximation scattering model.

    Implements the phase-tracking distorted-wave Born approximation
    model for calculating the acoustic backscatter from weakly scattering bodies.

    Parameters
    ----------
    volume : Numpy ndarray[int]
        The object to be modelled as a 3D volume of voxels. Array contents should be 0
        for the surrounding medium, then increasing by 1 for each additional material
        type (i.e., 1, 2, 3, etc). `volume` should be arranged as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems), where

        - axis 0 (rows) is the _x_-axis
        - axis 1 (columns) is the _y_-axis
        - axis 2: (slices) is the _z_-axis

        Increasing axes indices correspond to increasing _x_, _y_, and _z_ values.

    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].

    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].

    f : float
        Frequency to run the model at [Hz]

    voxel_size : iterable[float]
        The size of the voxels in `volume` [m], ordered (_x_, _y_, _z_).
        This code assumes that the voxels are cubes so _y_ and _z_ are currently irrelevant.

    rho : iterable[float]
        Densities of each material. Must have at least the same number of entries as unique
        integers in `volume` [kg/m³].

    c : iterable[float]
        Sound speed of each material. Must have at least the same number of entries as unique
        integers in `volume` [m/s].
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) [dB] of the target.

    Notes
    -----
    This class implements the method presented in Jones et. al. (2009). The code is
    based closely on the Matlab code in Jones (2006).

    References
    ----------
    Jones, B. A. (2006). Acoustic scattering of broadband echolocation signals
    from prey of Blainville's beaked whales: Modeling and analysis. Master of Science,
    Massachusetts Institute of Technology. <https://doi.org/10.1575/1912/1283>

    Jones, B. A., Lavery, A. C., & Stanton, T. K. (2009). Use of the distorted
    wave Born approximation to predict scattering by inhomogeneous objects:
    Application to squid. The Journal of the Acoustical Society of America,
    125(1), 73-88. <https://doi.org/10.1121/1.3021298>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    # Make sure things are numpy arrays
    rho = np.atleast_1d(rho)
    c = np.atleast_1d(c)
    voxel_size = np.array(voxel_size)

    # volume of the voxels [m^3]
    dv = voxel_size.prod()

    # input parameter checks
    if not len(volume.shape) == 3:
        raise TypeError('The volume input variable must be 3-dimensional.')

    if not voxel_size.shape[0] == 3:
        raise TypeError('The voxel_size input variable must contain 3 items.')

    if not np.any(voxel_size > 0):
        raise ValueError('All voxel_size values must be positive.')

    if f < 0.0:
        raise ValueError('The f input variable must contain only positive values.')

    if (theta < -0.0) or (theta > 180.0):
        raise ValueError('The theta (pitch) angle must be between -180.0 and +180.0')

    if (phi < -180.0) or (phi > 180.0):
        raise ValueError('The phi (roll) angle must be between -180.0 and +180.0')

    if volume.min() != 0:
        raise ValueError('The volume input variable must contain zeros.')

    categories = np.unique(volume)
    if not len(categories == (volume.max() + 1)):
        raise ValueError('The integers in volume must include all values in the series '
                         '(0, 1, 2, ..., n), where n is the largest integer in volume.')

    if not len(rho) >= len(categories):
        raise ValueError('The target_rho variable must contain at least as many values as '
                         'unique integers in the volume variable.')

    if not len(c) >= len(categories):
        raise ValueError('The target_c variable must contain at least as many values '
                         'as unique integers in the volume variable.')

    # density and sound speed ratios for all object materials
    g = rho[1:] / rho[0]
    h = c[1:] / c[0]

    # Do the pitch and roll rotations

    # Convert echoSMs rotation angles (which are intrinsic) into extrinsic as
    # that is what ndimage.rotate() below uses.
    if phi == 0.0:  # short circuit the coordinate transformation if we can
        pitch = theta-90
        roll = 0.0
    else:
        rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
        # for backscatter we don't care about yaw
        _, pitch, roll = rot.as_euler('zyz', degrees=True)

    v = ndimage.rotate(volume, pitch, axes=(0, 2), order=0)
    v = ndimage.rotate(v, roll, axes=(1, 2), order=0)

    categories = np.unique(v)  # or just take the max?

    # wavenumbers in the various media
    k = 2.0*np.pi * f / c

    # DWBA coefficients
    # amplitudes in media 1,2,...,n
    Cb = 1.0/(g * h**2) + 1.0/g - 2.0  # gamma_kappa - gamma_rho
    Ca = k[0]**2 * Cb / (4.0*np.pi)  # summation coefficient

    # Differential phase for each voxel.
    dph = np.zeros(v.shape)
    masks = []
    for i, category in enumerate(categories):
        masks.append(np.isin(v, category))
        dph[masks[i]] = k[i] * voxel_size[0]
    masks.pop(0)  # don't need to keep the category[0] mask

    # cumulative summation of phase along the z-direction
    phase = dph.cumsum(axis=2) - dph/2.0
    dA = np.zeros(phase.shape, dtype=np.complex128)

    # differential phases for each voxel
    for i, m in enumerate(masks):
        dA[m] = Ca[i] * np.exp(2.0*1j*phase[m]) * dv

    # Convert to TS
    return 20.0 * np.log10(np.abs(dA.sum()))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/ptdwbamodel.py 
23
24
25
26
27
28
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)

echosms.dwbautils

Miscellaneous functions for the DWBA models.

DWBAdata()

Example datasets for the SDWBA and DWBA models.

Source code in src/echosms/dwbautils.py 
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
def __init__(self):
    # Load in the shapes data
    self.file = Path(__file__).parent/Path('resources')/Path('DWBA_shapes.toml')
    with open(self.file, 'rb') as f:
        try:
            shapes = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Put the shapes into a dict of DWBAorganism().
    self.dwba_models = {}
    for s in shapes['shape']:
        mx = max(s['x'])
        s['x'] = np.array(s['x'])-mx  # put head at x=0
        # Estimate rv_tan from a spline through (x,y,z).
        tck, u = splprep([s['x'], s['y'], s['z']])
        rv_tan = np.vstack(splev(u, tck, der=1))
        # Make sure rv_tan holds only unit vectors
        n = np.linalg.norm(np.vstack(rv_tan), axis=0)
        rv_tan = (rv_tan / n).T

        # Convert the x, y, z lists into a 2D array with one row for each (x,y,z) point
        rv_pos = np.vstack((np.array(s['x']), np.array(s['y']), np.array(s['z']))).T

        organism = DWBAorganism(rv_pos, np.array(s['a']), np.array(s['g']), np.array(s['h']),
                                s['name'], s.get('source', ''), s.get('note', ''), rv_tan,
                                s['aphiaid'], s['length'], s['vernacular'],)
        self.dwba_models[s['name']] = organism
as_dict()

DWBA model shapes as a dict.

RETURNS DESCRIPTION
dict

All the DWBA model shapes. The dataset name is the dict key and the value is an instance of DWBAorganism.
Source code in src/echosms/dwbautils.py 
286
287
288
289
290
291
292
293
294
295
296
def as_dict(self) -> dict:
    """DWBA model shapes as a dict.

    Returns
    -------
    :
        All the DWBA model shapes. The dataset name is the dict key and the value is an
        instance of `DWBAorganism`.

    """
    return self.dwba_models
model(name)

DWBA model shape with requested name.

PARAMETER DESCRIPTION
name

The name of a DWBA model shape.

TYPE: str

RETURNS DESCRIPTION
DWBAorganism

An instance of DWBAorganism or None if there is no model with name.
Source code in src/echosms/dwbautils.py 
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
def model(self, name: str) -> DWBAorganism:
    """DWBA model shape with requested name.

    Parameters
    ----------
    name :
        The name of a DWBA model shape.

    Returns
    -------
    :
        An instance of `DWBAorganism` or None if there is no model with `name`.

    """
    try:
        return self.dwba_models[name]
    except KeyError:
        return None
names()

Available DWBA model names.

Source code in src/echosms/dwbautils.py 
282
283
284
def names(self):
    """Available DWBA model names."""
    return [*self.dwba_models]

DWBAorganism(rv_pos, a, g, h, name, source='', note='', rv_tan=None, aphiaid=1, length=0.0, vernacular_name='') dataclass

DWBA shape and property class to represent an organism.

ATTRIBUTE DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the target shape. Each vector should have three values corresponding to the x, y, and z coordinates [m] of the disc centre, as per the echoSMs coordinate system.

TYPE: iterable[ndarray]

a

The radii [m] of each disc.

TYPE: ndarray

g

The density contrast between medium and organism (organism divied by medium).

TYPE: ndarray

h

The sound speed contrast betwwen medium and organism (organism divied by medium).

TYPE: ndarray

source

A link to the source of the data.

TYPE: str

note

Information about the data.

TYPE: str

rv_tan

An iterable of unit vectors of the tangent to the body axis at the points given by (x, y, z). Each vector has three values corresponding to the x, y, and z components of the tangent vector [m]. If not given, unit vectors pointing along the positive x-axis are used.

TYPE: ndarray

plot()

Do a simple plot of the DWBA model data.

Source code in src/echosms/dwbautils.py 
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
def plot(self):
    """Do a simple plot of the DWBA model data."""
    import matplotlib.pyplot as plt
    fig, axs = plt.subplots(2, 1, layout='compressed')
    x = self.rv_pos[:, 0]*1e3
    outline1 = (-self.rv_pos[:, 1] + self.a)*1e3
    outline2 = (-self.rv_pos[:, 1] - self.a)*1e3
    axs[0].plot(x, -self.rv_pos[:, 1]*1e3, '.-', c='C0')
    axs[0].plot(x, outline1, c='C1')
    axs[0].plot(x, outline2, c='C1')
    axs[0].plot([x[0], x[0]], [outline1[0], outline2[0]], c='C1')
    axs[0].plot([x[-1], x[-1]], [outline1[-1], outline2[-1]], c='C1')
    axs[0].set_title('Dorsal', loc='left', fontsize=8)
    axs[0].axis('scaled')
    axs[0].xaxis.set_inverted(True)

    outline1 = (-self.rv_pos[:, 2] + self.a)*1e3
    outline2 = (-self.rv_pos[:, 2] - self.a)*1e3
    axs[1].plot(x, -self.rv_pos[:, 2]*1e3, '.-', c='C0')
    axs[1].plot(x, outline1, c='C1')
    axs[1].plot(x, outline2, c='C1')
    axs[1].plot([x[0], x[0]], [outline1[0], outline2[0]], c='C1')
    axs[1].plot([x[-1], x[-1]], [outline1[-1], outline2[-1]], c='C1')
    axs[1].set_title('Lateral', loc='left', fontsize=8)
    axs[1].axis('scaled')
    axs[1].xaxis.set_inverted(True)

    plt.suptitle(self.name)
    plt.show()

create_dwba_cylinder(radius, length, spacing=0.0001)

Create shape description variables for the DWBA model for cylinders.

The shape descriptions are essentially a set of discs and their orientation.

PARAMETER DESCRIPTION
radius

The radius [m] of the cylinder.

TYPE: float

length

The length [m] of the cylinder.

TYPE: float

spacing

The spacing [m] between successive discs.

TYPE: float DEFAULT: 0.0001

RETURNS DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the cylinder. Each vector has three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: iterable[ndarray]

rv_tan

An iterable of unit vectors of the tangent to the cylinder body axis at the points given in rv_pos. Each vector has three values corresponding to the x, y, and z components of the tangent vector.

TYPE: iterable[ndarray]

a

The radii [m] of the discs that define the spheroid.

TYPE: iterable

Source code in src/echosms/dwbautils.py 
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
def create_dwba_cylinder(radius: float, length: float, spacing: float = 0.0001):
    """Create shape description variables for the DWBA model for cylinders.

    The shape descriptions are essentially a set of discs and their orientation.

    Parameters
    ----------
    radius :
        The radius [m] of the cylinder.
    length :
        The length [m] of the cylinder.
    spacing :
        The spacing [m] between successive discs.

    Returns
    -------
    rv_pos : iterable[np.ndarray]
        An iterable of vectors of the 3D positions of the centre of each disc that
        defines the cylinder. Each vector has three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan : iterable[np.ndarray]
        An iterable of unit vectors of the tangent to the cylinder body axis at
        the points given in `rv_pos`. Each vector has three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    a : iterable
        The radii [m] of the discs that define the spheroid.
    """
    pos = np.linspace(0, length, int(round(length/spacing)))
    rv_pos = [np.array([x, 0, 0]) for x in pos]
    rv_tan = [np.array([1, 0, 0])] * len(pos)
    a = [radius] * len(pos)

    return rv_pos, rv_tan, a

create_dwba_from_xyza(x, y, z, a, name, g=1.0, h=1.0, source='', note='')

Create a DWBAorganism instance from shape data.

Converts a centreline and radius definiton of the DWBA shape into that required by the echoSMs implementation of the DWBA (centreline, tangential, and radii vectors).

PARAMETER DESCRIPTION
x

x-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

y

y-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

z

z-coordinates [m] of the centreline of the DWBA shape as per the echoSMs coordinate system.

TYPE: Iterable[float]

a

radius [m] of the DWBA shape at each centreline (x,y,z) position.

TYPE: Iterable[float]

name

A name for the organism.

TYPE: str

source

A link/URL/DOI or description for the source of the data.

TYPE: str DEFAULT: ''

note

Notes about the organism or data.

TYPE: str DEFAULT: ''

g

A single value of g, the ratio of organism density divided by the medium density. This is applied to all parts of the shape.

TYPE: float DEFAULT: 1.0

h

A single value of h, the ratio of organism sound speed divided by the medium sound speed. This is applied to all parts of the shape.

TYPE: float DEFAULT: 1.0

RETURNS DESCRIPTION
An instance of DWBAorganism.
Notes

Here is an example of how to use this function to read in .sat files from the ZooScatR package and convert them into the format required by the echoSMs DWBA model.

import pandas as pd
from echosms import create_dwba_from_xyza

filepath = 'shape.sat'

s = pd.read_csv(filepath, delimiter=' ', names=['y', 'x', 'a'])

# Adjust the data to match the echoSMs units
s /= 1000  # convert mm to m

# .sat files don't have a z coordinate but echoSMs requires it
s['z'] = 0.0

# An example of flipping left to right to match the echoSMS coordinate system
s['x'] = max(s['x']) - s['x']

shape = create_dwba_from_xyza(s['x'], s['y'], s['z'], s['a'], name=filepath, g=1.05, h=1.05)
Source code in src/echosms/dwbautils.py 
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
def create_dwba_from_xyza(x, y, z, a, name: str, g: float = 1.0, h: float = 1.0,
                          source: str = '', note: str = ''):
    """Create a DWBAorganism instance from shape data.

    Converts a centreline and radius definiton of the DWBA shape into
    that required by the echoSMs implementation of the DWBA (centreline, tangential, and
    radii vectors).

    Parameters
    ----------
    x : Iterable[float]
        x-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems).
    y : Iterable[float]
        y-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems).
    z : Iterable[float]
        z-coordinates [m] of the centreline of the DWBA shape as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems).
    a : Iterable[float]
        radius [m] of the DWBA shape at each centreline (x,y,z) position.
    name :
        A name for the organism.
    source :
        A link/URL/DOI or description for the source of the data.
    note :
        Notes about the organism or data.
    g :
        A single value of g, the ratio of organism density divided by the
        medium density. This is applied to all parts of the shape.
    h :
        A single value of h, the ratio of organism sound speed divided by
        the medium sound speed. This is applied to all parts of the shape.

    Returns
    -------
        An instance of DWBAorganism.

    Notes
    -----
    Here is an example of how to use this function to read in .sat files
    from the ZooScatR package and convert them into the format required
    by the echoSMs DWBA model.

    ```py
    import pandas as pd
    from echosms import create_dwba_from_xyza

    filepath = 'shape.sat'

    s = pd.read_csv(filepath, delimiter=' ', names=['y', 'x', 'a'])

    # Adjust the data to match the echoSMs units
    s /= 1000  # convert mm to m

    # .sat files don't have a z coordinate but echoSMs requires it
    s['z'] = 0.0

    # An example of flipping left to right to match the echoSMS coordinate system
    s['x'] = max(s['x']) - s['x']

    shape = create_dwba_from_xyza(s['x'], s['y'], s['z'], s['a'], name=filepath, g=1.05, h=1.05)

    ```
    """
    # Estimate rv_tan from a spline through (x,y).
    tck, u = splprep([x, y, z])
    rv_tan = np.vstack(splev(u, tck, der=1))
    # Make sure rv_tan holds only unit vectors
    n = np.linalg.norm(np.vstack(rv_tan), axis=0)
    rv_tan = (rv_tan / n).T

    # Convert the x, y, and z into a 2D array and get one row for each (x,y,z) point.
    rv_pos = np.vstack((x, y, z)).T

    gg = np.full((1, rv_pos.shape[0]), g)
    hh = np.full((1, rv_pos.shape[0]), h)

    return DWBAorganism(rv_pos, a, gg, hh, name, source, note, rv_tan)

create_dwba_spheroid(major_radius, minor_radius, spacing=0.0001)

Create shape description variables for the DWBA model for spheres and prolate spheroids.

The shape descriptions are essentially a set of discs and their orientation.

Notes

Currently only supports prolate spheroids and spheres (set major_radius and minor_radius to the same value to get a sphere).

PARAMETER DESCRIPTION
major_radius

The major radius [m] of the spheroid.

TYPE: float

minor_radius

The minor radius [m] of the spheroid.

TYPE: float

spacing

The spacing [m] between successive discs.

TYPE: float DEFAULT: 0.0001

RETURNS DESCRIPTION
rv_pos

An iterable of vectors of the 3D positions of the centre of each disc that defines the spheroid. Each vector has three values corresponding to the x, y, and z coordinates [m] of the disc centre.

TYPE: iterable[ndarray]

rv_tan

An iterable of unit vectors of the tangent to the target body axis at the points given in rv_pos. Each vector has three values corresponding to the x, y, and z components of the tangent vector.

TYPE: iterable[ndarray]

a

The radii [m] of the discs that define the spheroid.

TYPE: iterable

Source code in src/echosms/dwbautils.py 
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def create_dwba_spheroid(major_radius: float, minor_radius: float, spacing: float = 0.0001):
    """Create shape description variables for the DWBA model for spheres and prolate spheroids.

    The shape descriptions are essentially a set of discs and their orientation.

    Notes
    -----
    Currently only supports prolate spheroids and spheres (set `major_radius` and `minor_radius`
    to the same value to get a sphere).

    Parameters
    ----------
    major_radius :
        The major radius [m] of the spheroid.
    minor_radius :
        The minor radius [m] of the spheroid.
    spacing :
        The spacing [m] between successive discs.

    Returns
    -------
    rv_pos : iterable[np.ndarray]
        An iterable of vectors of the 3D positions of the centre of each disc that
        defines the spheroid. Each vector has three values corresponding to
        the _x_, _y_, and _z_ coordinates [m] of the disc centre.
    rv_tan : iterable[np.ndarray]
        An iterable of unit vectors of the tangent to the target body axis at
        the points given in `rv_pos`. Each vector has three values corresponding to
        the _x_, _y_, and _z_ components of the tangent vector.
    a : iterable
        The radii [m] of the discs that define the spheroid.
    """
    v = np.linspace(0, np.pi, int(round(2*major_radius/spacing)))
    a = minor_radius*np.sin(v)  # radius at points along the spheroid
    x = major_radius-major_radius*np.cos(v)  # shift so that origin is at one end

    # List of disc position vectors
    rv_pos = [np.array([i, 0, 0]) for i in x]

    # List of tangent vectors to sphere axis (all the same for spheroids)
    rv_tan = [np.array([1, 0, 0])] * len(x)

    return rv_pos, rv_tan, a

echosms.ESModel()

Bases: ScatterModelBase

Elastic sphere (ES) scattering model.

This class calculates acoustic backscatter from elastic spheres.

Source code in src/echosms/esmodel.py 
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'elastic sphere'
    self.short_name = 'es'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.elastic]
    self.shapes = ['sphere']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, f, target_longitudinal_c, target_transverse_c, target_rho, validate_parameters=True, **kwargs)

Calculate the backscatter from an elastic sphere for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the sphere [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the sphere [kg/m³].

TYPE: float

a

Radius of the sphere [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

target_longitudinal_c

Longitudinal sound speed in the material inside the sphere [m/s].

TYPE: float

target_transverse_c

Transverse sound speed in the material inside the sphere [m/s].

TYPE: float

target_rho

Density of the material inside the sphere [kg/m³].

TYPE: float

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the sphere [dB].
Notes

The class implements the code in MacLennan (1981).

References

MacLennan, D. N. (1981). The Theory of Solid Spheres as Sonar Calibration Targets. Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries for Scotland.

Source code in src/echosms/esmodel.py 
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def calculate_ts_single(self, medium_c, medium_rho, a, f,
                        target_longitudinal_c, target_transverse_c, target_rho,
                        validate_parameters=True,
                        **kwargs) -> float:
    """
    Calculate the backscatter from an elastic sphere for one set of parameters.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the sphere [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the sphere [kg/m³].
    a : float
        Radius of the sphere [m].
    f : float
        Frequency to calculate the scattering at [Hz].
    target_longitudinal_c : float
        Longitudinal sound speed in the material inside the sphere [m/s].
    target_transverse_c : float
        Transverse sound speed in the material inside the sphere [m/s].
    target_rho : float
        Density of the material inside the sphere [kg/m³].
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the sphere [dB].

    Notes
    -----
    The class implements the code in MacLennan (1981).

    References
    ----------
    MacLennan, D. N. (1981). The Theory of Solid Spheres as Sonar Calibration Targets.
    Scottish Fisheries Research Report Number 22. Department of Agriculture and Fisheries
    for Scotland.
    """
    if validate_parameters:
        self.validate_parameters(locals())

    q = wavenumber(medium_c, f)*a
    q1 = q*medium_c/target_longitudinal_c
    q2 = q*medium_c/target_transverse_c
    alpha = 2. * (target_rho/medium_rho) * (target_transverse_c/medium_c)**2
    beta = (target_rho/medium_rho) * (target_longitudinal_c/medium_c)**2 - alpha

    # Use n instead of l (ell) because l looks like 1.
    def S(n):
        A2 = (n**2 + n-2) * spherical_jn(n, q2) + q2**2 * spherical_jnpp(n, q2)
        A1 = 2*n*(n+1) * (q1*spherical_jn(n, q1, True) - spherical_jn(n, q1))
        B2 = A2*q1**2 * (beta*spherical_jn(n, q1) - alpha*spherical_jnpp(n, q1))\
            - A1*alpha * (spherical_jn(n, q2) - q2*spherical_jn(n, q2, True))
        B1 = q * (A2*q1*spherical_jn(n, q1, True) - A1*spherical_jn(n, q2))
        eta_n = atan(-(B2*spherical_jn(n, q, True) - B1*spherical_jn(n, q))
                     / (B2*spherical_yn(n, q, True) - B1*spherical_yn(n, q)))

        return (-1)**n * (2*n+1) * sin(eta_n) * exp(1j*eta_n)

    # Estimate the number of terms to use in the summation
    n_max = round(q+10)
    tol = 1e-10  # somewhat arbitrary
    while abs(S(n_max)) > tol:
        n_max += 10

    if n_max > 200:
        warn('TS results may be inaccurate because the modal series required a large '
             f'number ({n_max}) of terms to converge.')

    n = range(n_max)

    f_inf = -2.0/q * sum(map(S, n))

    return 10*log10(a**2 * abs(f_inf)**2 / 4.0)

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/esmodel.py 
26
27
28
29
30
31
32
33
34
35
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'medium_c', 'a', 'f',
                                      'target_longitudinal_c',
                                      'target_transverse_c', 'target_rho'])

echosms.HPModel()

Bases: ScatterModelBase

High-pass (HP) scattering model.

Source code in src/echosms/hpmodel.py 
11
12
13
14
15
16
17
18
def __init__(self):
    super().__init__()
    self.long_name = 'high pass'
    self.short_name = 'hp'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled, bt.elastic, bt.fixed_rigid]
    self.shapes = ['sphere', 'prolate spheroid', 'cylinder', 'bent cylinder']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(shape, medium_c, a, f, boundary_type, medium_rho=None, target_c=None, target_rho=None, theta=None, L=None, rho_c=None, irregular=False, validate_parameters=True, **kwargs)

Calculate the backscatter using the high pass model for one set of parameters.

PARAMETER DESCRIPTION
shape

The shape to model. Must be one of the shapes given in the shapes variable.

TYPE: str

medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³]. Not required when boundary_type is fixed_rigid.

TYPE: float DEFAULT: None

target_c

Longitudinal sound speed in the material inside the target [m/s]. Not required when boundary_type is fixed_rigid.

TYPE: float DEFAULT: None

target_rho

Density of the material inside the target [kg/m³]. Not required when boundary_type is fixed_rigid.

TYPE: float DEFAULT: None

a

Radius of the sphere, length of semi-minor axis of the prolate spheriod, or cylindrical radius of the straight or bent cylinder [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The boundary type for the model.

TYPE: boundary_type

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°]. Only required for the straight cylinder shape.

TYPE: float DEFAULT: None

L

Total length of the prolate spheroid and straight cylinder, or arc length of the bent cylinder [m]. Only required for prolate spheroid, cylinder, and bent cylinder shapes.

TYPE: float DEFAULT: None

rho_c

Radius of curvature of the axis of the bent cylinder [m]. Only required for the bent cylinder shape.

TYPE: float DEFAULT: None

irregular

Set to True if the modelled object is not exactly a sphere, prolate spheroid, straight or uniformly beny cylinder.

DEFAULT: False

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the sphere [dB].
Notes

The class implements the high-pass model in Stanton (1989) for spheres, prolate spheroids, cylinders, and bent cylinders with fluid filled, elastic, and rigid fixed boundary conditions. There are several restrictions on valid input parameters, so a careful reading of Stanton (1989) is recommended.

The theta angle convention used in Stanton (1989) is the same as the echoSMs coordinate system convention.

Stanton (1989) also provides parameters for gas-filled shapes, but more prior knowledge is required about the gas for useful results (e.g., damping characteristics of the gas and medium) and these have not been implemented here. There are many other models available that accurately model gas-filled shapes such that the lack of the high-pass gas-filled model should not be missed.

References

Stanton, T. K. (1989). Simple approximate formulas for backscattering of sound by spherical and elongated objects. The Journal of the Acoustical Society of America, 86(4), 1499-1510. https://doi.org/10.1121/1.398711

Source code in src/echosms/hpmodel.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
def calculate_ts_single(self, shape, medium_c, a, f, boundary_type: bt, medium_rho=None,
                        target_c=None, target_rho=None,
                        theta=None,
                        L=None, rho_c=None,
                        irregular=False,
                        validate_parameters=True, **kwargs) -> float:
    """
    Calculate the backscatter using the high pass model for one set of parameters.

    Parameters
    ----------
    shape : str
        The shape to model. Must be one of the shapes given in the `shapes` variable.
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³]. Not required when
        `boundary_type` is `fixed_rigid`.
    target_c : float
        Longitudinal sound speed in the material inside the target [m/s]. Not required when
        `boundary_type` is `fixed_rigid`.
    target_rho : float
        Density of the material inside the target [kg/m³]. Not required when
        `boundary_type` is `fixed_rigid`.
    a : float
        Radius of the sphere, length of semi-minor axis of the prolate spheriod, or cylindrical
        radius of the straight or bent cylinder [m].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The boundary type for the model.
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°]. Only required for the straight cylinder shape.
    L : float
        Total length of the prolate spheroid and straight cylinder, or arc length of
        the bent cylinder [m]. Only required for prolate spheroid, cylinder, and bent cylinder
        shapes.
    rho_c : float
        Radius of curvature of the axis of the bent cylinder [m]. Only required for the
        bent cylinder shape.
    irregular :
        Set to `True` if the modelled object is not exactly a sphere, prolate spheroid,
        straight or uniformly beny cylinder.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the sphere [dB].

    Notes
    -----
    The class implements the high-pass model in Stanton (1989) for spheres, prolate spheroids,
    cylinders, and bent cylinders with fluid filled, elastic, and rigid fixed boundary
    conditions. There are several restrictions on valid input parameters, so a careful
    reading of Stanton (1989) is recommended.

    The theta angle convention used in Stanton (1989) is the same as the echoSMs
    [coordinate system convention](https://ices-tools-dev.github.io/echoSMs/
    conventions/#coordinate-systems).

    Stanton (1989) also provides parameters for gas-filled shapes, but more
    prior knowledge is required about the gas for useful results (e.g., damping
    characteristics of the gas and medium) and these have not been implemented here. There
    are many other models available that accurately model gas-filled shapes such that the
    lack of the high-pass gas-filled model should not be missed.

    References
    ----------
    Stanton, T. K. (1989). Simple approximate formulas for backscattering of sound
    by spherical and elongated objects. The Journal of the Acoustical Society of
    America, 86(4), 1499-1510.
    <https://doi.org/10.1121/1.398711>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type == bt.fixed_rigid:
        # just need something large
        g = 1e20
        h = 1e20
    else:
        g = target_rho/medium_rho
        h = target_c/medium_c

    k = wavenumber(medium_c, f)

    G = 1.0
    F = 1.0

    def alpha_pic(g, h):
        return (1-g*h*h)/(2*g*h*h) + (1-g)/(1+g)

    match shape:
        case 'sphere':
            alpha_pis = (1-g*h*h)/(3*g*h*h) + (1-g)/(1+2*g)
            R = (g*h-1)/(g*h+1)
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 40 * (k*a)**(-0.4)
                        G = 1-0.8*exp(-2.5*(k*a-2.25)**2)
                    case bt.elastic:
                        F = 15 * (k*a)**(-1.9)
                    case bt.fixed_rigid:
                        F = 15 * (k*a)**(-1.9)

            sigma_bs = a*a * (k*a)**4 * alpha_pis**2 * G\
                / (1 + 4*(k*a)**4 * alpha_pis**2/(R**2 * F))
        case 'prolate spheroid':
            a_pic = alpha_pic(g, h)
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 2.5 * (k*a)**(1.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.3)**2)
                    case bt.elastic:
                        F = 1.8 * (k*a)**(-0.4)
                    case bt.fixed_rigid:
                        F = 1.8 * (k*a)**(-0.4)

            sigma_bs = 1/9 * L*L * (k*a)**4 * a_pic**2 * G\
                / (1 + 16/9*(k*a)**4 * a_pic**2/(R**2 * F))
        case 'cylinder':
            theta = radians(theta)
            a_pic = alpha_pic(g, h)
            s = sin(k*L*cos(theta)) / (k*L*cos(theta))
            Ka = k*sin(theta)*a
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 3 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case bt.elastic:
                        F = 3.5 * (k*a)**(-1.0)
                    case bt.fixed_rigid:
                        F = 3.5 * (k*a)**(-1.0)

            sigma_bs = 0.25 * L*L * (Ka)**4 * a_pic**2 * s*s * G\
                / (1 + pi*(Ka)**4 * a_pic**2/(R**2 * F))
        case 'bent cylinder':
            a_pic = alpha_pic(g, h)
            H = 1.
            if irregular:
                match boundary_type:
                    case bt.fluid_filled:
                        F = 3.0 * (k*a)**(0.65)
                        G = 1-0.8*exp(-2.5*(k*a-2.0)**2)
                    case bt.elastic:
                        F = 2.5 * (k*a)**(-1.0)
                    case bt.fixed_rigid:
                        F = 2.5 * (k*a)**(-1.0)

            sigma_bs = 0.25 * L*L * (k*a)**4 * a_pic**2 * H*H*G\
                / (1 + L*L*(k*a)**4 * a_pic**2 * H*H/(rho_c*a*R**2 * F))

    return 10*log10(sigma_bs)

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/hpmodel.py 
20
21
22
23
24
25
26
27
28
29
30
31
32
33
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_c', 'a', 'f'])

    if not p['shape'].isin(self.shapes).all():
        raise ValueError('The shape parameter must be one of: ' + ', '.join(self.shapes))

    if not p['boundary_type'].isin(self.boundary_types).all():
        raise ValueError('The boundary_type parameter must be one of: ' +
                         ', '.join(self.boundary_types))

echosms.KAModel()

Bases: ScatterModelBase

Kirchhoff approximation (KA) scattering model.

This class calculates acoustic scatter from arbitrary surfaces.

Source code in src/echosms/kamodel.py 
16
17
18
19
20
21
22
23
24
def __init__(self):
    super().__init__()
    self.long_name = 'Kirchhoff approximation'
    self.short_name = 'ka'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.pressure_release]
    self.shapes = ['closed surfaces']
    self.max_ka = 20  # [1]
    self.no_expand_parameters = ['mesh']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, theta, phi, f, mesh, boundary_type, validate_parameters=True, **kwargs)

Calculate the scatter using the ka model for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

phi

Roll angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

mesh

The triangular mesh that defines the scattering surface. This parameter must provide attributes with names of:

  • triangles_center (the position of the centre of each triangular face [m]),
  • face_normals (the outward-pointing unit normals for each triangular face),
  • area_faces (the area of each triangular face [m²]).

A suitable library for creating and manipulating triangular meshes is trimesh. Trimesh will accept the usual nodes/facets definition of a mesh and calculate the above attributes automatically.

TYPE: Any

boundary_type

The boundary type. Supported types are given in the boundary_types class variable.

TYPE: boundary_type

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Foote (1985).

References

Foote, K. G. (1985). Rather-high-frequency sound scattering of swimbladdered fish. The Journal of the Acoustical Society of America, 78(2), 688–700. https://doi.org/10.1121/1.392438

Source code in src/echosms/kamodel.py 
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
def calculate_ts_single(self, medium_c, theta, phi, f, mesh,
                        boundary_type: bt, validate_parameters=True, **kwargs) -> float:
    """
    Calculate the scatter using the ka model for one set of parameters.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    phi : float
        Roll angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    mesh : Any
        The triangular mesh that defines the scattering surface. This parameter must provide
        attributes with names of:

        - `triangles_center` (the position of the centre of each triangular face [m]),
        - `face_normals` (the outward-pointing unit normals for each triangular face),
        - `area_faces` (the area of each triangular face [m²]).

        A suitable library for creating and manipulating triangular meshes
        is [trimesh](https://trimesh.org). Trimesh will accept the usual nodes/facets
        definition of a mesh and calculate the above attributes automatically.
    boundary_type :
        The boundary type. Supported types are given in the `boundary_types` class variable.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Foote (1985).

    References
    ----------
    Foote, K. G. (1985). Rather-high-frequency sound scattering of swimbladdered fish.
    The Journal of the Acoustical Society of America, 78(2), 688–700.
    <https://doi.org/10.1121/1.392438>

    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type not in self.boundary_types:
        raise ValueError(f'The {self.long_name} model does not support '
                         f'a model type of "{boundary_type}".')

    # This model keeps the organism fixed and varies the incident wave vector. So need
    # to convert the theta and phi echoSMs coordinate sytem Tait-Bryan angles
    # into an (x,y,z) vector.

    # Acoustic wave incident vector and its' norm
    rot = R.from_euler('ZYX', (0, theta-90, -phi), degrees=True)
    k_norm = rot.as_matrix() @ np.array([[0, 0, 1]]).T
    k = k_norm * wavenumber(medium_c, f)

    r = mesh.triangles_center  # position vector of each surface element
    dS = mesh.area_faces.reshape((-1, 1))  # [m^2]

    kn_nn = mesh.face_normals @ k_norm

    fbs = 1./wavelength(medium_c, f)\
        * np.sum(np.exp(2j*r @ k) * np.heaviside(kn_nn, 0.5) * kn_nn * dS)

    return 10*log10(abs(fbs)**2)  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/kamodel.py 
26
27
28
29
30
31
32
33
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_c', 'f'])

KRM model & utilities

echosms.KRMModel()

Bases: ScatterModelBase

Kirchhoff ray mode (KRM) scattering model.

Source code in src/echosms/krmmodel.py 
30
31
32
33
34
35
36
37
38
def __init__(self):
    super().__init__()
    self.long_name = 'Kirchhoff ray mode'
    self.short_name = 'krm'
    self.analytical_type = 'approximate'
    self.boundary_types = [bt.fluid_filled]
    self.shapes = ['closed surfaces']
    self.max_ka = 20  # [1]
    self.no_expand_parameters = ['bodies']

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, theta, f, organism, high_ka_medium='body', low_ka_medium='body', validate_parameters=True, **kwargs)

Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.

Warning

The mode solution (low ka) part of this model has not yet been verified to give correct results.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density in the fluid medium surrounding the target [kg/m³]

TYPE: float

theta

Pitch angle to calculate the scattering at, as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

organism

The shapes that make up the model. This is typically a shape for the body and zero or more enclosed shapes that repesent internal parts of the organism.
high_ka_medium

If set to body the sound speed and density of the organism body is used for the fluid surrounding any inclusions. If set to anything else (e.g., water) the sound speed and density given by medium_c and medium_rho are used. This parameter applies to the Kirchhoff approximation part of the model (i.e., high ka) and corresponds to the use (or not) of the approximation given in Clay & Horne (1994) on the line immediately below Eqn (13): k_b ≈ k at low contrast.

DEFAULT: 'body'

low_ka_medium

If set to body the sound speed and density of the organism body is used for the fluid surrounding any inclusions. If set to anything else (e.g., water) the sound speed and density given by medium_c and medium_rho are used. This parameter applies to the mode solution part of the model (i.e., low ka) and corresponds to the use (or not) of the approximation given in Clay & Horne (1994) on the line immediately below Eqn (13): k_b ≈ k at low contrast.

DEFAULT: 'body'

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Clay & Horne (1994) and when ka < 0.15 uses Clay (1992).

The high_ka_medium and low_ka_medium parameters allow the user to select which medium surrounds the inclusions (e.g., the swimbladder) - the fish body or the water surrounding the fish body. The equations in Clay & Horne (1994) used the body but included a sentence saying that the water could be used at low contrast (between the water and the body). A later paper (Horne & Jech, 1999) used water for low ka (the mode solution) and the body for higher ka (the Kirchhoff approximation). Some open-source KRM model codes always use the water.

References

Clay, C. S. (1992). Composite ray-mode approximations for backscattered sound from gas-filled cylinders and swimbladders. The Journal of the Acoustical Society of America, 92(4), 2173–2180. https://doi.org/10.1121/1.405211

Clay, C. S., & Horne, J. K. (1994). Acoustic models of fish: The Atlantic cod (Gadus morhua). The Journal of the Acoustical Society of America, 96(3), 1661–1668. https://doi.org/10.1121/1.410245

Horne, J. K., & J. M. Jech. (1999). Multi–frequency estimates of fish abundance: constraints of rather high frequencies. ICES Journal of Marine Science, 56 (2), 184–199. https://doi.org/10.1006/jmsc.1998.0432

Source code in src/echosms/krmmodel.py 
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
def calculate_ts_single(self, medium_c, medium_rho, theta, f, organism,
                        high_ka_medium='body', low_ka_medium='body',
                        validate_parameters=True, **kwargs) -> float:
    """
    Calculate the scatter using the Kirchhoff ray mode model for one set of parameters.

    Warning
    --------
    The mode solution (low _ka_) part of this model has not yet been verified to give
    correct results.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density in the fluid medium surrounding the target [kg/m³]
    theta : float
        Pitch angle to calculate the scattering at, as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    organism: KRMorganism
        The shapes that make up the model. This is typically a shape for the body and zero or
        more enclosed shapes that repesent internal parts of the organism.
    high_ka_medium:
        If set to `body` the sound speed and density of the organism body is used for
        the fluid surrounding any inclusions. If set to anything else (e.g., `water`)
        the sound speed and density given by `medium_c` and `medium_rho` are used.
        This parameter applies to the Kirchhoff approximation part of
        the model (i.e., high _ka_) and corresponds to the use (or not) of the
        approximation given in Clay & Horne (1994) on the line immediately below Eqn (13):
        _k_b ≈ k at low contrast_.
    low_ka_medium:
        If set to `body` the sound speed and density of the organism body is used for
        the fluid surrounding any inclusions. If set to anything else (e.g., `water`)
        the sound speed and density given by `medium_c` and `medium_rho` are used.
        This parameter applies to the mode solution part of the model (i.e., low _ka_)
        and corresponds to the use (or not) of the approximation given in Clay & Horne (1994)
        on the line immediately below Eqn (13): _k_b ≈ k at low contrast_.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Clay & Horne (1994) and when _ka_ < 0.15 uses Clay (1992).

    The `high_ka_medium` and `low_ka_medium` parameters allow the user to select which
    medium surrounds the inclusions (e.g., the swimbladder) - the fish body or
    the water surrounding the fish body. The equations in Clay & Horne (1994) used the body
    but included a sentence saying that the water could be used at low contrast
    (between the water and the body). A later paper
    (Horne & Jech, 1999) used water for low _ka_ (the mode solution) and the body for
    higher _ka_ (the Kirchhoff approximation). Some open-source KRM model codes always
    use the water.

    References
    ----------
    Clay, C. S. (1992). Composite ray-mode approximations for backscattered sound from
    gas-filled cylinders and swimbladders. The Journal of the Acoustical Society of
    America, 92(4), 2173–2180.
    <https://doi.org/10.1121/1.405211>

    Clay, C. S., & Horne, J. K. (1994). Acoustic models of fish: The Atlantic cod
    (_Gadus morhua_). The Journal of the Acoustical Society of America, 96(3), 1661–1668.
    <https://doi.org/10.1121/1.410245>

    Horne, J. K., & J. M. Jech. (1999). Multi–frequency estimates of fish abundance:
    constraints of rather high frequencies. ICES Journal of Marine Science, 56 (2), 184–199.
    <https://doi.org/10.1006/jmsc.1998.0432>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    theta = radians(theta)

    body = organism.body

    k = wavenumber(medium_c, f)
    k_b = wavenumber(body.c, f)

    # Reflection coefficient between water and body
    R_wb = (body.rho*body.c - medium_rho*medium_c)\
        / (body.rho*body.c + medium_rho*medium_c)
    TwbTbw = 1-R_wb**2  # Eqn (15)

    sl = []  # scattering lengths for inclusions
    for incl in organism.inclusions:
        # Reflection coefficient between body and inclusion
        # The paper gives R_bc in terms of g & h, but it can also be done in the
        # same manner as R_wb above.
        gp = incl.rho / body.rho  # p is 'prime' to fit with paper notation
        hp = incl.c / body.c

        R_bc = (gp*hp-1) / (gp*hp+1)  # Eqn (9)

        # Equivalent radius of inclusion (as per Part A of paper)
        a_e = sqrt(incl.volume() / (pi * incl.length()))

        # Choose which modelling approach to use
        if k*a_e < 0.15:  # Do the mode solution for the inclusion
            if low_ka_medium != 'body':
                gp = incl.rho / medium_rho
                hp = incl.c / medium_c
            sl.append(self._mode_solution(1/gp, 1/hp, k, a_e, incl.length(), theta))
        elif incl.boundary == bt.pressure_release:
            kk = k_b if high_ka_medium == 'body' else k
            sl.append(self._soft_KA(incl, k, kk, R_bc, TwbTbw, theta))
        elif incl.boundary == bt.fluid_filled:
            kk = k_b if high_ka_medium == 'body' else k
            sl.append(self._fluid_KA(incl, k, kk, R_bc, TwbTbw, theta))
        else:
            raise ValueError(f'Unsupported boundary of "{incl.boundary}" for KRM inclusion')

    # Do the Kirchhoff-ray approximation for the body. This is always done as a fluid.
    body_sl = self._fluid_KA(body, k, k_b, R_wb, TwbTbw, theta)

    return 20*log10(abs(body_sl + sum(sl)))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/krmmodel.py 
40
41
42
43
44
45
46
47
48
49
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_positive(p, ['medium_c', 'f'])

    if np.any(np.atleast_1d(p['theta']) < 65) or np.any(np.atleast_1d(p['theta']) > 115):
        raise KeyError('Incidence angle(s) (theta) are outside 65 to 115°')

echosms.KRMdata()

Example datasets for the KRM model.

Source code in src/echosms/krmdata.py 
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
def __init__(self):
    # Load in the NOAA KRM shapes data
    self.file = Path(__file__).parent/Path('resources')/Path('KRM_shapes.toml')
    with open(self.file, 'rb') as f:
        try:
            shapes = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Put the shapes into a dict of KRMorganism(). Use some default values for sound speed and
    # density
    self.krm_models = {}
    for s in shapes['shape']:
        # These KRM data have the head pointing in the -ve x direction,
        # opposite to the echoSMs coordinate convetion, so fix the
        # x-coordinates here when ingesting the data. And set the posterior end
        # of the organism to have x=0
        m = max(s['x_b'])
        body = KRMshape(bt.fluid_filled, -np.array(s['x_b']), np.array(s['w_b']),
                        np.array(s['z_bU']), np.array(s['z_bL']),
                        s['body_c'], s['body_rho'])
        swimbladder = KRMshape(bt.pressure_release, -np.array(s['x_sb']), np.array(s['w_sb']),
                               np.array(s['z_sbU']), np.array(s['z_sbL']),
                               s['swimbladder_c'], s['swimbladder_rho'])
        self.krm_models[s['name']] = KRMorganism(s['name'], s['source'],
                                                 body, [swimbladder],
                                                 s['aphiaid'], s['length'],
                                                 s['vernacular'])

as_dict()

KRM model shapes as a dict.

RETURNS DESCRIPTION
dict

All the KRM model shapes. The dataset name is the dict key and the value is an instance of KRMorganism.
Source code in src/echosms/krmdata.py 
156
157
158
159
160
161
162
163
164
165
166
def as_dict(self) -> dict:
    """KRM model shapes as a dict.

    Returns
    -------
    :
        All the KRM model shapes. The dataset name is the dict key and the value is an instance
        of `KRMorganism`.

    """
    return self.krm_models

model(name)

KRM model shape with requested name.

PARAMETER DESCRIPTION
name

The name of a KRM model shape.

TYPE: str

RETURNS DESCRIPTION
KRMorganism

An instance of KRMorganism or None if there is no model with name.
Source code in src/echosms/krmdata.py 
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
def model(self, name: str) -> KRMorganism:
    """KRM model shape with requested name.

    Parameters
    ----------
    name :
        The name of a KRM model shape.

    Returns
    -------
    :
        An instance of `KRMorganism` or None if there is no model with `name`.

    """
    try:
        return self.krm_models[name]
    except KeyError:
        return None

names()

Available KRM model names.

Source code in src/echosms/krmdata.py 
152
153
154
def names(self):
    """Available KRM model names."""
    return [*self.krm_models]

ts(name) staticmethod

KRM model TS from model name.

PARAMETER DESCRIPTION
name

The name of a KRM model shape.

TYPE: str

RETURNS DESCRIPTION
ndarray

The TS (re 1 m²) for some default model parameters [dB] or None if no TS data are available.
Source code in src/echosms/krmdata.py 
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
@staticmethod
def ts(name: str) -> np.ndarray:
    """KRM model TS from model `name`.

    Parameters
    ----------
    name :
        The name of a KRM model shape.

    Returns
    -------
    :
        The TS (re 1 m²) for some default model parameters [dB] or None if no TS data
        are available.

    """
    # Sometimes there will be TS results for the model (available for testing of the
    # model), so load them in if present.
    tsfile = Path(__file__).parent/Path('resources')/Path('NOAA_KRM_ts_' + name + '.csv')

    if tsfile.exists():
        return pd.read_csv(tsfile)

    return None

echosms.KRMorganism(name, source, body, inclusions, aphiaid=1, length=0.0, vernacular_name='') dataclass

KRM body and inclusion shape(s).

ATTRIBUTE DESCRIPTION
name

A name for the organism.

TYPE: str

source

A link to or description of the source of the organism data.

TYPE: str

body

The shape that represents the organism's body.

TYPE: KRMshape

inclusions

The shapes that are internal to the organism (e.g., swimbladder, backbone, etc)

TYPE: List[KRMshape]

aphiaid

The aphiaID of the organism

TYPE: int

length

The length of the organism (m)

TYPE: float

vernacular_name

A vernacular name of the organism

TYPE: str

plot()

Plot of organism shape.

Source code in src/echosms/krmdata.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
def plot(self):
    """Plot of organism shape."""
    import matplotlib.pyplot as plt

    plt.plot(self.body.x*1e3, self.body.z_U*1e3, self.body.x*1e3, self.body.z_L*1e3, c='black')
    for i in [0, -1]:  # close the ends of the shape
        plt.plot([self.body.x[i]*1e3]*2, [self.body.z_U[i]*1e3, self.body.z_L[i]*1e3],
                 c='black')

    for s in self.inclusions:
        c = 'C0' if s.boundary == bt.fluid_filled else 'C1'
        plt.plot(s.x*1e3, s.z_U*1e3, s.x*1e3, s.z_L*1e3, c=c)
        for i in [0, -1]:  # close the ends of the shape
            plt.plot([s.x[i]*1e3]*2, [s.z_U[i]*1e3, s.z_L[i]*1e3], c=c)

    plt.gca().set_aspect('equal')
    plt.gca().xaxis.set_inverted(True)
    plt.title(self.name)
    plt.show()

echosms.KRMshape(boundary, x, w, z_U, z_L, c, rho) dataclass

KRM shape and property class.

ATTRIBUTE DESCRIPTION
boundary

The shape bounday condition - either pressure_release or fluid_filled.

TYPE: boundary_type

x

The x-axis coordinates [m].

TYPE: ndarray

w

Width of the shape [m].

TYPE: ndarray

z_U

Distance from the axis to the upper surface of the shape [m].

TYPE: ndarray

z_L

Distance from the axis to the lower surface of the shape [m].

TYPE: ndarray

c

Sound speed in the shape [m/s].

TYPE: float

rho

Density of the shape material [kg/m³].

TYPE: float

length()

Length of the shape.

RETURNS DESCRIPTION
float

The length of the shape [m].
Source code in src/echosms/krmdata.py 
58
59
60
61
62
63
64
65
66
def length(self) -> float:
    """Length of the shape.

    Returns
    -------
    :
        The length of the shape [m].
    """
    return self.x[-1] - self.x[0]

volume()

Volume of the shape.

RETURNS DESCRIPTION
float

The volume of the shape [m³].
Source code in src/echosms/krmdata.py 
46
47
48
49
50
51
52
53
54
55
56
def volume(self) -> float:
    """Volume of the shape.

    Returns
    -------
    :
        The volume of the shape [m³].
    """
    thickness = np.diff(self.x)
    thickness = np.append(thickness, thickness[1])
    return np.sum(np.pi * (self.z_U - self.z_L) * self.w * thickness)

echosms.MSSModel()

Bases: ScatterModelBase

Modal series solution (MSS) scattering model.

This class calculates acoustic scatter from spheres and shells with various boundary conditions, as listed in the boundary_types class attribute.

Source code in src/echosms/mssmodel.py 
17
18
19
20
21
22
23
24
25
26
def __init__(self):
    super().__init__()
    self.long_name = 'modal series solution'
    self.short_name = 'mss'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled,
                           bt.fluid_shell_fluid_interior,
                           bt.fluid_shell_pressure_release_interior]
    self.shapes = ['sphere']
    self.max_ka = 20  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, f, boundary_type, target_c=None, target_rho=None, shell_c=None, shell_rho=None, shell_thickness=None, validate_parameters=True, **kwargs)

Calculate the scatter using the mss model for one set of parameters.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Radius of the spherical target [m].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The boundary type. Supported types are given in the boundary_types class variable.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the sphere [m/s]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

target_rho

Density of the fluid inside the sphere [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

shell_c

Sound speed in the spherical shell [m/s]. Only required for boundary_types that include a fluid shell.

TYPE: float DEFAULT: None

shell_rho

Density in the spherical shell [kg/m³]. Only required for boundary_types that include a fluid shell.

TYPE: float DEFAULT: None

shell_thickness

Thickness of the spherical shell [m]. This value is subtracted from a to give the radius of the interior sphere. Only required for boundary_types that include a fluid shell.

TYPE: float DEFAULT: None

validate_parameters

Whether to validate the model parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The class implements the code in Section A.1 of Jech et al. (2015).

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742–3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/mssmodel.py 
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
def calculate_ts_single(self, medium_c, medium_rho, a, f, boundary_type: bt,
                        target_c=None, target_rho=None,
                        shell_c=None, shell_rho=None, shell_thickness=None,
                        validate_parameters=True,
                        **kwargs) -> float:
    """
    Calculate the scatter using the mss model for one set of parameters.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    a : float
        Radius of the spherical target [m].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The boundary type. Supported types are given in the `boundary_types` class variable.
    target_c : float, optional
        Sound speed in the fluid inside the sphere [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho : float, optional
        Density of the fluid inside the sphere [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    shell_c : float, optional
        Sound speed in the spherical shell [m/s].
        Only required for `boundary_type`s that include a fluid shell.
    shell_rho : float, optional
        Density in the spherical shell [kg/m³].
        Only required for `boundary_type`s that include a fluid shell.
    shell_thickness : float, optional
        Thickness of the spherical shell [m]. This value is subtracted from ``a`` to give
        the radius of the interior sphere.
        Only required for `boundary_type`s that include a fluid shell.
    validate_parameters : bool
        Whether to validate the model parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The class implements the code in Section A.1 of Jech et al. (2015).

    References
    ----------
    Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N.,
    Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015.
    Comparisons among ten models of acoustic backscattering used in aquatic ecosystem
    research. Journal of the Acoustical Society of America 138, 3742–3764.
    <https://doi.org/10.1121/1.4937607>
    """
    if validate_parameters:
        self.validate_parameters(locals())

    k0 = wavenumber(medium_c, f)
    ka = k0*a
    n = np.arange(0, round(ka+20))

    match boundary_type:
        case bt.fixed_rigid:
            A = list(map(lambda x: -spherical_jn(x, ka, True) / h1(x, ka, True), n))
        case bt.pressure_release:
            A = list(map(lambda x: -spherical_jn(x, ka) / h1(x, ka), n))
        case bt.fluid_filled:
            k1a = wavenumber(target_c, f)*a
            gh = target_rho/medium_rho * target_c/medium_c

            def Cn_fr(n):
                return\
                    ((spherical_jn(n, k1a, True)*spherical_yn(n, ka))
                        / (spherical_jn(n, k1a)*spherical_jn(n, ka, True))
                        - gh*(spherical_yn(n, ka, True)/spherical_jn(n, ka, True)))\
                    / ((spherical_jn(n, k1a, True)*spherical_jn(n, ka))
                       / (spherical_jn(n, k1a)*spherical_jn(n, ka, True))-gh)

            A = -1/(1 + 1j*np.asarray(list(map(Cn_fr, n)), dtype=complex))
        case bt.fluid_shell_fluid_interior:
            b = a - shell_thickness

            g21 = shell_rho / medium_rho
            h21 = shell_c / medium_c
            g32 = target_rho / shell_rho
            h32 = target_c / shell_c

            k1a = wavenumber(medium_c, f) * a
            k2 = wavenumber(shell_c, f)
            k3b = wavenumber(target_c, f) * b

            def Cn_fsfi(n):
                (b1, b2, a11, a21, a12, a22, a32, a13, a23, a33) =\
                    MSSModel.__eqn9(n, k1a, g21, h21, k2*a, k2*b, k3b, h32, g32)
                return (b1*a22*a33 + a13*b2*a32 - a12*b2*a33 - b1*a23*a32)\
                    / (a11*a22*a33 + a13*a21*a32 - a12*a21*a33 - a11*a23*a32)

            A = list(map(Cn_fsfi, n))
        case bt.fluid_shell_pressure_release_interior:
            b = a - shell_thickness

            g21 = shell_rho / medium_rho
            h21 = shell_c / medium_c

            k1a = wavenumber(medium_c, f) * a
            k2 = wavenumber(shell_c, f)
            ksa = k2 * a  # ksa is used in the paper, but isn't that the same as k2a?

            def Cn_fspri(n):
                (b1, b2, d1, d2, a11, a21) = MSSModel.__eqn10(n, k1a, g21, h21, ksa, k2*a, k2*b)
                return (b1*d2-d1*b2) / (a11*d2-d1*a21)

            A = list(map(Cn_fspri, n))
        case _:
            raise ValueError(f'The {self.long_name} model does not support '
                             f'a model type of "{boundary_type}".')

    fbs = -1j/k0 * np.sum((-1)**n * (2*n+1) * A)
    return 20*log10(abs(fbs))  # ts

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/mssmodel.py 
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.

    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_rho', 'a', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        mask = p['boundary_type'] == t
        match t:
            case bt.fluid_filled:
                super()._present_and_positive(p, ['target_c', 'target_rho'], mask=mask)
            case bt.fluid_shell_fluid_interior:
                super()._present_and_positive(p, ['target_c', 'target_rho', 'shell_c',
                                                  'shell_rho', 'shell_thickness'], mask=mask)
            case bt.fluid_shell_pressure_release_interior:
                super()._present_and_positive(p, ['shell_c', 'shell_rho', 'shell_thickness'],
                                              mask=mask)

echosms.PSMSModel()

Bases: ScatterModelBase

Prolate spheroidal modal series (PSMS) scattering model.

Note

The fluid filled boundary type implementation is under development and is of limited use at the moment.

Source code in src/echosms/psmsmodel.py 
18
19
20
21
22
23
24
25
def __init__(self):
    super().__init__()
    self.long_name = 'prolate spheroidal modal series'
    self.short_name = 'psms'
    self.analytical_type = 'exact'
    self.boundary_types = [bt.fixed_rigid, bt.pressure_release, bt.fluid_filled]
    self.shapes = ['prolate spheroid']
    self.max_ka = 10  # [1]

calculate_ts(data, expand=False, inplace=False, multiprocess=False, progress=False)

Calculate the target strength (TS) for many parameters.

PARAMETER DESCRIPTION
data

Requirements for the different input data types are:

  • DataFrame: column names must match the function parameter names in calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
  • DataArray: dimension names must match the function parameter names in calculate_ts_single(). TS values will be calculated for all combinations of the coordinate variables.
  • dict: keys must match the function parameters in calculate_ts_single(). TS values will be calculated for all combinations of the dict values.

TYPE: Pandas DataFrame, Xarray DataArray or dict

multiprocess

Split the ts calculation across CPU cores. Multiprocessing is currently provided by mapply. For more sophisticated uses it may be preferred to use a multiprocessing package of your choice directly on the calculate_ts_single() method. See the code in this method (calculate_ts()) for an example.

TYPE: bool DEFAULT: False

expand

Only applicable if data is a dict. If True, will use as_dataframe() to expand the dict into a DataFrame with one column per dict key and return that, adding a column named ts for the results.

TYPE: bool DEFAULT: False

inplace

Only applicable if data is a DataFrame. If True, the results will be added to the input DataFrame in a column named ts. If a ts column already exists, it is overwritten.

TYPE: bool DEFAULT: False

progress

If True, will produce a progress bar while running models

TYPE: bool DEFAULT: False

RETURNS DESCRIPTION
None, list[float], Series, or DataFrame

The return type and value are determined by the type of the input variable (data) and the expand and inplace parameters:

  • dict input and expand=False returns a list of floats.
  • dict input and expand=True returns a DataFrame.
  • DataFrame input and inplace=False returns a Series.
  • DataFrame input and inplace=True modifies data and returns None.
  • DataArray input always modifies data and returns None.
Source code in src/echosms/scattermodelbase.py 
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
def calculate_ts(self, data, expand=False, inplace=False, multiprocess=False, progress=False):
    """Calculate the target strength (TS) for many parameters.

    Parameters
    ----------
    data : Pandas DataFrame, Xarray DataArray or dict
        Requirements for the different input data types are:

        - **DataFrame**: column names must match the function parameter names in
          calculate_ts_single(). One TS value will be calculated for each row in the DataFrame.
        - **DataArray**: dimension names must match the function parameter names in
          calculate_ts_single(). TS values will be calculated for all combinations of the
          coordinate variables.
        - **dict**: keys must match the function parameters in calculate_ts_single().
          TS values will be calculated for all combinations of the dict values.

    multiprocess : bool
        Split the ts calculation across CPU cores. Multiprocessing is currently provided by
        [mapply](https://github.com/ddelange/mapply). For more
        sophisticated uses it may be preferred to use a multiprocessing package of your choice
        directly on the `calculate_ts_single()` method. See the code in this method
        (`calculate_ts()`) for an example.

    expand : bool
        Only applicable if `data` is a dict. If `True`, will use
        [`as_dataframe()`][echosms.utils.as_dataframe]
        to expand the dict into a DataFrame with one column per dict key
        and return that, adding a column named `ts` for the results.

    inplace : bool
        Only applicable if `data` is a DataFrame. If `True`, the results
        will be added to the input DataFrame in a column named `ts`. If a `ts` column
        already exists, it is overwritten.

    progress : bool
        If `True`, will produce a progress bar while running models

    Returns
    -------
    : None, list[float], Series, or DataFrame
        The return type and value are determined by the type of the input variable (`data`) and
        the `expand` and `inplace` parameters:

        - dict input and `expand=False` returns a list of floats.
        - dict input and `expand=True` returns a DataFrame.
        - DataFrame input and `inplace=False` returns a Series.
        - DataFrame input and `inplace=True` modifies `data` and returns `None`.
        - DataArray input always modifies `data` and returns `None`.

    """
    match data:
        case dict():
            data_df = as_dataframe(data, self.no_expand_parameters)
        case pd.DataFrame():
            data_df = data
        case xr.DataArray():
            data_df = data.to_dataframe().reset_index()
            data_df.attrs = data.attrs
        case _:
            raise ValueError(f'Data type of {type(data)} is not supported'
                             ' (only dictionaries, Pandas DataFrames and'
                             ' Xarray DataArrays are).')

    self.validate_parameters(data_df)

    # Get the non-expandable model parameters
    p = data_df.attrs['parameters'] if 'parameters' in data_df.attrs else {}

    # Note: the args argument in the apply call below requires a tuple. data_df.attrs is a
    # dict and the default behaviour is to make a tuple using the dict keys. The trailing comma
    # and parenthesis instead causes the tuple to have one entry of the dict.

    if multiprocess:
        from mapply.mapply import mapply
        ts = mapply(data_df, self.__ts_helper, args=(p,), axis=1, progressbar=progress)
    else:  # this uses just one CPU
        if progress:
            tqdm.pandas(desc=self.short_name, unit=' models',
                        bar_format='{l_bar}{bar} [{n_fmt}/{total_fmt}; {rate_noinv_fmt}]')
            ts = data_df.progress_apply(self.__ts_helper, args=(p,), axis=1)
        else:
            ts = data_df.apply(self.__ts_helper, args=(p,), axis=1)

    match data:
        case dict() if expand:
            data_df['ts'] = ts
            return data_df
        case dict():
            return ts.to_list()
        case pd.DataFrame() if inplace:
            data_df['ts'] = ts
            return None
        case pd.DataFrame():
            return ts.rename('ts', inplace=True)
        case xr.DataArray():
            data.values = ts.to_numpy().reshape(data.shape)
            return None
        case _:
            raise AssertionError('This code should never be reached - unsupported input data '
                                 f'type of {type(data)}.')

calculate_ts_single(medium_c, medium_rho, a, b, theta, f, boundary_type, target_c=None, target_rho=None, validate_parameters=True)

Prolate spheroid modal series (PSMS) solution model.

PARAMETER DESCRIPTION
medium_c

Sound speed in the fluid medium surrounding the target [m/s].

TYPE: float

medium_rho

Density of the fluid medium surrounding the target [kg/m³].

TYPE: float

a

Prolate spheroid major axis radius [m].

TYPE: float

b

Prolate spheroid minor axis radius [m].

TYPE: float

theta

Pitch angle to calculate the scattering as per the echoSMs coordinate system [°].

TYPE: float

f

Frequency to calculate the scattering at [Hz].

TYPE: float

boundary_type

The model type. Supported model types are given in the boundary_types class variable.

TYPE: boundary_type

target_c

Sound speed in the fluid inside the target [m/s]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

target_rho

Density of the fluid inside the target [kg/m³]. Only required for boundary_type of fluid_filled.

TYPE: float DEFAULT: None

validate_parameters

Whether to validate the input parameters.

TYPE: bool DEFAULT: True

RETURNS DESCRIPTION
float

The target strength (re 1 m²) of the target [dB].
Notes

The backscattered target strength of a pressure release or fluid-filled prolate spheroid is calculated using the PSMS method of Furusawa (1988) and corrections in Furusawa et al. (1994).

References

Furusawa, M. (1988). "Prolate spheroidal models for predicting general trends of fish target strength," J. Acoust. Soc. Jpn. 9, 13-24. Furusawa, M., Miyanohana, Y., Ariji, M., and Sawada, Y. (1994). “Prediction of krill target strength by liquid prolate spheroid model,” Fish. Sci., 60, 261-265.

Source code in src/echosms/psmsmodel.py 
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def calculate_ts_single(self, medium_c, medium_rho, a, b, theta, f, boundary_type: bt,
                        target_c=None, target_rho=None, validate_parameters=True):
    """Prolate spheroid modal series (PSMS) solution model.

    Parameters
    ----------
    medium_c : float
        Sound speed in the fluid medium surrounding the target [m/s].
    medium_rho : float
        Density of the fluid medium surrounding the target [kg/m³].
    a : float
        Prolate spheroid major axis radius [m].
    b : float
        Prolate spheroid minor axis radius [m].
    theta : float
        Pitch angle to calculate the scattering as per the echoSMs
        [coordinate system](https://ices-tools-dev.github.io/echoSMs/
        conventions/#coordinate-systems) [°].
    f : float
        Frequency to calculate the scattering at [Hz].
    boundary_type :
        The model type. Supported model types are given in the `boundary_types` class variable.
    target_c : float
        Sound speed in the fluid inside the target [m/s].
        Only required for `boundary_type` of ``fluid_filled``.
    target_rho : float
        Density of the fluid inside the target [kg/m³].
        Only required for `boundary_type` of ``fluid_filled``.
    validate_parameters : bool
        Whether to validate the input parameters.

    Returns
    -------
    : float
        The target strength (re 1 m²) of the target [dB].

    Notes
    -----
    The backscattered target strength of a pressure release or fluid-filled prolate spheroid
    is calculated using the PSMS method of Furusawa (1988) and corrections in
    Furusawa et al. (1994).

    References
    ----------
    Furusawa, M. (1988). "Prolate spheroidal models for predicting general
        trends of fish target strength," J. Acoust. Soc. Jpn. 9, 13-24.
    Furusawa, M., Miyanohana, Y., Ariji, M., and Sawada, Y. (1994).
        “Prediction of krill target strength by liquid prolate spheroid
        model,” Fish. Sci., 60, 261-265.
    """
    if validate_parameters:
        self.validate_parameters(locals())

    if boundary_type not in self.boundary_types:
        raise ValueError(f'The {self.long_name} model does not support '
                         f'a model type of "{boundary_type}".')

    xim = (1.0 - (b/a)**2)**(-0.5)
    q = a/xim  # semi-focal length

    km = wavenumber(medium_c, f)
    hm = km*q

    if boundary_type == bt.fluid_filled:
        g = target_rho / medium_rho
        ht = wavenumber(target_c, f)*q
        simplified = False
        # Note: we can implement simpler equations if sound speeds are
        # similar between the medium and the target. The simplified
        # equations are quicker, so it is worth to do.
        # But, it appears that target_c ≈ medium_c is not the only requirement for
        # the simplification to work well - see section 4.1.1 in:
        #
        # Gonzalez, J. D., Lavia, E. F., Blanc, S., & Prario, I. (2016).
        # Acoustic scattering by prolate and oblate liquid spheroids.
        # Proceedings of the 22nd International Congress on Acoustics.
        # Acoustics for the 21st Century, Buenos Aires.
        #
        # So, for the moment, the simplification is turned off.

        # if abs(1.0-target_c/medium_c) <= 0.01:
        #    simplified = True

    # Phi, the roll angle is fixed for this implementation
    phi_inc = np.pi  # incident direction
    phi_sca = np.pi + phi_inc  # scattered direction

    theta_inc = np.deg2rad(theta)  # incident direction
    theta_sca = np.pi - theta_inc  # scattered direction

    # Approximate limits on the summations. These come from Furusawa (1988), but other
    # implementations of this code use more complex functions to calculate the maximum orders
    # of spheroidal wave functions to calculate. It is also feasible to work to a defined
    # convergence level. This is a potenital future improvement.
    m_max = int(np.ceil(2*km*b))  # +1
    n_max = int(m_max + np.ceil(hm/2))  # +3

    f_sca = 0.0
    for m in range(m_max+1):
        epsilon_m = Neumann(m)
        cos_term = np.cos(m*(phi_sca - phi_inc))

        if boundary_type == bt.fluid_filled and not simplified:
            Am = PSMSModel._fluidfilled(m, n_max, hm, ht, xim, g, theta_inc)

        for n_i, n in enumerate(range(m, n_max+1)):
            # Use the normalisation offered by spheroidalwavefunctions.pro_ang1() because
            # it removes the need to calculate a normalisation factor (N_mn in Furusawa, 1998).
            # This is because the pro_ang1(norm=True) norm is unity.
            Smn_inc, _ = pro_ang1(m, n, hm, np.cos(theta_inc), norm=True)
            Smn_sca, _ = pro_ang1(m, n, hm, np.cos(theta_sca), norm=True)

            # FYI, the code used to use the Meixner-Schäfke normalisation scheme for the
            # angular function of the first kind (eqn 21.7.11 in Abramowitz, M., and Stegun,
            # I. A. (1964). Handbook of Mathematical Functions with Formulas, Graphs, and
            # Mathematical Tables # (Dover, New York), 10th ed.
            # N_mn = 2/(2*n+1) * factorial(n+m) / factorial(n-m)

            R1m, dR1m = pro_rad1(m, n, hm, xim)
            R2m, dR2m = pro_rad2(m, n, hm, xim)

            match boundary_type:
                case bt.fluid_filled:
                    if simplified:
                        E1, E3 = PSMSModel._fluidfilled_Emn(m, n, n, hm, ht, xim, g)
                        Amn = -E1/E3
                    else:
                        Amn = Am[n_i][0]
                case bt.pressure_release:
                    Amn = -R1m/(R1m + 1j*R2m)
                case bt.fixed_rigid:
                    Amn = -dR1m/(dR1m + 1j*dR2m)

            f_sca += epsilon_m * Smn_inc * Amn * Smn_sca * cos_term

    return 20*np.log10(np.abs(-2j / km * f_sca))

validate_parameters(params)

Validate the model parameters.

See here for calling details.

Source code in src/echosms/psmsmodel.py 
27
28
29
30
31
32
33
34
35
36
37
38
39
40
def validate_parameters(self, params):
    """Validate the model parameters.

    See [here][echosms.ScatterModelBase.validate_parameters] for calling details.
    """
    p = as_dict(params)
    super()._present_and_in(p, ['boundary_type'], self.boundary_types)
    super()._present_and_positive(p, ['medium_c', 'medium_rho', 'a', 'b', 'f'])

    types = np.unique(np.atleast_1d(p['boundary_type']))
    for t in types:
        if t == bt.fluid_filled:
            super()._present_and_positive(p, ['target_c', 'target_rho'],
                                          mask=p['boundary_type'] == t)

Utilities

echosms.BenchmarkData()

Convenient interface to the benchmark dataset.

This dataset contains the benchmark TS results from Jech et al. (2015).

Notes

The column names in the source benchmark files have been changed to be the same as those used in the ReferenceModels model definitions.

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742-3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/benchmarkdata.py 
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
def __init__(self):

    data_directory = Path(__file__).parent/Path('resources')/Path('BenchMark_Data')

    angle_data_file = data_directory/'Benchmark_Angle_TS.csv'
    freq_data_file = data_directory/'Benchmark_Frequency_TS.csv'

    self.angle_dataset = pd.read_csv(angle_data_file)
    self.freq_dataset = pd.read_csv(freq_data_file)

    # Change the column names to match the reference model names used in ReferenceModels
    self.angle_dataset.rename(columns=BenchmarkData.a_rename, inplace=True)
    self.freq_dataset.rename(columns=BenchmarkData.f_rename, inplace=True)

    self.freq_dataset['frequency (kHz)'] *= 1e3  # want Hz not kHz

    # Remove units from the column names (we have the echoSMs units convention instead)
    self.freq_dataset.rename(columns={'frequency (kHz)': 'frequency'}, inplace=True)
    self.angle_dataset.rename(columns={'angle (deg)': 'angle'}, inplace=True)

    self.angle_dataset.set_index('angle', inplace=True)
    self.freq_dataset.set_index('frequency', inplace=True)

angle_as_dataframe()

Provide the angle benchmark dataset as a Pandas DataFrame.

RETURNS DESCRIPTION
DataFrame

Dataframe containing the benchmark data.
Source code in src/echosms/benchmarkdata.py 
133
134
135
136
137
138
139
140
141
def angle_as_dataframe(self) -> pd.DataFrame:
    """Provide the angle benchmark dataset as a Pandas DataFrame.

    Returns
    -------
    :
        Dataframe containing the benchmark data.
    """
    return self.angle_dataset

angle_data(name)

Provide the benchmark TS values verses angle for the name model.

PARAMETER DESCRIPTION
name

The name of the benchmark model (available from angle_names()).

TYPE: str

RETURNS DESCRIPTION
tuple

Tuple containing the angles (°) and TS (dB) for the requested benchmark model.
Source code in src/echosms/benchmarkdata.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
def angle_data(self, name: str) -> tuple:
    """Provide the benchmark TS values verses angle for the `name` model.

    Parameters
    ----------
    name :
        The name of the benchmark model (available from `angle_names()`).

    Returns
    -------
    :
        Tuple containing the angles (°) and TS (dB) for the requested benchmark model.
    """
    if name not in self.angle_names():
        raise ValueError(f'The requested model ({name}) is not in the angle benchmark dataset.')
    return (self.angle_dataset.index.values, self.angle_dataset[name].values)

angle_names()

Provide the model names for the angle benchmark data.

RETURNS DESCRIPTION
list

List of model names.
Source code in src/echosms/benchmarkdata.py 
77
78
79
80
81
82
83
84
85
86
def angle_names(self) -> list:
    """Provide the model names for the angle benchmark data.

    Returns
    -------
    :
        List of model names.

    """
    return self.angle_dataset.columns.values.tolist()

freq_as_dataframe()

Provide the frequency benchmark dataset as a Pandas DataFrame.

RETURNS DESCRIPTION
DataFrame

Dataframe containing the benchmark data.
Source code in src/echosms/benchmarkdata.py 
143
144
145
146
147
148
149
150
151
def freq_as_dataframe(self) -> pd.DataFrame:
    """Provide the frequency benchmark dataset as a Pandas DataFrame.

    Returns
    -------
    :
        Dataframe containing the benchmark data.
    """
    return self.freq_dataset

freq_data(name)

Provide the benchmark TS values verses frequency for the name model.

PARAMETER DESCRIPTION
name

The name of the benchmark model (available from freq_names()).

TYPE: str

RETURNS DESCRIPTION
tuple

Tuple containing the frequencies (Hz) and TS (dB) for the requested benchmark model.
Source code in src/echosms/benchmarkdata.py 
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def freq_data(self, name: str) -> tuple:
    """Provide the benchmark TS values verses frequency for the `name` model.

    Parameters
    ----------
    name :
        The name of the benchmark model (available from `freq_names()`).

    Returns
    -------
    :
        Tuple containing the frequencies (Hz) and TS (dB) for the requested benchmark model.
    """
    if name not in self.freq_names():
        raise ValueError(f'The requested model ({name}) '
                         'is not in the frequency benchmark dataset.')
    return (self.freq_dataset.index.values, self.freq_dataset[name].values)

freq_names()

Provide the model names for the frequency benchmark data.

RETURNS DESCRIPTION
list

List of model names.
Source code in src/echosms/benchmarkdata.py 
88
89
90
91
92
93
94
95
96
def freq_names(self) -> list:
    """Provide the model names for the frequency benchmark data.

    Returns
    -------
    :
        List of model names.
    """
    return self.freq_dataset.columns.values.tolist()

echosms.ReferenceModels()

Provide access to reference scattering model parameters.

Reference models are the models and parameters defined in Jech et al. (2015). The parameters are stored in a TOML-formatted file in the echoSMs repository and this class provides easy access to the data in that file. Additional reference models may be defined in the future and added to the TOML file (for example, entries have been added for all known calibration sphere sizes).

ATTRIBUTE DESCRIPTION
definitions

A dict representation of the target definitions.toml file.

TYPE: dict

RAISES DESCRIPTION
TOMLDecodeError

If the target definitions.toml file is not valid TOML.
KeyError

If the target definitions.toml file has multiple target entries with the same name.
References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742–3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/referencemodels.py 
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def __init__(self):
    self.defs_filename = Path(__file__).parent/Path('resources')/Path('target definitions.toml')

    self.definitions = []

    with open(self.defs_filename, 'rb') as f:
        try:
            self.definitions = tomllib.load(f)
        except tomllib.TOMLDecodeError as e:
            raise SyntaxError(f'Error while parsing file "{self.defs_filename.name}"') from e

    # Flag duplicate target names
    pda = pd.Series(self.names())
    duplicates = list(pda[pda.duplicated()])
    if duplicates:
        raise KeyError(f'The "{self.defs_filename.name}" file has multiple targets '
                       f'with the same name: '+', '.join(duplicates))

    # Substitute parameters names in the target section by the values of those
    # parameters.
    for t in self.definitions['target']:
        for k, v in t.items():
            try:
                t[k] = self.definitions['parameters'][v]
            except (KeyError, TypeError):
                pass

names()

Names of all model definitions.

RETURNS DESCRIPTION
iterable of str

All model names in the target definitions.toml file.
Source code in src/echosms/referencemodels.py 
69
70
71
72
73
74
75
76
77
def names(self):
    """Names of all model definitions.

    Returns
    -------
    : iterable of str
        All model names in the ``target definitions.toml`` file.
    """
    return [n['name'] for n in self.definitions['target']]

parameters(name)

Model parameters for a particular model.

Model parameters are a subset of the model specification where the metadata items have been removed.

PARAMETER DESCRIPTION
name

The name of a model in the target definitions.toml file.

TYPE: str

RETURNS DESCRIPTION
dict

The model parameters for the requested model or an empty dict if no model with that name.
Source code in src/echosms/referencemodels.py 
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
def parameters(self, name: str) -> dict:
    """Model parameters for a particular model.

    Model parameters are a subset of the model specification where the metadata items have
    been removed.

    Parameters
    ----------
    name :
        The name of a model in the ``target definitions.toml`` file.

    Returns
    -------
    :
        The model parameters for the requested model or an empty dict if no model with
        that name.

    """
    s = self.specification(name)

    if not s:
        return {}

    # Remove the entries that are not parameters
    p = s.copy()
    for k in ['name', 'shape', 'description', 'source', 'benchmark_model']:
        p.pop(k, None)
    return p

specification(name)

Model definitions for a particular model.

PARAMETER DESCRIPTION
name

The name of a model in the target definitions.toml file.

TYPE: str

RETURNS DESCRIPTION
dict

The model definitions for the requested model or an empty dict if no model with that name.
Source code in src/echosms/referencemodels.py 
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
def specification(self, name: str) -> dict:
    """Model definitions for a particular model.

    Parameters
    ----------
    name :
        The name of a model in the ``target definitions.toml`` file.

    Returns
    -------
    :
        The model definitions for the requested model or an empty dict if no model
        with that name.
    """
    s = [t for t in self.definitions['target'] if t['name'] == name]
    if not s:
        return {}

    return s[0]

echosms.JechEtAlData()

Simple access to all model results in Jech et al (2015).

ATTRIBUTE DESCRIPTION
data

One entry in the dict for each model results file.

TYPE: dict[DataFrame]

data_directory

The directory containing the model results files.

TYPE: Path

References

Jech, J.M., Horne, J.K., Chu, D., Demer, D.A., Francis, D.T.I., Gorska, N., Jones, B., Lavery, A.C., Stanton, T.K., Macaulay, G.J., Reeder, D.B., Sawada, K., 2015. Comparisons among ten models of acoustic backscattering used in aquatic ecosystem research. Journal of the Acoustical Society of America 138, 3742-3764. https://doi.org/10.1121/1.4937607

Source code in src/echosms/jechetaldata.py 
25
26
27
def __init__(self):
    self.data_directory = Path(__file__).parent/Path('resources')/Path('Jechetal_allmodels')
    self.data = {f.stem: pd.read_csv(f) for f in self.data_directory.glob('*.csv')}

echosms.utils

Miscellaneous utility functions.

boundary_type

Bases: StrEnum

Scattering model boundary types.

fixed_rigid = 'fixed-rigid' class-attribute instance-attribute

Fixed-rigid boundary condition, where the normal velocity is zero at the object's boundary.

pressure_release = 'pressure-release' class-attribute instance-attribute

Pressure-release boundary condition, where the acoustic pressure is zero at the object's boundary.

fluid_filled = 'fluid-filled' class-attribute instance-attribute

Fluid-filled boundary condition, where the acoustic pressure and normal velocity are both non-zero at the object's boundary.

elastic = 'elastic' class-attribute instance-attribute

The scattering object supports compressional and shear waves.

fluid_shell_fluid_interior = 'fluid shell fluid interior' class-attribute instance-attribute

The object has a fluid interior surrounded by a fluid shell.

fluid_shell_pressure_release_interior = 'fluid shell pressure release interior' class-attribute instance-attribute

The object has a pressure release interior surrounded by a fluid shell.

hard = 'fixed-rigid' class-attribute instance-attribute

A synonym for fixed_rigid.

soft = 'pressure-release' class-attribute instance-attribute

A synonym for pressure_release.

fluid = 'fluid-filled' class-attribute instance-attribute

A synonym for fluid_filled.

theoretical_Sa(ts, eba, r, nautical=False)

Theoretical area backscattering strength (Sₐ) of a target.

PARAMETER DESCRIPTION
ts

The target strength of the object [dB re 1 m²].

TYPE: float | ndarray

eba

Ten times the logarithm to the base 10 of the transducer's equivalent two-way beam angle (ψ, sr). In formula form this is: eba = 10 log10(ψ) dB (re 1 sr).

TYPE: float

r

The range from the transducer to the target [m]. Used for acoustic beam spreading.

TYPE: float

nautical

If True, the nautical area scattering strength (SA) is returned instead of the area backscattering strength (Sₐ).

DEFAULT: False

RETURNS DESCRIPTION
float | ndarray

The theoretical Sₐ [dB re 1 m² m⁻²] or SA [dB re 1 m² nmi⁻²] of the input ts.
RAISES DESCRIPTION
ValueError

If an input value is out of bounds.
Notes

While the calculation is valid for any target, the theoretical area strengths are most relevant when calibrating an echosounder using a sphere. The difference between the theoretical and measured can be used to calculate the calibration gain for an echosounder (when the sphere is on-axis).

Source code in src/echosms/utils.py 
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
def theoretical_Sa(ts: float | np.ndarray, eba: float, r: float, nautical=False)\
                   -> float | np.ndarray:
    """Theoretical area backscattering strength (Sₐ) of a target.

    Parameters
    ----------
    ts :
        The target strength of the object [dB re 1 m²].
    eba :
        Ten times the logarithm to the base 10 of the transducer's equivalent two-way beam angle (ψ, sr).
        In formula form this is: eba = 10 log<sub>10</sub>(ψ) dB (re 1 sr).
    r :
        The range from the transducer to the target [m]. Used for acoustic beam spreading.
    nautical :
        If `True`, the nautical area scattering strength (S<sub>A</sub>) is returned instead of the
        area backscattering strength (Sₐ).

    Returns
    -------
    :
        The theoretical Sₐ [dB re 1 m² m⁻²] or S<sub>A</sub> [dB re 1 m² nmi⁻²] of the input `ts`.

    Raises
    ------
    ValueError
        If an input value is out of bounds.

    Notes
    -----
    While the calculation is valid for any target, the theoretical area strengths are most relevant
    when calibrating an echosounder using a sphere. The difference between
    the theoretical and measured can be used to calculate the calibration gain for an
    echosounder (when the sphere is on-axis).
    """
    if eba > 0.0:
        raise ValueError('A positive eba value is not supported.')
    if r <= 0.0:
        raise ValueError('An r value less than or equal to 0 is not supported.')

    factor = 10.0*np.log10(4.0*π*1852.0**2) if nautical else 0.0
    return ts - eba - 20.0*np.log10(r) + factor

Neumann(m)

Neumann number.

PARAMETER DESCRIPTION
m

The input integer.

TYPE: int

RETURNS DESCRIPTION
int

The Neumann number.
Source code in src/echosms/utils.py 
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
def Neumann(m: int) -> int:
    """Neumann number.

    Parameters
    ----------
    m :
        The input integer.

    Returns
    -------
    :
        The Neumann number.
    """
    if m == 0:
        return 1
    return 2

wavenumber(c, f)

Calculate the acoustic wavenumber.

PARAMETER DESCRIPTION
c

Sound speed [m/s]

TYPE: float

f

Frequency [Hz]

TYPE: float

RETURNS DESCRIPTION
float

The acoustic wavenumber [m⁻¹].
Source code in src/echosms/utils.py 
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
def wavenumber(c: float, f: float) -> float:
    """Calculate the acoustic wavenumber.

    Parameters
    ----------
    c :
        Sound speed [m/s]

    f :
        Frequency [Hz]

    Returns
    -------
    :
        The acoustic wavenumber [m⁻¹].
    """
    return 2*π*f/c

wavelength(c, f)

Calculate the acoustic wavelength.

PARAMETER DESCRIPTION
c

Sound speed [m/s]

TYPE: float

f

Frequency [Hz]

TYPE: float

RETURNS DESCRIPTION
float

The acoustic wavelength [m].
Source code in src/echosms/utils.py 
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
def wavelength(c: float, f: float) -> float:
    """Calculate the acoustic wavelength.

    Parameters
    ----------
    c :
        Sound speed [m/s]

    f :
        Frequency [Hz]

    Returns
    -------
    :
        The acoustic wavelength [m].
    """
    return c/f

h1(n, z, derivative=False)

Spherical Hankel function of the first kind or its' derivative.

PARAMETER DESCRIPTION
n

Order (n ≥ 0).

TYPE: int

z

Argument of the Hankel function.

TYPE: float

derivative

if True, the value of the derivative (rather than the function itself) is returned.

DEFAULT: False

RETURNS DESCRIPTION
complex

Value of the spherical Hankel function
RAISES DESCRIPTION
ValueError

For negative n values.
Notes

The value of the Hankel function is calculated from spherical Bessel functions [1].

The derivative is computed from spherical Hankel functions [2].

References

[1] https://dlmf.nist.gov/10.47.E10

[2] https://dlmf.nist.gov/10.51.E2

Source code in src/echosms/utils.py 
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
def h1(n: int, z: float, derivative=False) -> complex:
    """Spherical Hankel function of the first kind or its' derivative.

    Parameters
    ----------
    n :
        Order (n ≥ 0).
    z :
        Argument of the Hankel function.
    derivative :
        if True, the value of the derivative (rather than the function itself) is returned.

    Returns
    -------
    :
        Value of the spherical Hankel function

    Raises
    ------
    ValueError
        For negative n values.

    Notes
    -----
    The value of the Hankel function is calculated from spherical Bessel functions [1].

    The derivative is computed from spherical Hankel functions [2].

    References
    ----------
    [1] <https://dlmf.nist.gov/10.47.E10>

    [2] <https://dlmf.nist.gov/10.51.E2>
    """
    if n < 0:
        raise ValueError('Negative n values are not supported for spherical Hankel functions.')

    if not derivative:
        return spherical_jn(n, z) + 1j*spherical_yn(n, z)
    return -h1(n+1, z) + (n/z) * h1(n, z)

spherical_jnpp(n, z)

Second derivative of the spherical Bessel function.

PARAMETER DESCRIPTION
n

Order (n ≥ 0)

TYPE: int

z

Argument of the Bessel function.

TYPE: float

RETURNS DESCRIPTION
float

The second derivative of the spherical Bessel function.
Source code in src/echosms/utils.py 
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
def spherical_jnpp(n: int, z: float) -> float:
    """Second derivative of the spherical Bessel function.

    Parameters
    ----------
    n :
        Order (n ≥ 0)
    z :
        Argument of the Bessel function.

    Returns
    -------
    :
        The second derivative of the spherical Bessel function.

    """
    return 1./z**2 * ((n**2-n-z**2)*spherical_jn(n, z) + 2.*z*spherical_jn(n+1, z))

split_dict(d, s)

Split a dict into two dicts based on a list of keys.

PARAMETER DESCRIPTION
d

Dict to be split.

TYPE: dict

s

List of dict keys to use for splitting d.

TYPE: list

RETURNS DESCRIPTION
tuple(dict, dict)

The input dict split into two dicts based on the keys in s. The first tuple item contains the items that do not have keys in s.
Source code in src/echosms/utils.py 
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
def split_dict(d: dict, s: list) -> tuple[dict, dict]:
    """Split a dict into two dicts based on a list of keys.

    Parameters
    ----------
    d : dict
        Dict to be split.

    s: list
        List of dict keys to use for splitting `d`.

    Returns
    -------
    : tuple(dict, dict)
        The `input` dict split into two dicts based on the keys in `s`. The first tuple item
        contains the items that do not have keys in `s`.
    """
    contains = {k: v for k, v in d.items() if k in s}
    ncontains = {k: v for k, v in d.items() if k not in s}
    return ncontains, contains

as_dataarray(params, no_expand=[])

Convert model parameters from dict form to a Xarray DataArray.

PARAMETER DESCRIPTION
params

The model parameters.

TYPE: dict

no_expand

Key values of the non-expandable model parameters in params.

TYPE: list DEFAULT: []

RETURNS DESCRIPTION
DataArray

Returns a multi-dimensional DataArray generated from the Cartesian product of all expandable items in the input dict. Non-expandable items are added to the DataArray attrs property. Expandable items are those that can be sensibly expanded into DataArray coordinates. Not all models have non-expandable items. The array is named ts, the values are initialised to nan, the dimension names are the dict keys, and the coordinate variables are the dict values.
Source code in src/echosms/utils.py 
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
def as_dataarray(params: dict, no_expand: list = []) -> xr.DataArray:
    """Convert model parameters from dict form to a Xarray DataArray.

    Parameters
    ----------
    params :
        The model parameters.

    no_expand :
        Key values of the non-expandable model parameters in `params`.

    Returns
    -------
    :
        Returns a multi-dimensional DataArray generated from the Cartesian product of all
        expandable items in the input dict. Non-expandable items are added to the DataArray
        attrs property. Expandable items are those that can be sensibly expanded into
        DataArray coordinates. Not all models have non-expandable items.
        The array is named `ts`, the values are initialised to `nan`, the
        dimension names are the dict keys, and the coordinate variables are the dict values.

    """
    expand, nexpand = split_dict(params, no_expand)

    # Convert scalars to iterables so xarray is happy
    for k, v in expand.items():
        if not isinstance(v, Iterable) or isinstance(v, str):
            expand[k] = [v]

    sz = [len(v) for k, v in expand.items()]
    return xr.DataArray(data=np.full(sz, np.nan), coords=expand, name='ts',
                        attrs={'units': 'dB', 'dB_reference': '1 m^2',
                               'parameters': nexpand})

as_dataframe(params, no_expand=[])

Convert model parameters from dict form to a Pandas DataFrame.

PARAMETER DESCRIPTION
params

The model parameters.

TYPE: dict

no_expand

Key values of the non-expandable model parameters in params.

TYPE: list DEFAULT: []

RETURNS DESCRIPTION
DataFrame

Returns a Pandas DataFrame generated from the Cartesian product of all expandable items in the input dict. DataFrame column names are obtained from the dict keys. Non-expandable items are added to the DataFrame attrs property. Expandable items are those that can be sensibly expanded into DataFrame columns. Not all models have non-expandable items.
Source code in src/echosms/utils.py 
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
def as_dataframe(params: dict, no_expand: list = []) -> pd.DataFrame:
    """Convert model parameters from dict form to a Pandas DataFrame.

    Parameters
    ----------
    params :
        The model parameters.

    no_expand :
        Key values of the non-expandable model parameters in `params`.

    Returns
    -------
    :
        Returns a Pandas DataFrame generated from the Cartesian product of all expandable
        items in the input dict. DataFrame column names are obtained from the dict keys.
        Non-expandable items are added to the DataFrame attrs property. Expandable items are
        those that can be sensibly expanded into DataFrame columns. Not all models have
        non-expandable items.

    """
    expand, nexpand = split_dict(params, no_expand)

    # Use meshgrid to do the Cartesian product then create a Pandas DataFrame from that, having
    # flattened the multidimensional arrays and using a dict to provide column names.
    # This preserves the differing dtypes in each column compared to other ways of
    # constructing the DataFrame).
    df = pd.DataFrame({k: t.flatten()
                       for k, t in zip(expand.keys(), np.meshgrid(*tuple(expand.values())))})
    df.attrs = {'parameters': nexpand}
    return df

as_dict(params)

Convert model parameters from DataFrame or DataArray to dict.

PARAMETER DESCRIPTION
params

The model parameters

TYPE: dict | DataFrame | DataArray

RAISES DESCRIPTION
TypeError:

If the input data type is not supported.
RETURNS DESCRIPTION
dict

A dict containing the model parameters.
Source code in src/echosms/utils.py 
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
def as_dict(params: dict | pd.DataFrame | xr.DataArray) -> dict:
    """Convert model parameters from DataFrame or DataArray to dict.

    Parameters
    ----------
    params:
        The model parameters

    Raises
    ------
    TypeError:
        If the input data type is not supported.

    Returns
    -------
    :
        A dict containing the model parameters.
    """
    if isinstance(params, dict):
        return params

    # Get the non-expandable model parameters
    p = params.attrs['parameters'] if 'parameters' in params.attrs else {}

    if isinstance(params, xr.DataArray):
        return dict(zip(params.coords, params.indexes.values())) | p
    elif isinstance(params, pd.DataFrame):
        # params.attrs = {}  # otherwise to_dict() exposes a bug in pandas?
        return params.to_dict(orient='series') | p

    raise TypeError('Only dict, DataFrame, or DataArray are accepted.')

pro_ang1(m, n, c, eta, norm=False)

Prolate spheroidal angular function of the first kind and derivative.

Calculates the prolate spheroidal angular function of the first kind and its' derivative with respect to eta.

PARAMETER DESCRIPTION
m

The order parameter (≥ 0)

TYPE: int

n

The degree parameter (≥ m).

TYPE: int

c

The size parameter.

TYPE: float

eta

The angular coordinate, η, where |η| ≤ 1.

TYPE: float

norm

If False, returned values are not normalised (i.e., the Meixner-Schäfke normlalisation scheme is used). For large m this norm becomes very large. If True the returned values are scaled by the square root of the normalisation of the corresponding Legendre function. This avoids the large values that occur when norm is False.

DEFAULT: False

RETURNS DESCRIPTION
tuple[float, float]

The value of the prolate spheroidal angular function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
def pro_ang1(m: int, n: int, c: float, eta: float, norm=False) -> tuple[float, float]:
    """Prolate spheroidal angular function of the first kind and derivative.

    Calculates the prolate spheroidal angular function of the first kind and its'
    derivative with respect to `eta`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0)
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    eta :
        The angular coordinate, η, where |η| ≤ 1.
    norm :
        If `False`, returned values are not normalised (i.e., the Meixner-Schäfke normlalisation
        scheme is used). For large `m` this norm becomes very large. If `True` the returned values
        are scaled by the square root of the normalisation of the corresponding Legendre function.
        This avoids the large values that occur when `norm` is `False`.

    Returns
    -------
    :
        The value of the prolate spheroidal angular function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if abs(eta) > 1.0:
        raise ValueError('The eta parameter must be less than or equal to 1')

    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=0.0, ioprad=0, iopang=2,
                           iopnorm=int(norm), arg=[eta])
    p = swf_t._make(a)
    if np.isnan(p.s1c[n-m]) or np.isnan(p.s1dc[n-m]):
        # print('Root - trying again.')
        a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=0.0, ioprad=0, iopang=2,
                               iopnorm=int(norm), arg=[eta+sys.float_info.epsilon])
        p = swf_t._make(a)

    s = p.s1c[n-m] * np.float_power(10.0, p.is1e[n-m])
    sp = p.s1dc[n-m] * np.float_power(10.0, p.is1de[n-m])

    return s[0], sp[0]

pro_rad1(m, n, c, xi)

Prolate spheroidal radial function of the first kind and derivative.

Calculates the prolate spheroidal radial function of the first kind and its' derivative with respect to xi.

PARAMETER DESCRIPTION
m

The order parameter (≥ 0).

TYPE: int

n

The degree parameter (≥ m).

TYPE: int

c

The size parameter.

TYPE: float

xi

The radial coordinate, ξ, where ξ ≥ 1.

TYPE: float

RETURNS DESCRIPTION
tuple[float, float]

The value of the prolate spheroidal radial function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
def pro_rad1(m: int, n: int, c: float, xi: float) -> tuple[float, float]:
    """Prolate spheroidal radial function of the first kind and derivative.

    Calculates the prolate spheroidal radial function of the first kind and its'
    derivative with respect to `xi`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0).
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    xi :
        The radial coordinate, ξ, where ξ ≥ 1.

    Returns
    -------
    :
        The value of the prolate spheroidal radial function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if xi < 1.0:
        raise ValueError('The xi parameter must be greater than or equal to 1')

    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=xi-1.0, ioprad=1, iopang=0, iopnorm=0, arg=[0])
    p = swf_t._make(a)
    s = p.r1c * np.float_power(10.0, p.ir1e)
    sp = p.r1dc * np.float_power(10.0, p.ir1de)

    return s[n-m], sp[n-m]

pro_rad2(m, n, c, xi)

Prolate spheroidal radial function of the second kind and derivative.

Calculates the prolate spheroidal radial function of the second kind and its' derivative with respect to xi.

PARAMETER DESCRIPTION
m

The order parameter (≥ 0).

TYPE: int

n

The degree parameter (≥ m).

TYPE: int

c

The size parameter.

TYPE: float

xi

The radial coordinate, ξ, where ξ ≥ 1.

TYPE: float

RETURNS DESCRIPTION
tuple[float, float]

The value of the prolate spheroidal radial function and its' derivative.
Notes

This method uses the prolate spheroidal wave function code for non complex arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on github. This code is in Fortran90 and was interfaced to Python using numpy.f2py then wrapped with the current method to provide a similar calling convention as the Scipy function of the same name.

References

Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal radial functions of the first kind and their first derivatives. Quarterly of Applied Mathematics, 60(3), 589-599. https://doi.org/10.1090/qam/1914443

Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied Mathematics, 62(3), 493-507. https://doi.org/10.1090/qam/2086042

Source code in src/echosms/utils.py 
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
def pro_rad2(m: int, n: int, c: float, xi: float) -> tuple[float, float]:
    """Prolate spheroidal radial function of the second kind and derivative.

    Calculates the prolate spheroidal radial function of the second kind and its'
    derivative with respect to `xi`.

    Parameters
    ----------
    m :
        The order parameter (≥ 0).
    n :
        The degree parameter (≥ `m`).
    c :
        The size parameter.
    xi :
        The radial coordinate, ξ, where ξ ≥ 1.

    Returns
    -------
    :
        The value of the prolate spheroidal radial function and its' derivative.

    Notes
    -----
    This method uses the prolate spheroidal wave function code for non complex
    arguments (van Buren & Boisvert, 2002, and van Buren & Boisvert, 2024), available on
    [github](https://github.com/MathieuandSpheroidalWaveFunctions). This code is in Fortran90
    and was interfaced to Python using `numpy.f2py` then wrapped with the current method to
    provide a similar calling convention as the Scipy function of the same name.

    References
    ----------
    Van Buren, A. L., & Boisvert, J. E. (2002). Accurate calculation of prolate spheroidal
    radial functions of the first kind and their first derivatives. Quarterly of Applied
    Mathematics, 60(3), 589-599. <https://doi.org/10.1090/qam/1914443>

    Van Buren, A. L., & Boisvert, J. E. (2004). Improved Calculation of Prolate Spheroidal
    Radial Functions of the Second Kind and Their First Derivatives. Quarterly of Applied
    Mathematics, 62(3), 493-507. <https://doi.org/10.1090/qam/2086042>
    """
    if m < 0:
        raise ValueError('The m parameter must be positive.')
    if n < m:
        raise ValueError('The n parameter must be greater than or equal to the m parameter.')
    if xi < 1.0:
        raise ValueError('The xi parameter must be greater than or equal to 1')

    ioprad = 1 if xi-1.0 < 1e-10 else 2

    # Add +2 to lnum instead of +1 as it exposes a bug in the Fortran code - if n = 0, zeros
    # are returned instead of the correct value.
    a = prolate_swf.profcn(c=c, m=m, lnum=n-m+2, x1=xi-1.0,
                           ioprad=ioprad, iopang=0, iopnorm=0, arg=[0])
    p = swf_t._make(a)

    if ioprad == 1:
        s = np.inf
        sp = np.inf
    else:
        s = p.r2c * np.float_power(10.0, p.ir2e)
        sp = p.r2dc * np.float_power(10.0, p.ir2de)

    return s[n-m], sp[n-m]

echosms.shape_conversions

Functions that convert between different echoSMs datastore shape representations.

surface_to_outline(shape, slice_thickness=0.005)

Convert a surface shape to an outline shape.

PARAMETER DESCRIPTION
shape

An echoSMs surface shape.

TYPE: dict

slice_thickness

The slice thickness [m] used when generating the outline.

TYPE: float DEFAULT: 0.005

RETURNS DESCRIPTION
dict

An echoSMs outline shape with shape metadata as per the input shape.
Notes

The conversion projects the surface shape to get dorsal and lateral outlines and then slices along the x-axis at a configurable resolution to produce the outline shape.

Source code in src/echosms/shape_conversions.py 
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
def surface_to_outline(shape: dict, slice_thickness: float=5e-3) -> dict:
    """Convert a surface shape to an outline shape.

    Parameters
    ----------
    shape :
        An echoSMs surface shape.
    slice_thickness :
        The slice thickness [m] used when generating the outline.

    Returns
    -------
    :
        An echoSMs outline shape with shape metadata as per the input shape.

    Notes
    -----
    The conversion projects the surface shape to get dorsal and lateral outlines and then
    slices along the _x_-axis at a configurable resolution to produce the outline shape.
    """

    # Put the shape into a trimesh mesh
    v = np.array([shape['x'], shape['y'], shape['z']]).T
    f = np.array([shape['facets_0'], shape['facets_1'], shape['facets_2']]).T
    mesh = trimesh.Trimesh(vertices=v, faces=f)

    # Project the surface mesh onto dorsal and lateral planes
    dorsal = projected(mesh, normal=[0, 0, 1], ignore_sign=True)
    lateral = projected(mesh, normal=[0, -1, 0], ignore_sign=True)

    # Get bounds for a centreline on the x-axis that extends the full length of the organism
    bounds = mesh.bounding_box
    xmin = bounds.vertices[0, 0]
    xmax = bounds.vertices[7, 0]

    # calculate the shape heights, widths, and y and z coordinates of the centreline
    widths = []
    heights = []
    centreline_x = []
    centreline_y = []
    centreline_z = []

    for x in np.arange(xmin, xmax, slice_thickness):
        # Get the points where a line perpendicular to the x-axis intersects
        # the dorsal and lateral shapes

        # A line perpendicular to the x-axis that extends off a long way (1000 m)
        line = LineString([[x, -1000], [x, 1000]])
        # and the intersection of that line with the dorsal shape
        intersect = intersection(dorsal, line)

        # If there is no intersection, go to the next x-point. This can happen
        # at the start or end of the bounding box
        if not intersect:
            continue

        # The length of that line is the width of the shape at this x position
        w = intersect.length
        # and the centre point of that line is the y coordinate of the centreline
        centre = intersect.interpolate(0.5, normalized=True)
        y = -centre.y

        # Do similar for the lateral outline
        line = LineString([[-1000, x], [1000, x]])
        intersect = intersection(lateral, line)
        if not intersect:
            continue

        heights.append(intersect.length)
        centre = intersect.interpolate(0.5, normalized=True)

        widths.append(w)
        centreline_x.append(x)
        centreline_y.append(y)
        centreline_z.append(centre.x)

    # Create an echoSMs shape dict using the metadata from the input surface shape
    to_remove = ['x', 'y', 'z', 'facets_0', 'facets_1', 'facets_2',
                 'normals_x', 'normals_y', 'normals_z']
    outline_shape = {k: v for k, v in shape.items() if k not in to_remove}

    # Add the outline shape data
    outline_shape['x'] = centreline_x
    outline_shape['y'] = centreline_y
    outline_shape['z'] = centreline_z
    outline_shape['height'] = heights
    outline_shape['width'] = widths

    return outline_shape

Anatomical datastore

echosms.utils_datastore

Utilities for working with the echoSMs anatomical datastore.

dwbaorganism_from_datastore(shape)

Create a DWBAorganism instance from an echoSMs datastore shape.

Converts the centreline and width/height definition of a shape into that required by the echoSMs implementation of the DWBA (centreline, tangential, and radii vectors).

PARAMETER DESCRIPTION
shape

The shape to convert, in the echoSMs datastore outline shape data structure.

TYPE: dict

RETURNS DESCRIPTION
An instance of DWBAorganism
Notes

The DWBA simulates a circular shape but the echoSMs datastore shape can store non- circular shapes (via the height and width). This function uses the height information and ignores the width information.

If mass_density_ratio and sound_speed_ratio are present into the shape dict, these are used. If not, default values are used by DWBorganism().

Source code in src/echosms/utils_datastore.py 
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
def dwbaorganism_from_datastore(shape: dict):
    """Create a DWBAorganism instance from an echoSMs datastore shape.

    Converts the centreline and width/height definition of a shape into that
    required by the echoSMs implementation of the DWBA (centreline, tangential, and
    radii vectors).

    Parameters
    ----------
    shape :
        The shape to convert, in the echoSMs datastore `outline` shape data structure.

    Returns
    -------
        An instance of DWBAorganism

    Notes
    -----
    The DWBA simulates a circular shape but the echoSMs datastore shape can store non-
    circular shapes (via the height and width). This function uses the height information
    and ignores the width information.

    If `mass_density_ratio` and `sound_speed_ratio` are present into the shape dict,
    these are used. If not, default values are used by DWBorganism().
    """
    from echosms import create_dwba_from_xyza  # here to avoid a circular import
    a = np.array(shape['height']) * 0.5
    if 'mass_density_ratio' in shape and 'sound_speed_ratio' in shape:
        return create_dwba_from_xyza(shape['x'], shape['y'], shape['z'], a,
                                     shape['name'], shape['mass_density_ratio'],
                                     shape['sound_speed_ratio'])

    return create_dwba_from_xyza(shape['x'], shape['y'], shape['z'], a, shape['name'])

krmorganism_from_datastore(shapes)

Create a KRMorganism instance from an echoSMs datastore shape.

Converts the centreline and width/height definition of a shape into that required by the echoSMs implementation of the KRM (straight centreline, width, upper and lower heights from the centreline).

PARAMETER DESCRIPTION
shapes

The shapes to convert, in the echoSMs datastore outline shape data structure.

TYPE: list[dict]

RETURNS DESCRIPTION
Instances of KRMorganism
Notes

The shape with name body becomes the main organism body and all other shapes become inclusions. If there is no shape with name of body, the first shape is used for the body.

The KRM uses just one sound speed and density per shape, but datastore shapes can have values per x-axis value. The mean of the sound speed and density values are used if so.

Datastore shapes can have non-zero y-axis values but these are ignored when creating a KRMorganism instance.

Source code in src/echosms/utils_datastore.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
def krmorganism_from_datastore(shapes: list[dict]) -> list:
    """Create a KRMorganism instance from an echoSMs datastore shape.

    Converts the centreline and width/height definition of a shape into that
    required by the echoSMs implementation of the KRM (straight centreline, width, upper and
    lower heights from the centreline).

    Parameters
    ----------
    shapes :
        The shapes to convert, in the echoSMs datastore `outline` shape data structure.

    Returns
    -------
        Instances of KRMorganism

    Notes
    -----
    The shape with name `body` becomes the main organism body and all other shapes become
    inclusions. If there is no shape with name of `body`, the first shape is used for the body.

    The KRM uses just one sound speed and density per shape, but datastore shapes can have values
    per _x_-axis value. The mean of the sound speed and density values are used if so.

    Datastore shapes can have non-zero _y_-axis values but these are ignored when creating
    a KRMorganism instance.

    """
    from echosms import KRMorganism, KRMshape  # here to avoid a circular import

    def _to_KRMshape(s: dict):
        """Convert echoSMs datstore shape into a KRMshape."""
        # Take mean of sound speed and density in case there is more than one value.
        c = sum(s['sound_speed_compressional'])/len(s['sound_speed_compressional'])
        rho = sum(s['mass_density'])/len(s['mass_density'])

        height2 = np.array(s['height'])/2.0
        return KRMshape(s['boundary'], np.array(s['x']), np.array(s['width']),
                        s['z'] + height2, s['z'] - height2, c, rho)

    if len(shapes) == 0:
        return KRMorganism('', '', [], [])

    KRMshapes = [_to_KRMshape(s) for s in shapes]

    # get the index of the first shape with name == 'body' (if any)
    idx = [i for i, s in enumerate(shapes) if s['name'] == 'body']
    if not idx:
        idx = [0]  # No shape with name of body so we use the first shape as the body

    body = KRMshapes.pop(idx[0])
    inclusions = KRMshapes

    return KRMorganism('', '', body, inclusions)

mesh_from_datastore(shapes)

Create trimesh instances from an echoSMs datastore surface shape.

PARAMETER DESCRIPTION
shapes

The shapes to convert, in the echoSMs datastore surface shape data structure.

TYPE: list[dict]

RETURNS DESCRIPTION
The shapes in trimesh form, in the same order as the input.
Source code in src/echosms/utils_datastore.py 
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
def mesh_from_datastore(shapes: list[dict]) -> list[trimesh]:
    """Create trimesh instances from an echoSMs datastore surface shape.

    Parameters
    ----------
    shapes :
        The shapes to convert, in the echoSMs datastore `surface` shape data structure.

    Returns
    -------
        The shapes in trimesh form, in the same order as the input.

    """

    def _to_trimesh(s: dict) -> trimesh.Trimesh:
        """Put echoSMs datstore shape into a trimesh instance."""
        faces = [f for f in zip(s['facets_0'], s['facets_1'], s['facets_2'])]
        vertices = [v for v in zip(s['x'], s['y'], s['z'])]

        return trimesh.Trimesh(vertices=vertices, faces=faces, process=False)

    return [_to_trimesh(s) for s in shapes]

outline_from_dwba(x, z, radius, anatomical_type='body', boundary='pressure-release')

Convert DWBA shape to the echoSMs outline shape representation.

PARAMETER DESCRIPTION
x

The x values of the centreline
z

The distance of the centreline from z = 0. Positive values are towards the dorsal surface and negative values towards the ventral surface.
radius

The radius of the shape at each x coordinate
anatomical_type

The anatomical type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'body'

boundary

The boundary type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'pressure-release'

RETURNS DESCRIPTION
An echoSMs outline shape representation.
Source code in src/echosms/utils_datastore.py 
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def outline_from_dwba(x, z, radius, anatomical_type: str = "body",
                      boundary: str = 'pressure-release') -> dict:
    """
    Convert DWBA shape to the echoSMs outline shape representation.

    Parameters
    ----------
    x :
        The _x_ values of the centreline
    z :
        The distance of the centreline from _z_ = 0. Positive values are towards
        the dorsal surface and negative values towards the ventral surface.
    radius :
        The radius of the shape at each _x_ coordinate
    anatomical_type :
        The anatomical type for this shape, as per the echoSMs datastore schema.
    boundary :
        The boundary type for this shape, as per the echoSMs datastore schema.

    Returns
    -------
     An echoSMs outline shape representation.

    """
    return {'anatomical_type': anatomical_type,
            'boundary': boundary,
            'shape_units': 'm',
            'x': np.array(x).tolist(),
            'y': np.zeros(len(x)).tolist(),
            'z': (-np.array(z)).tolist(),
            'height': (2*np.array(radius)).tolist(),
            'width': (2*np.array(radius)).tolist()}

outline_from_krm(x, height_u, height_l, width, anatomical_type='body', boundary='pressure-release')

Convert KRM shape representation to the echoSMs outline shape representation.

PARAMETER DESCRIPTION
x

The x values of the centreline

TYPE: ArrayLike

height_u

The distance from z = 0 to the upper part of the shape at each x coordinate. Positive values are towards the dorsal surface and negative values towards the ventral surface.

TYPE: ArrayLike

height_l

The distance from z = 0 to the lower part of the shape at each x coordinate. Positive values are towards the dorsal surface and negative values towards the ventral surface.

TYPE: ArrayLike

width

The width of the shape at each x coordinate

TYPE: ArrayLike

anatomical_type

The anatomical type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'body'

boundary

The boundary type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'pressure-release'

RETURNS DESCRIPTION
An echoSMs outline shape representation.
Source code in src/echosms/utils_datastore.py 
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
def outline_from_krm(x: npt.ArrayLike, height_u: npt.ArrayLike, height_l: npt.ArrayLike,
                     width: npt.ArrayLike,
                     anatomical_type: str = "body",
                     boundary: str = 'pressure-release') -> dict:
    """
    Convert KRM shape representation to the echoSMs outline shape representation.

    Parameters
    ----------
    x :
        The _x_ values of the centreline
    height_u :
        The distance from _z_ = 0 to the upper part of the shape at each _x_ coordinate.
        Positive values are towards the dorsal surface and negative values towards the ventral
        surface.
    height_l :
        The distance from _z_ = 0 to the lower part of the shape at each _x_ coordinate.
        Positive values are towards the dorsal surface and negative values towards the ventral
        surface.
    width :
        The width of the shape at each _x_ coordinate
    anatomical_type :
        The anatomical type for this shape, as per the echoSMs datastore schema.
    boundary :
        The boundary type for this shape, as per the echoSMs datastore schema.

    Returns
    -------
     An echoSMs outline shape representation.
    """
    y = np.zeros(len(x))
    height = np.array(height_u) - np.array(height_l)
    z = -(np.array(height_l) + height / 2.0)

    return {'anatomical_type': anatomical_type, 'boundary': boundary,
            'shape_units': 'm',
            'x': np.array(x).tolist(),
            'y': y.tolist(),
            'z': z.tolist(),
            'height': height.tolist(),
            'width': np.array(width).tolist()}

plot_shape_categorised_voxels(s, title='')

Plot the specimen's categorised voxels.

Normally called via plot_specimen().

PARAMETER DESCRIPTION
s

The categorised voxel shape data structure as per the echoSMs datastore.
title

Title for the plot.

DEFAULT: ''

Source code in src/echosms/utils_datastore.py 
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
def plot_shape_categorised_voxels(s, title=''):
    """Plot the specimen's categorised voxels.

    Normally called via [plot_specimen()][echosms.utils_datastore.plot_specimen].

    Parameters
    ----------
    s :
        The categorised voxel shape data structure as per the echoSMs datastore.
    title :
        Title for the plot.
    """
    d = np.array(s['categories'])
    voxel_size = np.array(s['voxel_size'])
    shape = d.shape

    cats = np.unique(d)
    norm = colors.Normalize(vmin=min(cats), vmax=max(cats))

    row_dim = np.linspace(0, voxel_size[0]*1e3*shape[0], shape[0]+1)
    slice_dim = np.linspace(0, voxel_size[2]*1e3*shape[2], shape[2]+1)

    cmap = colormaps['Dark2']

    # Create 25 plots along the organism's echoSMs x-axis
    fig, axs = plt.subplots(5, 5, sharex=True, sharey=True)
    cols = np.linspace(0, shape[1]-1, num=25)

    for col, ax in zip(cols, axs.flat):
        c = int(floor(col))
        # The [::-1] and .invert_ axis calls give the appropriate
        # x and y axes directions in the plots.
        ax.pcolormesh(slice_dim[::-1], row_dim[::-1], d[:,c,:],
                      norm=norm, cmap=cmap)

        ax.set_aspect('equal')
        ax.invert_xaxis()
        ax.invert_yaxis()

        ax.text(0.05, .86, f'{col*1e3*voxel_size[1]:.0f} mm',
                transform=ax.transAxes, fontsize=6, color='white')

    fig.supxlabel('y [mm]')
    fig.supylabel('z [mm]')
    fig.suptitle(title)

plot_shape_outline(shapes, axs)

Plot an echoSMs anatomical outline shape.

Normally called via plot_specimen().

PARAMETER DESCRIPTION
shapes

Outline shapes to be plotted

TYPE: list[dict]

axs

Two matplotlib axes - one for the dorsal view and one for the lateral view

TYPE: list

Source code in src/echosms/utils_datastore.py 
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
def plot_shape_outline(shapes: list[dict], axs: list) -> None:
    """Plot an echoSMs anatomical outline shape.

    Normally called via [plot_specimen()][echosms.utils_datastore.plot_specimen].

    Parameters
    ----------
    shapes :
        Outline shapes to be plotted
    axs :
        Two matplotlib axes - one for the dorsal view and one for the
        lateral view
    """
    for s in shapes:
        c = 'C0' if s['boundary'] == bt.fluid_filled else 'C1'
        x = np.array(s['x'])*1e3
        z = np.array(s['z'])*1e3
        y = np.array(s['y'])*1e3
        width_2 = np.array(s['width'])/2*1e3
        zU = (z + np.array(s['height'])/2*1e3)
        zL = (z - np.array(s['height'])/2*1e3)

        # Dorsal view
        axs[0].plot(x, y, c='grey', linestyle='--', linewidth=1)  # centreline
        axs[0].plot(x, y+width_2, c=c)
        axs[0].plot(x, y-width_2, c=c)

        # Lateral view
        axs[1].plot(x, z, c='grey', linestyle='--', linewidth=1)  # centreline
        axs[1].plot(x, zU, c=c)
        axs[1].plot(x, zL, c=c)

        # close the ends of the shapes
        for i in [0, -1]:
            axs[1].plot([x[i], x[i]], [zU[i], zL[i]], c=c)
            axs[0].plot([x[i], x[i]], [(y+width_2)[i], (y-width_2)[i]], c=c)
            axs[i].xaxis.set_inverted(True)
            axs[i].yaxis.set_inverted(True)

plot_shape_surface(shapes, ax)

Plot an echoSMs anatomical surface shape.

Normally called via plot_specimen().

PARAMETER DESCRIPTION
shapes

Surface shapes to be plotted
ax

A matplotlib axis.
Source code in src/echosms/utils_datastore.py 
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
def plot_shape_surface(shapes, ax):
    """Plot an echoSMs anatomical surface shape.

    Normally called via [plot_specimen()][echosms.utils_datastore.plot_specimen].

    Parameters
    ----------
    shapes :
        Surface shapes to be plotted
    ax :
        A matplotlib axis.
    """
    for s in shapes:
        # c = 'C0' if s['boundary'] == bt.fluid_filled else 'C1'
        facets = np.array([s['facets_0'], s['facets_1'], s['facets_2']]).transpose()
        x = 1e3 * np.array(s['x'])
        y = 1e3 * np.array(s['y'])
        z = 1e3 * np.array(s['z'])

        ax.plot_trisurf(x, y, z, triangles=facets)
        ax.view_init(elev=210, azim=-60, roll=0)
        ax.set_xlabel('x')
        ax.set_ylabel('y')
        ax.set_zlabel('z')

        ax.set_aspect('equal')
        ax.xaxis.set_inverted(True)
        ax.yaxis.set_inverted(True)

plot_shape_voxels(s, title='')

Plot the specimen's voxels.

Normally called via plot_specimen().

PARAMETER DESCRIPTION
s

The voxel shape data structure as per the echoSMs datastore.
title

Title for the plot.

DEFAULT: ''

Source code in src/echosms/utils_datastore.py 
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
def plot_shape_voxels(s, title=''):
    """Plot the specimen's voxels.

    Normally called via [plot_specimen()][echosms.utils_datastore.plot_specimen].

    Parameters
    ----------
    s :
        The voxel shape data structure as per the echoSMs datastore.

    title :
        Title for the plot.

    """
    # Show density. Could do sound speed or some impedance proxy.
    d = np.array(s['sound_speed_compressional'])
    voxel_size = np.array(s['voxel_size'])
    shape = d.shape

    # Make the colours ignore extreme high value and the lowest low values
    norm = colors.Normalize(vmin=np.percentile(d.flat, 1),
                            vmax=np.percentile(d.flat, 99))

    row_dim = np.linspace(0, voxel_size[0]*1e3*shape[0], shape[0]+1)
    slice_dim = np.linspace(0, voxel_size[2]*1e3*shape[2], shape[2]+1)

    cmap = colormaps['viridis']

    # Create 25 plots along the organism's echoSMs x-axis
    fig, axs = plt.subplots(5, 5, sharex=True, sharey=True)
    cols = np.linspace(0, shape[1]-1, num=25)

    for col, ax in zip(cols, axs.flat):
        c = int(floor(col))
        # The [::-1] and .invert_ axis calls give the appropriate
        # x and y axes directions in the plots.
        im = ax.pcolormesh(slice_dim[::-1], row_dim[::-1], d[:,c,:],
                           norm=norm, cmap=cmap)

        ax.set_aspect('equal')
        ax.invert_xaxis()
        ax.invert_yaxis()

        ax.text(0.05, .86, f'{col*1e3*voxel_size[1]:.0f} mm',
                transform=ax.transAxes, fontsize=6, color='white')

    # A single colorbar in the plot
    cbar = fig.colorbar(im, ax=axs, orientation='vertical',
                        fraction=0.1, extend='both', cmap=cmap)
    cbar.ax.set_ylabel('[kg m$^{-3}$]')
    fig.supxlabel('y [mm]')
    fig.supylabel('z [mm]')
    fig.suptitle(title)

plot_specimen(specimen, dataset_id='', title='', savefile=None, dpi=150)

Plot the specimen shape.

Produces a relevant plot for all echoSMs anatomical datastore shape types.

PARAMETER DESCRIPTION
specimen

Specimen data as per the echoSMs anatomical datastore schema.

TYPE: dict

dataset_id

Used to form a plot title if title is an empty string.

TYPE: str DEFAULT: ''

title

A title for the plot.

TYPE: str DEFAULT: ''

savefile

Filename to save the plot to. If None, generate the plot in the interactive terminal (if that's supported).

TYPE: str | None DEFAULT: None

dpi

The resolution of the figure in dots per inch.

TYPE: float DEFAULT: 150

Source code in src/echosms/utils_datastore.py 
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
def plot_specimen(specimen: dict, dataset_id: str='', title: str='',
                  savefile: str|None=None, dpi: float=150) -> None:
    """Plot the specimen shape.

    Produces a relevant plot for all echoSMs anatomical datastore shape types.

    Parameters
    ----------
    specimen :
        Specimen data as per the echoSMs anatomical datastore schema.
    dataset_id :
        Used to form a plot title if `title` is an empty string.
    title :
        A title for the plot.
    savefile :
        Filename to save the plot to. If None, generate the plot in the
        interactive terminal (if that's supported).
    dpi :
        The resolution of the figure in dots per inch.

    """
    labels = ['Dorsal', 'Lateral']
    t = title if title else dataset_id + ' ' + specimen['specimen_id']

    match specimen['shape_type']:
        case 'outline':
            fig, axs = plt.subplots(2, 1, sharex=True, layout='tight', dpi=dpi)
            fig.set_tight_layout({'h_pad': 1, 'w_pad': 1})
            plot_shape_outline(specimen['shapes'], axs)
            for label, a in zip(labels, axs):
                a.set_title(label, loc='left', fontsize=8)
                a.axis('scaled')
            axs[0].set_title(t)
            _fit_to_axes(fig)

        case 'surface':
            fig, ax = plt.subplots(subplot_kw={'projection': '3d'})
            plot_shape_surface(specimen['shapes'], ax)
            plt.tight_layout()
            ax.set_title(t)
        case 'voxels':
            plot_shape_voxels(specimen['shapes'][0], t)
        case 'categorised voxels':
            plot_shape_categorised_voxels(specimen['shapes'][0], t)

    if savefile:
        plt.savefig(savefile, format='png', dpi=dpi)
        plt.close()
    else:
        plt.show()

surface_from_stl(stl_file, dim_scale=1.0, anatomical_type='body', boundary='pressure-release')

Create an echoSMs surface shape from an .stl file.

PARAMETER DESCRIPTION
stl_file

An .stl file

TYPE: str | Path

dim_scale

Scaling factor applied to the node positions. Use to convert from one length unit to another (e.g., 1e-3 will convert from mm to m).

TYPE: float DEFAULT: 1.0

anatomical_type

The anatomical type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'body'

boundary

The boundary type for this shape, as per the echoSMs datastore schema.

TYPE: str DEFAULT: 'pressure-release'

RETURNS DESCRIPTION
dict

An echoSMs surface shape representation.
Notes

This function uses a call to load_mesh() from the trimesh library to read the .stl file. If there are problems with loading your .stl file, please refer to the trimesh documentation.

Source code in src/echosms/utils_datastore.py 
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
def surface_from_stl(stl_file: str | Path,
                     dim_scale: float = 1.0,
                     anatomical_type: str = 'body',
                     boundary: str = 'pressure-release') -> dict:
    """Create an echoSMs surface shape from an .stl file.

    Parameters
    ----------
    stl_file :
        An .stl file
    dim_scale :
        Scaling factor applied to the node positions. Use to convert from one
        length unit to another (e.g., 1e-3 will convert from mm to m).
    anatomical_type :
        The anatomical type for this shape, as per the echoSMs datastore schema.
    boundary :
        The boundary type for this shape, as per the echoSMs datastore schema.

    Returns
    -------
    :
        An echoSMs surface shape representation.

    Notes
    -----
    This function uses a call to `load_mesh()` from the `trimesh` library to read the
    .stl file. If there are problems with loading your .stl file, please refer to the
    `trimesh` documentation.
    """
    mesh = trimesh.load_mesh(stl_file)

    # Bundle up into a dict as per the echoSMs schema for a surface
    return {'anatomical_type': anatomical_type, 'boundary': boundary,
            'shape_units': 'm',
            'x': (mesh.vertices[:, 0]*dim_scale).tolist(),
            'y': (mesh.vertices[:, 1]*dim_scale).tolist(),
            'z': (mesh.vertices[:, 2]*dim_scale).tolist(),
            'facets_0': mesh.faces[:, 0].tolist(),
            'facets_1': mesh.faces[:, 1].tolist(),
            'facets_2': mesh.faces[:, 2].tolist(),
            'normals_x': mesh.face_normals[:, 0].tolist(),
            'normals_y': mesh.face_normals[:, 1].tolist(),
            'normals_z': mesh.face_normals[:, 2].tolist()}

volume_from_datastore(voxels)

Create a 3D numpy array from nested lists.

PARAMETER DESCRIPTION
voxels

The datastore 3D voxel structure (list of list of list)

TYPE: list

RETURNS DESCRIPTION
A numpy 3D array.
Source code in src/echosms/utils_datastore.py 
126
127
128
129
130
131
132
133
134
135
136
137
138
def volume_from_datastore(voxels: list):
    """Create a 3D numpy array from nested lists.

    Parameters
    ----------
    voxels :
        The datastore 3D voxel structure (list of list of list)

    Returns
    -------
        A numpy 3D array.
    """
    return np.array(voxels)  # TODO - check ordering is correct!