Skip to content

Utilities

This section provides documentation on the various echoSMs utility functions and classes. Miscellaneous utility code

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.

outline_to_surface(outline, num_pts=20, mesh_len=2.0)

Convert an outline shape to a surface shape.

PARAMETER DESCRIPTION
outline

An echoSMs outline shape.

TYPE: dict

num_pts

The number of points to place on each cross-sectional ellipse.

TYPE: int DEFAULT: 20

mesh_len

The desired typical mesh length as a percentage of overall object size

TYPE: float DEFAULT: 2.0

RETURNS DESCRIPTION
dict[str, list]

An echoSMs surface shape with shape metadata as per the input shape.
Notes

Each outline cross-sectional ellipse is represented by a polygon with num_pts vertices. Triangles are created that join the vertices on adjacent polygons. The two ends are meshed using an ear slicing algorithm (using the mapbox_earcut package, a Python binding to a C++ implementation of the algorithm).

The mesh is then remeshed using the pymeshlab isotropic remeshing algorithm.

Source code in src/echosms/shape_conversions.py 
 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
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
def outline_to_surface(outline: dict, num_pts:int = 20, mesh_len:float = 2.0) -> dict:
    """Convert an outline shape to a surface shape.

    Parameters
    ----------
    outline :
        An echoSMs outline shape.
    num_pts :
        The number of points to place on each cross-sectional ellipse.
    mesh_len :
        The desired typical mesh length as a percentage of overall object size

    Returns
    -------
    : dict[str, list]
        An echoSMs surface shape with shape metadata as per the input shape.

    Notes
    -----
    Each outline cross-sectional ellipse is represented by a polygon with num_pts
    vertices. Triangles are created that join the vertices on adjacent polygons.
    The two ends are meshed using an ear slicing algorithm (using the `mapbox_earcut` package, a
    Python binding to a [C++ implementation](https://github.com/mapbox/earcut.hpp) of the
    algorithm).

    The mesh is then remeshed using the pymeshlab isotropic remeshing algorithm.
    """
    num_discs = len(outline['x'])

    # Create points around each ellipse cross-section of the outline shape
    t = np.linspace(0, 2*np.pi, num=num_pts, endpoint=False)
    pts = []
    # Could vectorise this, but then the code is harder to understand and the number of
    # discs that this will iterate over is fairly small so speed isn't a concern
    for i in range(num_discs):
        pts_y = outline['y'][i] + outline['width'][i]/2 * np.cos(t)
        pts_z = outline['z'][i] + outline['height'][i]/2 * np.sin(t)
        pts_x = np.full(pts_y.shape, outline['x'][i])
        pts.extend(np.c_[pts_x, pts_y, pts_z].tolist())

    # Create triangles connecting respective points on each ellipse
    # Same vectorisation/speed comment here as above
    faces = []
    for disc in range(num_discs-1):
        for pt_i in range(num_pts):
            face = [disc*num_pts + pt_i,
                    (disc+1)*num_pts + pt_i,
                    disc*num_pts + (pt_i+1) % num_pts]
            faces.append(face)

            face = [disc*num_pts + (pt_i+1) % num_pts,
                    (disc+1)*num_pts + pt_i,
                    (disc+1)*num_pts + (pt_i+1) % num_pts]
            faces.append(face)

    # Create triangles for the two end surfaces
    pts2d = [[p[1], p[2]] for p in pts]  # shapely.Polygon wants a 2D polygon, so remove the x coord

    _, endcap1_faces = triangulate_polygon(Polygon(pts2d[:num_pts]), engine='earcut')
    # the order of the nodes is inverted for endcap2 to have the normals point outwards
    _, endcap2_faces = triangulate_polygon(Polygon(pts2d[:-(num_pts+1):-1]), engine='earcut')
    # Get the right facet indices for endcap2
    endcap2_faces = [f - 1 + num_pts * (num_discs-1) for f in endcap2_faces]

    faces.extend(endcap1_faces)
    faces.extend(endcap2_faces)

    # Tidy the mesh using pymeshlab
    ms = pymeshlab.MeshSet()
    ms.add_mesh(pymeshlab.Mesh(pts, faces))
    # ms.save_current_mesh('test_before.stl')
    #ms.meshing_merge_close_vertices()
    #ms.meshing_remove_t_vertices()
    #ms.meshing_close_holes()
    #ms.meshing_repair_non_manifold_edges()
    #ms.meshing_re_orient_faces_coherently()
    #ms.meshing_isotropic_explicit_remeshing(targetlen=pymeshlab.PercentageValue(mesh_len),
    #                                        adaptive=True)
    # ms.save_current_mesh('test_after.stl')
    m = ms.current_mesh()

    # Put into trimesh to get the face normals and do some checks
    mesh = Trimesh(vertices=m.vertex_matrix(), faces=m.face_matrix())

    errors = []
    if not mesh.is_watertight:
        errors.append('Mesh is not watertight')
    if not mesh.is_winding_consistent:
        errors.append('Mesh winding is not consistent')
    if errors:
        raise ValueError(', '.join(errors))

    # structure as an echoSMs surface dict
    surface = {'x': mesh.vertices[:, 0].tolist(),
               'y': mesh.vertices[:, 1].tolist(),
               'z': mesh.vertices[:, 2].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(),
               }

    # Copy across other attributes from the outline shape
    attrs = {k:v for k, v in outline.items() if k not in ['x', 'y', 'z', 'height', 'width']}

    return attrs | surface

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 
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
195
196
197
198
199
200
201
202
203
204
205
206
207
208
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