diff --git a/echopype/calibrate/__init__.py b/echopype/calibrate/__init__.py index 94613c73b..8fb4e732e 100644 --- a/echopype/calibrate/__init__.py +++ b/echopype/calibrate/__init__.py @@ -1,3 +1,3 @@ -from .api import compute_Sv, compute_TS +from .api import compute_Sp, compute_Sv, compute_Sv_spectrum, compute_TS, compute_TS_spectrum -__all__ = ["compute_Sv", "compute_TS"] +__all__ = ["compute_Sv", "compute_TS", "compute_Sp", "compute_Sv_spectrum", "compute_TS_spectrum"] diff --git a/echopype/calibrate/api.py b/echopype/calibrate/api.py index d12b27b7c..058b4e96e 100644 --- a/echopype/calibrate/api.py +++ b/echopype/calibrate/api.py @@ -1,3 +1,5 @@ +import warnings + import numpy as np import xarray as xr @@ -21,7 +23,7 @@ def _compute_cal( - cal_type, + cal_type: str, echodata: EchoData, env_params=None, cal_params=None, @@ -30,9 +32,20 @@ def _compute_cal( encode_mode=None, assume_single_filter_time=None, drop_last_hanning_zero=False, + **kwargs, ): - # Make waveform_mode "FM" equivalent to "BB" - waveform_mode = "BB" if waveform_mode == "FM" else waveform_mode + # Make waveform_mode "FM" equivalent to "BB". + # Accept legacy "BB" for backward compatibility. + # Ref: https://github.com/echostack-org/echopype/issues/1651 + if waveform_mode == "BB": + warnings.warn( + "'BB' is deprecated and will be removed in a future release. " + "Please use 'FM' instead.", + DeprecationWarning, + stacklevel=2, + ) + + waveform_mode = "BB" if waveform_mode in ("FM", "BB") else waveform_mode # TODO: consolidate the below block with simrad.py::check_input_args_combination() # Check on waveform_mode, encode_mode inputs, and assumption on single filter time @@ -80,13 +93,27 @@ def _compute_cal_ds(echodata, slice_dict): # Check Echodata backscatter data size and recommend chunking if data is too large cal_obj._check_echodata_backscatter_size() - # Perform calibration - if cal_type == "Sv": - cal_ds = cal_obj.compute_Sv() - else: - cal_ds = cal_obj.compute_TS() + compute_methods = { + "Sp": "compute_Sp", + "TS": "compute_TS", + "Sv": "compute_Sv", + # add Sp_spectrum?? + "TS_spectrum": "compute_TS_spectrum", + } + + try: + method_name = compute_methods[cal_type] + except KeyError: + raise ValueError(f"Unsupported calibration type: {cal_type}") from None + + compute_method = getattr(cal_obj, method_name, None) + + if compute_method is None: + raise ValueError( + f"{cal_type} calibration is not supported for " f"{echodata.sonar_model} data." + ) - return cal_ds + return compute_method(**kwargs) # Calibrate as a single dataset if not Ex80 if echodata.sonar_model not in ["EK80", "ES80", "EA640"]: @@ -161,7 +188,10 @@ def _compute_cal_ds(echodata, slice_dict): # Calibrate and drop filter_time cal_ds_iteration = _compute_cal_ds(echodata, slice_dict) - cal_ds_list.append(cal_ds_iteration.drop_vars("filter_time")) + if "filter_time" in cal_ds_iteration: + cal_ds_iteration = cal_ds_iteration.drop_vars("filter_time") + + cal_ds_list.append(cal_ds_iteration) # # Alternative? # for channel in echodata[ed_beam_group]["channel"].values: @@ -201,12 +231,27 @@ def _add_attrs(cal_type, ds): """Add attributes to backscattering strength dataset. cal_type: Sv or TS """ - ds["range_sample"].attrs = {"long_name": "Along-range sample number, base 0"} - ds["echo_range"].attrs = {"long_name": "Range distance", "units": "m"} + if "range_sample" in ds: + ds["range_sample"].attrs = {"long_name": "Along-range sample number, base 0"} + + if "echo_range" in ds: + ds["echo_range"].attrs = { + "long_name": "Range distance", + "units": "m", + } + + if "frequency" in ds: + ds["frequency"].attrs = { + "long_name": "Frequency", + "units": "Hz", + } + ds[cal_type].attrs = { "long_name": { - "Sv": "Volume backscattering strength (Sv re 1 m-1)", + "Sp": "Point scattering strength (Sp re 1 m^2)", "TS": "Target strength (TS re 1 m^2)", + "Sv": "Volume backscattering strength (Sv re 1 m-1)", + "TS_spectrum": "Frequency-dependent target strength spectrum (TS(f) re 1 m^2)", }[cal_type], "units": "dB", } @@ -345,7 +390,133 @@ def compute_Sv(echodata: EchoData, **kwargs) -> xr.Dataset: return _compute_cal(cal_type="Sv", echodata=echodata, **kwargs) -def compute_TS(echodata: EchoData, **kwargs): +def compute_Sv_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute frequency-dependent volume backscattering strength Sv(f) + from broadband EK80 complex data. + + Notes + ----- + This functionality is not yet implemented. + """ + raise NotImplementedError("compute_Sv_spectrum is not yet implemented.") + + +def compute_Sp(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute point scattering strength (Sp) from raw data. + + For CW data, Sp is computed from received power samples on the range grid. + For EK80 broadband/FM complex data, Sp is computed after pulse compression + and represents a band-averaged point-scattering-strength echogram. + """ + return _compute_cal(cal_type="Sp", echodata=echodata, **kwargs) + + +def _compute_TS_from_Sp( + source_Sp: xr.Dataset, + point_locations: xr.Dataset, +) -> xr.Dataset: + """Compute single-target TS values from an Sp dataset.""" + + target_dim = "single_target" + + if target_dim not in point_locations.dims: + raise ValueError("point_locations must use the 'single_target' dimension.") + + if "channel" not in point_locations: + raise ValueError("point_locations must contain a 'channel' variable.") + + if "beam_comp_db" not in point_locations: + raise ValueError("point_locations must contain a 'beam_comp_db' variable.") + + if point_locations.sizes[target_dim] == 0: + return point_locations.assign( + uncompensated_TS=( + target_dim, + np.array([], dtype=np.float64), + ), + compensated_TS=( + target_dim, + np.array([], dtype=np.float64), + ), + ) + + target_dsets = [] + + for channel in np.unique(point_locations["channel"].values): + targets_channel = point_locations.where( + point_locations["channel"] == channel, + drop=True, + ) + + source_channel = source_Sp.sel(channel=channel) + + ping_index = targets_channel["ping_index"].values.astype(int) + range_sample = targets_channel["range_sample"].values.astype(int) + + uncompensated_ts_raw = source_channel["Sp"].isel( + ping_time=xr.DataArray( + ping_index, + dims=target_dim, + ), + range_sample=xr.DataArray( + range_sample, + dims=target_dim, + ), + ) + + uncompensated_ts = xr.DataArray( + uncompensated_ts_raw.values, + dims=(target_dim,), + coords={ + target_dim: targets_channel[target_dim], + }, + name="uncompensated_TS", + ) + + compensated_ts = xr.DataArray( + (uncompensated_ts.values + targets_channel["beam_comp_db"].values), + dims=(target_dim,), + coords={ + target_dim: targets_channel[target_dim], + }, + name="compensated_TS", + ) + + targets_channel = targets_channel.assign( + uncompensated_TS=uncompensated_ts, + compensated_TS=compensated_ts, + ) + + target_dsets.append(targets_channel) + + result = xr.concat( + target_dsets, + dim=target_dim, + ) + + result["uncompensated_TS"].attrs = { + "long_name": ("Calculated target strength (re 1 m2) " "uncompensated for off-axis angle"), + "units": "dB", + } + + result["compensated_TS"].attrs = { + "long_name": ( + "Calculated target strength (re 1 m2) " "after compensation for off-axis angle" + ), + "units": "dB", + } + + return result + + +def compute_TS( + echodata: EchoData | xr.Dataset, + *, + point_locations: xr.Dataset | None = None, + **kwargs, +) -> xr.Dataset: """ Compute target strength (TS) from raw data. @@ -357,11 +528,15 @@ def compute_TS(echodata: EchoData, **kwargs): ---------- echodata : EchoData An `EchoData` object created by using `open_raw` or `open_converted` + point_locations : xr.Dataset, optional + + Single-target locations produced by ``detect_from_Sp``. + Required when ``echodata`` is an Sp dataset rather than an EchoData object. env_params : dict, optional Environmental parameters needed for calibration. Users can supply `"sound speed"` and `"absorption"` directly, - or specify other variables that can be used to compute them, + or specify other variables that can be used to compok what nute them, including `"temperature"`, `"salinity"`, and `"pressure"`. For EK60 and EK80 echosounders, by default echopype uses @@ -446,4 +621,62 @@ def compute_TS(echodata: EchoData, **kwargs): symbols in fisheries acoustics. ICES J. Mar. Sci. 59: 365-369. https://doi.org/10.1006/jmsc.2001.1158 """ + + if isinstance(echodata, xr.Dataset): + if point_locations is None: + raise ValueError( + "point_locations must be provided when computing " + "single-target TS from an Sp dataset." + ) + + return _compute_TS_from_Sp( + source_Sp=echodata, + point_locations=point_locations, + ) + return _compute_cal(cal_type="TS", echodata=echodata, **kwargs) + + +def compute_TS_spectrum(echodata: EchoData, **kwargs) -> xr.Dataset: + """ + Compute broadband frequency-dependent target strength spectrum, TS(f), + from EK80 broadband/FM complex data. + + Parameters + ---------- + point_locations : xr.Dataset + Locations of targets for which TS(f) should be computed. + Must contain ``channel``, ``ping_time``, and ``target_range`` for each + ``target_id``. If ``target_range_min`` and ``target_range_max`` are + provided, they define the target echo segment. Otherwise, the segment + is built around ``target_range`` using ``NFFT`` and ``split_front``. + + NFFT : int, optional + Number of FFT points used to compute the target spectrum. If not + provided, a value is inferred from the output frequency grid. + + n_f_points : int, optional + Number of frequency points in the output TS(f) spectrum. Used when + ``frequency_resolution`` is not provided. + + split_front : float, default 0.25 + Each echo spectrum is computed from a segment of the complex echo signal. + This parameter specifies how to position that segment around the target location + when only ``target_range`` is provided. For example, if ``split_front=0.25``, + then 25% of the NFFT window is placed before ``target_range`` and the remaining + 75% after it. + + window : str, tuple, float or None, default None + Window passed directly to ``scipy.signal.get_window``. If ``None``, + a rectangular/boxcar window is used. + + frequency_resolution : float, optional + Desired spacing of the output frequency grid in Hz. Used to define the + frequency grid on which TS(f) is evaluated. + + Returns + ------- + xr.Dataset + Dataset containing beam-compensated frequency-dependent target strength, TS(f). + """ + return _compute_cal(cal_type="TS_spectrum", echodata=echodata, **kwargs) diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index b58be2ca7..046b5616e 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -2,17 +2,26 @@ import numpy as np import xarray as xr +from scipy.signal import get_window from ..echodata import EchoData from ..echodata.simrad import retrieve_correct_beam_group +from ..utils import uwa from ..utils.log import _init_logger from .cal_params import _get_interp_da, get_cal_params_EK from .calibrate_base import CalibrateBase from .ecs import conform_channel_order, ecs_ds2dict, ecs_ev2ep from .ek80_complex import ( - compress_pulse, + _align_autocorrelation, + _compute_power_from_complex_signal, + _compute_ts_spectrum, + _compute_ts_spectrum_calibrated, + _compute_ts_spectrum_power, + _get_autocorrelation, + _get_average_signal, + _get_pulse_compressed_signal, + _get_splitbeam_angles, get_filter_coeff, - get_norm_fac, get_tau_effective, get_transmit_signal, ) @@ -171,17 +180,38 @@ def _cal_power_samples(self, cal_type: str) -> xr.Dataset: ) out.name = "Sv" - elif cal_type == "TS": - # Calc gain + elif cal_type in ("Sp", "TS"): CSp = ( 10 * np.log10(self.beam["transmit_power"]) + 2 * self.cal_params["gain_correction"] + 10 * np.log10(wavelength**2 / (16 * np.pi**2)) ) - # Calibration and echo integration - out = self.beam["backscatter_r"] + spreading_loss * 2 + absorption_loss - CSp - out.name = "TS" + sp = self.beam["backscatter_r"] + spreading_loss * 2 + absorption_loss - CSp + + if cal_type == "TS": + angle_alongship = self.beam["angle_alongship"] + angle_athwartship = self.beam["angle_athwartship"] + + angle_offset_alongship = self.cal_params["angle_offset_alongship"] + angle_offset_athwartship = self.cal_params["angle_offset_athwartship"] + beamwidth_alongship = self.cal_params["beamwidth_alongship"] + beamwidth_athwartship = self.cal_params["beamwidth_athwartship"] + + beam_correction_db = self._get_beam_correction( + theta=angle_alongship, + phi=angle_athwartship, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, + ) + + out = sp + beam_correction_db + else: + out = sp + + out.name = cal_type # Attach calculated range (with units meter) into data set out = out.to_dataset() @@ -262,7 +292,13 @@ def compute_Sv(self, **kwargs): return self._cal_power_samples(cal_type="Sv") def compute_TS(self, **kwargs): - return self._cal_power_samples(cal_type="TS") + raise NotImplementedError( + "Gridded TS is deprecated. Use compute_Sp() first, detect targets, " + "then call echopype.calibrate.compute_TS(source_Sp, point_locations=...)." + ) + + def compute_Sp(self, **kwargs): + return self._cal_power_samples(cal_type="Sp") class CalibrateEK80(CalibrateEK): @@ -480,25 +516,19 @@ def _get_power_from_complex( Power computed from complex samples """ - def _get_prx(sig): - return ( - beam["beam"].size # number of transducer sectors - * np.abs(sig.mean(dim="beam")) ** 2 - / (2 * np.sqrt(2)) ** 2 - * (np.abs(z_er + z_et) / z_er) ** 2 - / z_et - ) - - # Compute power if self.waveform_mode == "BB": - pc = compress_pulse( - backscatter=beam["backscatter_r"] + 1j * beam["backscatter_i"], chirp=chirp - ) # has beam dim - pc = pc / get_norm_fac(chirp=chirp) # normalization for each channel - prx = _get_prx(pc) # ensure prx is xr.DataArray + signal = _get_pulse_compressed_signal( + beam=beam, + matched_filter=chirp, + ) else: - bs_cw = beam["backscatter_r"] + 1j * beam["backscatter_i"] - prx = _get_prx(bs_cw) + signal = beam["backscatter_r"] + 1j * beam["backscatter_i"] + + prx = _compute_power_from_complex_signal( + signal=signal, + z_et=z_et, + z_er=z_er, + ) prx.name = "received_power" @@ -535,8 +565,10 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: Parameters ---------- cal_type : str - 'Sv' for calculating volume backscattering strength, or - 'TS' for calculating target strength + 'Sv' for volume backscattering strength, 'Sp' for point scattering + strength, or 'TS' for beam-compensated target strength. For BB/FM + complex data, TS is a band-averaged, center-frequency approximation; + use compute_TS_spectrum for frequency-dependent TS(f). Returns ------- @@ -626,16 +658,47 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: out.name = "Sv" # out = out.rename_vars({list(out.data_vars.keys())[0]: "Sv"}) + elif cal_type in ("Sp", "TS"): + range_safe = self._safe_range_for_log(range_meter) + spreading_loss_safe = 20 * np.log10(range_safe) - elif cal_type == "TS": - out = ( + sp = ( 10 * np.log10(prx) - + 2 * spreading_loss + + 2 * spreading_loss_safe + absorption_loss - 10 * np.log10(wavelength**2 * transmit_power / (16 * np.pi**2)) - 2 * gain ) - out.name = "TS" + + if cal_type == "TS": + # For BB/FM, use pulse-compressed sector signals for split-beam angles. + # For CW complex, use the recorded complex sector signals directly. + signal_for_angles = ( + _get_pulse_compressed_signal(beam=self.beam, matched_filter=tx) + if self.waveform_mode == "BB" + else self.beam["backscatter_r"] + 1j * self.beam["backscatter_i"] + ) + + angle_alongship, angle_athwartship = _get_splitbeam_angles( + pc=signal_for_angles, + gamma_alongship=self.cal_params["angle_sensitivity_alongship"], + gamma_athwartship=self.cal_params["angle_sensitivity_athwartship"], + ) + + beam_correction_db = self._get_beam_correction( + theta=angle_alongship, + phi=angle_athwartship, + angle_offset_alongship=self.cal_params["angle_offset_alongship"], + angle_offset_athwartship=self.cal_params["angle_offset_athwartship"], + beamwidth_alongship=self.cal_params["beamwidth_alongship"], + beamwidth_athwartship=self.cal_params["beamwidth_athwartship"], + ) + + out = sp + beam_correction_db + else: + out = sp + + out.name = cal_type # Attach calculated range (with units meter) into data set out = out.to_dataset().merge(range_meter) @@ -658,53 +721,560 @@ def _cal_complex_samples(self, cal_type: str) -> xr.Dataset: return out - def _compute_cal(self, cal_type) -> xr.Dataset: + def _select_param( + self, + da: xr.DataArray, + channel=None, + ping_idx=None, + ): + """Select channel- and ping-specific parameter values when dimensions exist.""" + if channel is not None and "channel" in da.dims: + da = da.sel(channel=channel) + if ping_idx is not None and "ping_time" in da.dims: + da = da.isel(ping_time=ping_idx) + return da + + def _safe_range_for_log(self, range_meter: xr.DataArray) -> xr.DataArray: + """Avoid log10(0) in range-dependent calibration terms.""" + return range_meter.where(range_meter > 0, 1e-20) + + def _compute_absorption_f( + self, + frequency: np.ndarray, + channel: str, + ping_idx: int, + sound_speed: float, + ) -> np.ndarray: + """Compute frequency-dependent absorption for broadband spectra. + + Uses the Francois & Garrison (1982) formulation, matching the + absorption model used by CRIMAC ``calc_alpha``. + """ + return uwa.calc_absorption( + frequency=frequency, + temperature=float( + self._select_param(self.env_params["temperature"], channel, ping_idx) + ), + salinity=float(self._select_param(self.env_params["salinity"], channel, ping_idx)), + pressure=float(self._select_param(self.env_params["pressure"], channel, ping_idx)), + pH=float(self._select_param(self.env_params["pH"], channel, ping_idx)), + sound_speed=sound_speed, + formula_source="FG", + ) + + def _get_beam_correction( + self, + theta, + phi, + angle_offset_alongship, + angle_offset_athwartship, + beamwidth_alongship, + beamwidth_athwartship, + ): + """Compute empirical split-beam beam correction in dB.""" + fac_along = (np.abs(theta - angle_offset_alongship) / (beamwidth_alongship / 2)) ** 2 + fac_athwart = (np.abs(phi - angle_offset_athwartship) / (beamwidth_athwartship / 2)) ** 2 + + return 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart) + + def _get_beam_compensated_gain( + self, + channel: str, + theta: float, + phi: float, + frequency: np.ndarray, + ) -> np.ndarray: + """Calculate beam-compensated gain. + + Equivalent to CRIMAC ``calc_g(theta, phi, f)``. + + # NOTE: + # For TS(f), theta and phi are expected to be CRIMAC-style split-beam + # angles computed directly from the pulse-compressed sector signals. + # These angles include the sensitivity conversion but do not include + # the angle-offset correction. + # + # The frequency-dependent beam compensation is therefore applied here, + # following CRIMAC ``calc_g(theta, phi, f)``: gain, beamwidths, and + # angle offsets are interpolated onto the TS(f) frequency grid, and the + # target angular distance from the beam axis is evaluated independently + # at each frequency. + # + # This differs from ``add_splitbeam_angle()``, which returns conventional + # echopype mechanical angles with the angle offsets already removed. + """ - Private method to compute Sv or TS from EK80 data, called by compute_Sv or compute_TS. + gain_db = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["gain"].sel(cal_channel_id=channel).values, + ) + + angle_offset_alongship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_alongship"].sel(cal_channel_id=channel).values, + ) + + angle_offset_athwartship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["angle_offset_athwartship"].sel(cal_channel_id=channel).values, + ) + + beamwidth_alongship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["beamwidth_alongship"].sel(cal_channel_id=channel).values, + ) + + beamwidth_athwartship = np.interp( + frequency, + self.vend["cal_frequency"].values, + self.vend["beamwidth_athwartship"].sel(cal_channel_id=channel).values, + ) + + beam_correction_db = self._get_beam_correction( + theta=theta, + phi=phi, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, + ) + + return 10 ** ((gain_db - beam_correction_db) / 10) + + def _cal_complex_samples_TS_spectrum( + self, + pc: xr.DataArray, + matched_filter: Dict, + point_locations: xr.Dataset, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str | None = None, + frequency_resolution: float | None = None, + ) -> xr.Dataset: + """Compute frequency-dependent target strength spectrum.""" + + if point_locations is None: + raise ValueError("TS_spectrum requires a point_locations dataset.") + + if not 0 <= split_front <= 1: + raise ValueError("split_front must be between 0 and 1.") + + pc_avg = _get_average_signal(pc) + + out_by_channel = [] + + for channel in self.beam["channel"].values: + if channel not in self.vend["cal_channel_id"].values: + continue + + points_ch = point_locations.where( + point_locations["channel"] == channel, + drop=True, + ) + + if points_ch.sizes.get("target_id", 0) == 0: + continue + + pc_avg_ch = pc_avg.sel(channel=channel) + range_ch = self.range_meter.sel(channel=channel) + + ts_list = [] + target_range_list = [] + theta_list = [] + phi_list = [] + ping_time_list = [] + target_id_list = [] + frequency_ref = None + + for point_idx in range(points_ch.sizes["target_id"]): + point_ping_time = points_ch["ping_time"].isel(target_id=point_idx).values + target_range = float(points_ch["target_range"].isel(target_id=point_idx)) + + target_id = points_ch["target_id"].isel(target_id=point_idx).values + + ping_idx = int( + np.argmin( + np.abs( + self.beam["ping_time"].values.astype("datetime64[ns]") + - np.datetime64(point_ping_time, "ns") + ) + ) + ) + + pc_avg_1d = pc_avg_ch.isel(ping_time=ping_idx) + range_1d = range_ch.isel(ping_time=ping_idx) + + valid = np.isfinite(pc_avg_1d) & np.isfinite(range_1d) + + if not bool(valid.any()): + continue + + pc_avg_1d = pc_avg_1d.where(valid, drop=True).values + range_1d = range_1d.where(valid, drop=True).values + + sound_speed = float( + self._select_param(self.env_params["sound_speed"], channel, ping_idx) + ) + transmit_power = float( + self._select_param(self.beam["transmit_power"], channel, ping_idx) + ) + sample_interval = float( + self._select_param(self.beam["sample_interval"], channel, ping_idx) + ) + f_start = float( + self._select_param(self.beam["transmit_frequency_start"], channel, ping_idx) + ) + f_stop = float( + self._select_param(self.beam["transmit_frequency_stop"], channel, ping_idx) + ) + + if frequency_resolution is not None: + n_f_points_local = int(np.floor((f_stop - f_start) / frequency_resolution)) + 1 + else: + n_f_points_local = n_f_points if n_f_points is not None else 1000 + + frequency = np.linspace(f_start, f_stop, n_f_points_local) + + if NFFT is None: + n_fft = int(2 ** np.ceil(np.log2(n_f_points_local))) + else: + n_fft = NFFT + + z_et = float( + self._select_param(self.cal_params["impedance_transducer"], channel, ping_idx) + ) + z_er = float( + self._select_param(self.cal_params["impedance_transceiver"], channel, None) + ) + + fs_dec = 1 / sample_interval + + absorption_f = self._compute_absorption_f( + frequency=frequency, + channel=channel, + ping_idx=ping_idx, + sound_speed=sound_speed, + ) + + if {"target_range_min", "target_range_max"}.issubset(points_ch.data_vars): + target_range_min = float( + points_ch["target_range_min"].isel(target_id=point_idx) + ) + target_range_max = float( + points_ch["target_range_max"].isel(target_id=point_idx) + ) + target_mask = (range_1d >= target_range_min) & (range_1d <= target_range_max) + else: # splitfront is added + idx_target = int(np.nanargmin(np.abs(range_1d - target_range))) + n_before = int(np.floor(split_front * n_fft)) + n_after = n_fft - n_before + idx_start = max(0, idx_target - n_before) + idx_stop = min(range_1d.size, idx_target + n_after) + target_mask = np.zeros(range_1d.size, dtype=bool) + target_mask[idx_start:idx_stop] = True + + if not np.any(target_mask): + continue + + pc_target = pc_avg_1d[target_mask] + + gamma_alongship = float( + self._select_param( + self.cal_params["angle_sensitivity_alongship"], + channel, + ping_idx, + ) + ) + gamma_athwartship = float( + self._select_param( + self.cal_params["angle_sensitivity_athwartship"], + channel, + ping_idx, + ) + ) + + theta_raw_da, phi_raw_da = _get_splitbeam_angles( + pc=pc.sel(channel=channel).isel(ping_time=ping_idx), + gamma_alongship=gamma_alongship, + gamma_athwartship=gamma_athwartship, + ) + + theta_raw = theta_raw_da.where(valid, drop=True).values + phi_raw = phi_raw_da.where(valid, drop=True).values + + idx_peak = int(np.nanargmax(np.abs(pc_avg_1d[target_mask]) ** 2)) + theta_t = float(theta_raw[target_mask][idx_peak]) + phi_t = float(phi_raw[target_mask][idx_peak]) + + if window is None: + win = np.ones(pc_target.size) + else: + win = get_window(window, pc_target.size) + + pc_target = pc_target * win + + mf_auto, _ = _get_autocorrelation( + matched_filter=matched_filter[channel], + n_window=n_fft, + ) + + mf_auto_red = _align_autocorrelation( + mf_auto=mf_auto, + pc_target=pc_target, + ) + + if mf_auto_red.size < n_fft: + mf_pad = np.zeros(n_fft, dtype=complex) + mf_pad[: mf_auto_red.size] = mf_auto_red + mf_auto_red = mf_pad + elif mf_auto_red.size > n_fft: + mf_auto_red = mf_auto_red[:n_fft] + + _, _, normalized_spectrum = _compute_ts_spectrum( + pc_target=pc_target, + mf_auto_red=mf_auto_red, + NFFT=n_fft, + frequency=frequency, + fs_dec=fs_dec, + ) + + power_spectrum = _compute_ts_spectrum_power( + normalized_spectrum=normalized_spectrum, + n_beams=self.beam["beam"].size, + z_et=z_et, + z_er=z_er, + ) + + gain_f = self._get_beam_compensated_gain( + channel=channel, + theta=theta_t, + phi=phi_t, + frequency=frequency, + ) + + ts = _compute_ts_spectrum_calibrated( + power_spectrum=power_spectrum, + target_range=target_range, + frequency=frequency, + sound_speed=sound_speed, + absorption_f=absorption_f, + transmit_power=transmit_power, + gain_f=gain_f, + ) + + if frequency_ref is None: + frequency_ref = frequency + elif frequency.shape != frequency_ref.shape: + continue + + ts_list.append(ts) + target_range_list.append(target_range) + theta_list.append(theta_t) + phi_list.append(phi_t) + ping_time_list.append(point_ping_time) + target_id_list.append(target_id) + + if not ts_list: + continue + + ds_ch = xr.Dataset( + { + "TS_spectrum": ( + ["channel", "target_id", "frequency"], + np.asarray(ts_list)[None, :, :], + ), + "target_range": ( + ["channel", "target_id"], + np.asarray(target_range_list)[None, :], + ), + "angle_alongship": ( + ["channel", "target_id"], + np.asarray(theta_list)[None, :], + ), + "angle_athwartship": ( + ["channel", "target_id"], + np.asarray(phi_list)[None, :], + ), + "ping_time": ( + ["channel", "target_id"], + np.asarray(ping_time_list)[None, :], + ), + }, + coords={ + "channel": [channel], + "target_id": np.asarray(target_id_list), + "frequency": frequency_ref, + }, + ) + + out_by_channel.append(ds_ch) + + if not out_by_channel: + raise ValueError("No valid TS_spectrum targets produced.") + + return xr.concat(out_by_channel, dim="channel") + + def _cal_complex_samples_f( + self, + cal_type: str, + frequency_resolution: float | None = None, + point_locations: xr.Dataset | None = None, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str | None = None, + ) -> xr.Dataset: + """Calibrate EK80 FM complex data to frequency-dependent TS(f).""" + + if self.waveform_mode not in ("FM", "BB") or self.encode_mode != "complex": + raise ValueError(f"{cal_type} is only supported for EK80 FM complex data.") + + tx_coeff = get_filter_coeff(self.vend) + fs = self.cal_params["receiver_sampling_frequency"] + + matched_filter, _ = get_transmit_signal( + self.beam, + tx_coeff, + self.waveform_mode, + fs, + self.drop_last_hanning_zero, + ) + + pc = _get_pulse_compressed_signal( + beam=self.beam, + matched_filter=matched_filter, + ) + + if cal_type == "TS_spectrum": + return self._cal_complex_samples_TS_spectrum( + pc=pc, + matched_filter=matched_filter, + point_locations=point_locations, + NFFT=NFFT, + n_f_points=n_f_points, + split_front=split_front, + window=window, + frequency_resolution=frequency_resolution, + ) + + raise ValueError(f"Unsupported calibration type: {cal_type}") + + def _compute_cal(self, cal_type, **kwargs) -> xr.Dataset: + """ + Private dispatcher for EK80 calibration. + + This routes calibration to one of three paths: + + 1. Power-sample CW calibration: + Used for EK60 and EK80 power-encoded CW data. + + 2. Complex-sample calibration: + Used for EK80 complex data in CW or BB/FM mode. For BB/FM data, + this path computes the conventional broadband-averaged Sv, Sp, or + band-averaged gridded TS product: pulse compression is applied, + transducer sectors are averaged, received power is computed, and the + standard calibration equation is used. + + 3. Frequency-dependent complex-sample calibration: + Used for EK80 FM complex data when computing TS(f). + This path keeps the pulse-compressed signal before final calibration + so that FFT-based spectral processing can be applied. Parameters ---------- cal_type : str - 'Sv' for calculating volume backscattering strength, or - 'TS' for calculating target strength + Calibration type. Supported values include ``"Sv"``, ``"Sp"``, + ``"TS"``, and ``"TS_spectrum"``. Returns ------- xr.Dataset - An xarray Dataset containing either Sv or TS. + Dataset containing the calibrated output requested by ``cal_type``. """ # Set flag_complex: True-complex cal, False-power cal flag_complex = ( True if self.waveform_mode == "BB" or self.encode_mode == "complex" else False ) - if flag_complex: - # Complex samples can be BB or CW + if cal_type in ("TS_spectrum",): + # Frequency-dependent broadband calibration: keep pulse-compressed + # signal for FFT-based Sv(f) or TS(f) processing. + ds_cal = self._cal_complex_samples_f(cal_type=cal_type, **kwargs) + elif flag_complex: + # Complex-sample calibration: BB/FM data are pulse-compressed and + # averaged over transducer sectors before computing + # Sv/Sp/band-averaged gridded TS. CW complex data are calibrated directly + # from complex samples. ds_cal = self._cal_complex_samples(cal_type=cal_type) else: - # Power samples only make sense for CW mode data + # Power-sample calibration: applies to power-encoded CW data. ds_cal = self._cal_power_samples(cal_type=cal_type) return ds_cal + ### Public API + def compute_Sv(self): """Compute volume backscattering strength (Sv). Returns ------- - Sv : xr.DataSet + Sv : xr.Dataset A DataSet containing volume backscattering strength (``Sv``) and the corresponding range (``echo_range``) in units meter. """ return self._compute_cal(cal_type="Sv") - def compute_TS(self): - """Compute target strength (TS). + def compute_Sp(self): + """ + Compute point scattering strength (Sp) from raw data. + + For CW data, Sp is computed on the range grid from power samples + or from complex samples converted to received power. For EK80 + broadband/FM complex data, Sp is computed after pulse compression + and represents a band-averaged point-scattering-strength echogram. Returns ------- - TS : xr.DataSet - A DataSet containing target strength (``TS``) - and the corresponding range (``echo_range``) in units meter. + Sp : xr.Dataset + Dataset containing point scattering strength (``Sp``) and the + corresponding range (``echo_range``) in metres. + """ + return self._compute_cal(cal_type="Sp") + + def compute_TS(self, **kwargs): + raise NotImplementedError( + "Gridded TS is deprecated. Use compute_Sp() first, detect targets, " + "then call echopype.calibrate.compute_TS(source_Sp, point_locations=...)." + ) + + def compute_TS_spectrum( + self, + point_locations: xr.Dataset, + NFFT: int | None = None, + n_f_points: int | None = None, + split_front: float = 0.25, + window: str = None, + frequency_resolution: float | None = None, + ): + """Compute frequency-dependent target strength spectrum TS(f). + + Returns + ------- + TS_f : xr.Dataset + A Dataset containing frequency-dependent target strength. """ - return self._compute_cal(cal_type="TS") + return self._compute_cal( + cal_type="TS_spectrum", + point_locations=point_locations, + NFFT=NFFT, + n_f_points=n_f_points, + split_front=split_front, + window=window, + frequency_resolution=frequency_resolution, + ) diff --git a/echopype/calibrate/ek80_complex.py b/echopype/calibrate/ek80_complex.py index 1717083b0..1ba65a848 100644 --- a/echopype/calibrate/ek80_complex.py +++ b/echopype/calibrate/ek80_complex.py @@ -9,6 +9,236 @@ from ..convert.set_groups_ek80 import DECIMATION, FILTER_IMAG, FILTER_REAL +def _get_transducer_halves( + pc: xr.DataArray, +) -> tuple[xr.DataArray, xr.DataArray, xr.DataArray, xr.DataArray]: + """Calculate half-transducer pulse-compressed signals. + + Equivalent to CRIMAC ``calcTransducerHalves`` for 4-sector transducers. + """ + if pc.sizes["beam"] != 4: + raise NotImplementedError( + "Transducer halves are only defined for 4-sector split-beam data." + ) + + pc_fore = 0.5 * (pc.isel(beam=2) + pc.isel(beam=3)) + pc_aft = 0.5 * (pc.isel(beam=0) + pc.isel(beam=1)) + pc_star = 0.5 * (pc.isel(beam=0) + pc.isel(beam=3)) + pc_port = 0.5 * (pc.isel(beam=1) + pc.isel(beam=2)) + + return pc_fore, pc_aft, pc_star, pc_port + + +def _get_splitbeam_angles( + pc: xr.DataArray, + gamma_alongship, + gamma_athwartship, +) -> tuple[xr.DataArray, xr.DataArray]: + """Calculate raw split-beam physical angles before angle-offset correction. + + For 4-sector data this follows CRIMAC ``calcAngles``. For 3-sector data, + the sector geometry follows the same convention used by ``add_splitbeam_angle``. + Angle offsets are not applied here because TS(f) beam compensation applies + frequency-dependent offsets later. + """ + if pc.sizes["beam"] == 4: + pc_fore, pc_aft, pc_star, pc_port = _get_transducer_halves(pc) + + y_theta = pc_fore * np.conj(pc_aft) + y_phi = pc_star * np.conj(pc_port) + + theta = np.rad2deg( + np.arcsin(np.arctan2(np.imag(y_theta), np.real(y_theta)) / gamma_alongship) + ) + + phi = np.rad2deg(np.arcsin(np.arctan2(np.imag(y_phi), np.real(y_phi)) / gamma_athwartship)) + else: + raise NotImplementedError( + f"Split-beam angle calculation is not implemented for {pc.sizes['beam']} sectors." + ) + + theta.name = "angle_alongship" + phi.name = "angle_athwartship" + + return theta, phi + + +def _compute_power_from_complex_signal( + signal: xr.DataArray, + z_et, + z_er, +) -> xr.DataArray: + """Calculate received electrical power from sector-level complex samples. + + The input is expected to retain the ``beam`` dimension. The function + averages over transducer sectors internally before converting the + complex signal to received electrical power. + + Equivalent to CRIMAC ``calcPower``. + """ + prx = ( + signal["beam"].size + * np.abs(signal.mean(dim="beam")) ** 2 + / (2 * np.sqrt(2)) ** 2 + * (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 + / np.abs(z_et) + ) + + prx = prx.where(prx > 0, 1e-20) + prx.name = "received_power" + + return prx + + +def _align_autocorrelation( + mf_auto: np.ndarray, + pc_target: np.ndarray, +) -> np.ndarray: + """Align matched-filter autocorrelation to target echo. + + Equivalent to CRIMAC ``alignAuto``. + """ + idx_peak_auto = np.argmax(np.abs(mf_auto)) + idx_peak_target = np.argmax(np.abs(pc_target)) + + left_samples = idx_peak_target + right_samples = len(pc_target) - idx_peak_target + + idx_start = max(0, idx_peak_auto - left_samples) + idx_stop = min(len(mf_auto), idx_peak_auto + right_samples) + + return mf_auto[idx_start:idx_stop] + + +def _compute_ts_spectrum( + pc_target: np.ndarray, + mf_auto_red: np.ndarray, + NFFT: int, + frequency: np.ndarray, + fs_dec: float, +): + """Compute target, autocorrelation, and normalised DFTs for TS spectrum. + + Equivalent to CRIMAC ``calcDFTforTS``, with explicit NFFT. + """ + frequency_index = np.mod( + np.floor(frequency / fs_dec * NFFT).astype(int), + NFFT, + ) + + pc_target_spectrum = np.fft.fft(pc_target, n=NFFT)[frequency_index] + mf_auto_red_spectrum = np.fft.fft(mf_auto_red, n=NFFT)[frequency_index] + + normalized_spectrum = pc_target_spectrum / mf_auto_red_spectrum + + return pc_target_spectrum, mf_auto_red_spectrum, normalized_spectrum + + +def _compute_ts_spectrum_calibrated( + power_spectrum: np.ndarray, + target_range: float, + frequency: np.ndarray, + sound_speed: float, + absorption_f: np.ndarray, + transmit_power: float, + gain_f: np.ndarray, +): + """Apply CRIMAC-style TS(f) calibration equation. + + Equivalent to CRIMAC ``calcTSf``. + """ + wavelength_f = sound_speed / frequency + + return ( + 10 * np.log10(power_spectrum) + + 40 * np.log10(target_range) + + 2 * absorption_f * target_range + - 10 * np.log10(transmit_power * wavelength_f**2 * gain_f**2 / (16 * np.pi**2)) + ) + + +def _compute_ts_spectrum_power( + normalized_spectrum: np.ndarray, + n_beams: int, + z_et: float, + z_er: float, +): + """Convert normalised TS(f) spectrum to received power spectrum. + + Equivalent to CRIMAC ``calcPowerFreqTS``. + """ + return _compute_complex_power( + normalized_spectrum=normalized_spectrum, + n_beams=n_beams, + z_et=z_et, + z_er=z_er, + ) + + +def _get_autocorrelation( + matched_filter: np.ndarray, + n_window: int, +): + """Get matched-filter autocorrelation spectrum. + + Equivalent to CRIMAC ``calcAutoCorrelation``. + """ + mf_auto = ( + np.convolve( + matched_filter, + np.conj(matched_filter[::-1]), + mode="full", + ) + / np.linalg.norm(matched_filter) ** 2 + ) + + mf_auto_spectrum = np.fft.fft(mf_auto, n=n_window) + + return mf_auto, mf_auto_spectrum + + +def _get_pulse_compressed_signal( + beam: xr.Dataset, + matched_filter: Dict, +) -> xr.DataArray: + """Calculate pulse-compressed complex samples for each transducer sector. + + Equivalent to CRIMAC ``calcPulseCompressedSignals``. + """ + pc = compress_pulse( + backscatter=beam["backscatter_r"] + 1j * beam["backscatter_i"], + chirp=matched_filter, + ) + pc = pc / get_norm_fac(chirp=matched_filter) + pc.name = "pulse_compressed_signal" + + return pc + + +def _get_average_signal( + signal: xr.DataArray, +) -> xr.DataArray: + """Average complex signal over transducer sectors. + + Equivalent to CRIMAC ``calcAverageSignal``. + """ + out = signal.mean(dim="beam") + out.name = "average_signal" + + return out + + +def _compute_complex_power( + normalized_spectrum: np.ndarray, + n_beams: int, + z_et: float, + z_er: float, +): + impedance_factor = (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 / np.abs(z_et) + + return n_beams * (np.abs(normalized_spectrum) / (2 * np.sqrt(2))) ** 2 * impedance_factor + + def tapered_chirp( fs, transmit_duration_nominal, diff --git a/echopype/mask/__init__.py b/echopype/mask/__init__.py index 7f0a27be8..b63fc26ba 100644 --- a/echopype/mask/__init__.py +++ b/echopype/mask/__init__.py @@ -1,3 +1,17 @@ -from .api import apply_mask, detect_seafloor, detect_shoal, frequency_differencing, regrid_mask +from .api import ( + apply_mask, + detect_seafloor, + detect_shoal, + detect_single_targets, + frequency_differencing, + regrid_mask, +) -__all__ = ["frequency_differencing", "apply_mask", "detect_seafloor", "detect_shoal", "regrid_mask"] +__all__ = [ + "frequency_differencing", + "apply_mask", + "detect_seafloor", + "detect_shoal", + "detect_single_targets", + "regrid_mask", +] diff --git a/echopype/mask/api.py b/echopype/mask/api.py index f5279f514..9322d5e53 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -11,9 +11,9 @@ import xarray as xr from flox.xarray import xarray_reduce -# for seafloor detection from echopype.mask.seafloor_detection.bottom_basic import bottom_basic from echopype.mask.seafloor_detection.bottom_blackwell import bottom_blackwell +from echopype.mask.single_target_detection.detect_from_Sp import detect_from_Sp from ..commongrid.utils import ( _convert_bins_to_interval_index, @@ -994,3 +994,76 @@ def detect_shoal( raise ValueError(f"Unsupported shoal detection method: {method}") return METHODS_SHOAL[method](ds, **params) + + +# Registry of supported methods for single_target_detection +METHODS_SINGLE_TARGET = { + "from_Sp": detect_from_Sp, +} + + +def detect_single_targets( + ds: xr.Dataset, + method: str, + params: dict, + waveform_mode: Literal["CW", "FM"] = "CW", +) -> xr.Dataset: + """ + Run single-target detection using the selected method. + + Parameters + ---------- + ds : xr.Dataset + Acoustic dataset containing the fields required by the selected method. + method : str + Name of the detection method to use. Currently supported for CW: + ``"from_Sp"`` + params : dict + Method-specific parameters. This argument is required and no defaults + are assumed. + waveform_mode : {"CW", "FM"}, default "CW" + Transmit waveform mode. FM single-target detection is not yet implemented. + + Returns + ------- + xr.Dataset + Per-target detection results with dimension ``target``. + """ + + if waveform_mode == "FM": + raise NotImplementedError("FM single-target detection is not yet implemented. ") + + if waveform_mode != "CW": + raise ValueError("waveform_mode must be 'CW' or 'FM'.") + + if method not in METHODS_SINGLE_TARGET: + raise ValueError(f"Unsupported single-target method: {method}") + + if params is None: + raise ValueError("No parameters given.") + + if "beam_type" not in ds: + raise ValueError("beam_type variable is missing from dataset.") + + beam_vals = np.unique(ds["beam_type"].values) + + if not np.all(np.isin(beam_vals, [1, 65])): + raise ValueError(f"Only split-beam data supported (beam_type 1 or 65). Found: {beam_vals}") + + out = METHODS_SINGLE_TARGET[method](ds, params) + + if not isinstance(out, xr.Dataset) or "single_target" not in out.dims: + raise TypeError(f"{method} must return an xr.Dataset with a 'single_target' dimension.") + + required = ("ping_time", "range_sample", "frequency_nominal") + missing = [v for v in required if v not in out] + if missing: + raise ValueError( + f"{method} output missing required field(s): {missing} (expected {list(required)})." + ) + + bad_dims = [v for v in required if out[v].dims != ("single_target",)] + if bad_dims: + raise ValueError(f"{method} field(s) must have dims ('single_target',): {bad_dims}.") + + return out diff --git a/echopype/mask/single_target_detection/detect_from_Sp.py b/echopype/mask/single_target_detection/detect_from_Sp.py new file mode 100644 index 000000000..5c0e9cf3e --- /dev/null +++ b/echopype/mask/single_target_detection/detect_from_Sp.py @@ -0,0 +1,619 @@ +from __future__ import annotations + +import numpy as np +import xarray as xr + +REQUIRED_PARAMS = { + "pldl_db", + "min_norm_pulse", + "max_norm_pulse", + "beam_comp_model", + "max_beam_comp_db", + "max_sd_minor_deg", + "max_sd_major_deg", +} + +OPTIONAL_PARAMS = { + "dec_tir_samples", + "bottom_offset_m", + "exclude_above_m", + "exclude_below_m", + "allow_nans_inside_envelope", +} + + +def _validate_params(params: dict) -> dict: + if params is None: + raise ValueError("No parameters given.") + + unknown = set(params.keys()) - (REQUIRED_PARAMS | OPTIONAL_PARAMS) + if unknown: + raise ValueError(f"Unknown parameters: {sorted(unknown)}") + + missing = REQUIRED_PARAMS - set(params.keys()) + if missing: + raise ValueError(f"Missing required parameters: {sorted(missing)}") + + if params["min_norm_pulse"] > params["max_norm_pulse"]: + raise ValueError("min_norm_pulse must be <= max_norm_pulse") + + return params + + +def _validate_from_Sp_dataset(ds_sp: xr.Dataset) -> xr.Dataset: + must = [ + "Sp", + "echo_range", + "angle_alongship", + "angle_athwartship", + "sound_absorption", + "sample_interval", + "tau_effective", + ] + + for v in must: + if v not in ds_sp: + raise ValueError(f"ds_sp missing required variable: {v}") + + sp_mat = ds_sp["Sp"] + + if sp_mat.dims != ("ping_time", "range_sample"): + raise ValueError("Expected ds_sp['Sp'] dims exactly ('ping_time', 'range_sample').") + + for v in ["echo_range", "angle_alongship", "angle_athwartship"]: + if ds_sp[v].dims != ("ping_time", "range_sample"): + raise ValueError(f"Expected ds_sp['{v}'] dims exactly ('ping_time', 'range_sample').") + + alpha = ds_sp["sound_absorption"] + + if alpha.ndim == 1: + ds_sp = ds_sp.assign(sound_absorption=alpha.broadcast_like(sp_mat)) + elif alpha.ndim == 0: + ds_sp = ds_sp.assign(sound_absorption=xr.zeros_like(sp_mat) + alpha) + elif alpha.ndim == 2: + pass + else: + raise ValueError("sound_absorption must be scalar, 1D(ping_time), or 2D like Sp.") + + return ds_sp + + +def _plike_from_sp( + sp_db: xr.DataArray, + r_m: xr.DataArray, + alpha_db_m: xr.DataArray, +) -> xr.DataArray: + r = xr.where(r_m > 0, r_m, np.nan) + return sp_db - 40.0 * xr.apply_ufunc(np.log10, r) - 2.0 * alpha_db_m * r_m + + +def _local_max_first_plateau(plike_mat: xr.DataArray) -> xr.DataArray: + prev = plike_mat.shift(range_sample=1) + nxt = plike_mat.shift(range_sample=-1) + + peak = (plike_mat > prev) & (plike_mat >= nxt) & ~(plike_mat == prev) + peak = peak & xr.apply_ufunc(np.isfinite, prev) + peak = peak & xr.apply_ufunc(np.isfinite, nxt) + peak = peak & xr.apply_ufunc(np.isfinite, plike_mat) + + return peak + + +def _nech_p_samples(ds_sp: xr.Dataset) -> np.ndarray: + tau = ds_sp["tau_effective"] + dt = ds_sp["sample_interval"] + + if tau.ndim == 0: + tau_vec = np.full(ds_sp.sizes["ping_time"], float(tau.values), dtype=float) + elif tau.ndim == 1 and tau.dims == ("ping_time",): + tau_vec = tau.values.astype(float) + else: + tau_vec = tau.isel(range_sample=0).values.astype(float) + + if dt.ndim == 0: + dt_vec = np.full(ds_sp.sizes["ping_time"], float(dt.values), dtype=float) + elif dt.ndim == 1 and dt.dims == ("ping_time",): + dt_vec = dt.values.astype(float) + else: + dt_vec = dt.isel(range_sample=0).values.astype(float) + + return tau_vec / dt_vec + + +def _envelope_bounds_1d( + plike_row: np.ndarray, + p: int, + thr: float, + allow_nans: bool, +) -> tuple[int | None, int | None]: + m = p + while m > 0: + v = plike_row[m - 1] + if not np.isfinite(v): + return (None, None) if not allow_nans else (m, p) + if v >= thr: + m -= 1 + else: + break + + last = p + while last < plike_row.size - 1: + v = plike_row[last + 1] + if not np.isfinite(v): + return (None, None) if not allow_nans else (m, last) + if v >= thr: + last += 1 + else: + break + + return m, last + + +def _beam_comp_db(ds_sp: xr.Dataset, params: dict) -> xr.DataArray: + model = params["beam_comp_model"] + + if model == "none": + return xr.zeros_like(ds_sp["Sp"]) + + if model == "provided": + if "beam_comp_db" not in ds_sp: + raise ValueError("beam_comp_db must exist if beam_comp_model='provided'") + return ds_sp["beam_comp_db"].broadcast_like(ds_sp["Sp"]) + + if model == "simrad_lobe": + th_al = ds_sp["angle_alongship"] + th_at = ds_sp["angle_athwartship"] + + bw_al = ds_sp["beamwidth_alongship"].broadcast_like(th_al) + bw_at = ds_sp["beamwidth_athwartship"].broadcast_like(th_at) + + # Angles from add_splitbeam_angle() are already offset-corrected. + x = 2.0 * th_al / bw_al + y = 2.0 * th_at / bw_at + + beam_comp_db = 6.0206 * (x**2 + y**2 - 0.18 * x**2 * y**2) + return beam_comp_db.broadcast_like(ds_sp["Sp"]) + + raise ValueError(f"Unknown beam_comp_model: {model}") + + +def _phase1_simple( + ds_sp: xr.Dataset, + params: dict, + beam_comp_db: xr.DataArray, +) -> xr.Dataset: + # Echoview Method 2 detects peaks on power-like data obtained by + # removing TVG/range and absorption from the TS operand. + # + # In echopype, Sp is already the uncompensated TS-like quantity, + # so reconstruct the power-like signal directly from Sp. + # + # Beam compensation is NOT part of the detection signal; it is only + # used as a peak-selection criterion and later for the final TS. + plike_mat = _plike_from_sp( + ds_sp["Sp"], + ds_sp["echo_range"], + ds_sp["sound_absorption"], + ) + + cand_mask = _local_max_first_plateau(plike_mat) + cand_mask = cand_mask & (beam_comp_db <= float(params["max_beam_comp_db"])) + + if params.get("dec_tir_samples") is not None: + dec_tir = int(params["dec_tir_samples"]) + idx = xr.DataArray( + np.arange(plike_mat.sizes["range_sample"]), + dims=("range_sample",), + coords={"range_sample": plike_mat["range_sample"]}, + ) + cand_mask = cand_mask & (idx >= dec_tir) + + if params.get("exclude_above_m") is not None: + cand_mask = cand_mask & (ds_sp["echo_range"] >= float(params["exclude_above_m"])) + + if params.get("exclude_below_m") is not None: + cand_mask = cand_mask & (ds_sp["echo_range"] <= float(params["exclude_below_m"])) + + if "bottom" in ds_sp and params.get("bottom_offset_m") is not None: + off = float(params["bottom_offset_m"]) + bottom2d = ds_sp["bottom"].broadcast_like(ds_sp["Sp"]) + cand_mask = cand_mask & (ds_sp["echo_range"] <= (bottom2d - off)) + + plike_np = plike_mat.values + cand_np = cand_mask.values + al_np = ds_sp["angle_alongship"].values + ath_np = ds_sp["angle_athwartship"].values + range_np = ds_sp["echo_range"].values + beam_comp_np = beam_comp_db.values + + nech_p = _nech_p_samples(ds_sp) + + pldl_db = float(params["pldl_db"]) + min_norm_pulse = float(params["min_norm_pulse"]) + max_norm_pulse = float(params["max_norm_pulse"]) + max_sd_minor_deg = float(params["max_sd_minor_deg"]) + max_sd_major_deg = float(params["max_sd_major_deg"]) + allow_nans = bool(params.get("allow_nans_inside_envelope", False)) + + ping_index_list = [] + range_sample_list = [] + iinf_list = [] + isup_list = [] + pulse_len_samples_list = [] + norm_pulse_len_list = [] + plike_peak_list = [] + target_range_list = [] + beam_comp_db_list = [] + angle_minor_sd_deg_list = [] + angle_major_sd_deg_list = [] + + for it in range(plike_np.shape[0]): + peaks = np.where(cand_np[it])[0] + if peaks.size == 0: + continue + + nch = nech_p[it] + if not np.isfinite(nch) or nch <= 0: + continue + + plike_row = plike_np[it] + ali = al_np[it] + athi = ath_np[it] + + for p in peaks: + plike_peak = plike_row[p] + if not np.isfinite(plike_peak): + continue + + iinf, isup = _envelope_bounds_1d( + plike_row, + int(p), + plike_peak - pldl_db, + allow_nans=allow_nans, + ) + if iinf is None: + continue + + pulse_len_samples = isup - iinf + 1 + norm_pulse_len = pulse_len_samples / nch + + if norm_pulse_len < min_norm_pulse or norm_pulse_len > max_norm_pulse: + continue + + seg_al = ali[iinf : isup + 1] * 180.0 / np.pi + seg_ath = athi[iinf : isup + 1] * 180.0 / np.pi + + angle_minor_sd_deg = float(np.nanstd(seg_ath)) + angle_major_sd_deg = float(np.nanstd(seg_al)) + + if angle_minor_sd_deg > max_sd_minor_deg: + continue + + if angle_major_sd_deg > max_sd_major_deg: + continue + + ping_index_list.append(it) + range_sample_list.append(int(p)) + iinf_list.append(int(iinf)) + isup_list.append(int(isup)) + pulse_len_samples_list.append(int(pulse_len_samples)) + norm_pulse_len_list.append(float(norm_pulse_len)) + plike_peak_list.append(float(plike_peak)) + target_range_list.append(float(range_np[it, p])) + beam_comp_db_list.append(float(beam_comp_np[it, p])) + angle_minor_sd_deg_list.append(angle_minor_sd_deg) + angle_major_sd_deg_list.append(angle_major_sd_deg) + + return xr.Dataset( + data_vars=dict( + ping_index=( + "single_target", + np.array(ping_index_list, dtype=np.int64), + ), + range_sample=( + "single_target", + np.array(range_sample_list, dtype=np.int64), + ), + iinf=( + "single_target", + np.array(iinf_list, dtype=np.int64), + ), + isup=( + "single_target", + np.array(isup_list, dtype=np.int64), + ), + pulse_len_samples=( + "single_target", + np.array(pulse_len_samples_list, dtype=np.int64), + ), + norm_pulse_len=( + "single_target", + np.array(norm_pulse_len_list, dtype=np.float64), + ), + plike_peak=( + "single_target", + np.array(plike_peak_list, dtype=np.float64), + ), + single_target_range=( + "single_target", + np.array(target_range_list, dtype=np.float64), + ), + beam_comp_db=( + "single_target", + np.array(beam_comp_db_list, dtype=np.float64), + ), + single_target_athwartship_angle_sd=( + "single_target", + np.array(angle_minor_sd_deg_list, dtype=np.float64), + ), + single_target_alongship_angle_sd=( + "single_target", + np.array(angle_major_sd_deg_list, dtype=np.float64), + ), + ), + coords={ + "single_target": np.arange( + len(range_sample_list), + dtype=np.int64, + ) + }, + ) + + +def _reject_overlaps_per_ping(feats: xr.Dataset) -> xr.Dataset: + if feats.sizes.get("single_target", 0) <= 1: + return feats + + ping_idx = feats["ping_index"].values + range_sample = feats["range_sample"].values + iinf = feats["iinf"].values + isup = feats["isup"].values + plike_peak = feats["plike_peak"].values + + keep = np.ones(feats.sizes["single_target"], dtype=bool) + + for it in np.unique(ping_idx): + ii = np.where(ping_idx == it)[0] + if ii.size <= 1: + continue + + # Echoview: screen pulses from low to high range/depth. + # Use peak sample order, not envelope-start order. + order = ii[np.argsort(range_sample[ii])] + accepted = [] + + for j in order: + if not accepted: + accepted.append(j) + continue + + k = accepted[-1] + + if iinf[j] > isup[k]: + accepted.append(j) + continue + + # If pulses overlap, reject the lower-power / lower-TS one. + if plike_peak[j] >= plike_peak[k]: + keep[k] = False + accepted[-1] = j + else: + keep[j] = False + + return feats.isel(single_target=keep) + + +def _pack_targets(feats: xr.Dataset, ds_sp: xr.Dataset) -> xr.Dataset: + n_targets = feats.sizes.get("single_target", 0) + channel_value = ds_sp["channel"].item() + + if n_targets == 0: + return xr.Dataset( + data_vars=dict( + channel=( + "single_target", + np.array([], dtype=object), + ), + ping_time=( + "single_target", + np.array([], dtype=ds_sp["ping_time"].dtype), + ), + range_sample=( + "single_target", + np.array([], dtype=np.int64), + ), + frequency_nominal=( + "single_target", + np.array([], dtype=np.float64), + ), + ping_index=( + "single_target", + np.array([], dtype=np.int64), + ), + iinf=( + "single_target", + np.array([], dtype=np.int64), + ), + isup=( + "single_target", + np.array([], dtype=np.int64), + ), + pulse_len_samples=( + "single_target", + np.array([], dtype=np.int64), + ), + norm_pulse_len=( + "single_target", + np.array([], dtype=np.float64), + ), + single_target_range=( + "single_target", + np.array([], dtype=np.float64), + ), + single_target_alongship_angle=( + "single_target", + np.array([], dtype=np.float64), + ), + single_target_athwartship_angle=( + "single_target", + np.array([], dtype=np.float64), + ), + single_target_athwartship_angle_sd=( + "single_target", + np.array([], dtype=np.float64), + ), + single_target_alongship_angle_sd=( + "single_target", + np.array([], dtype=np.float64), + ), + beam_comp_db=( + "single_target", + np.array([], dtype=np.float64), + ), + plike_peak=( + "single_target", + np.array([], dtype=np.float64), + ), + ), + coords={ + "single_target": np.arange( + n_targets, + dtype=np.int64, + ) + }, + attrs=dict(method="from_Sp"), + ) + + it = feats["ping_index"].values.astype(np.int64) + p = feats["range_sample"].values.astype(np.int64) + + ping_time = ds_sp["ping_time"].values[it] + single_target_alongship_angle = ds_sp["angle_alongship"].values[it, p] * 180.0 / np.pi + single_target_athwartship_angle = ds_sp["angle_athwartship"].values[it, p] * 180.0 / np.pi + + fn = ds_sp["frequency_nominal"] + if fn.ndim == 0: + freq_val = float(fn.values) + else: + freq_val = float(fn.values[0]) + + frequency_nominal = np.full(n_targets, freq_val, dtype=np.float64) + channel = np.full( + n_targets, + channel_value, + dtype=object, + ) + + return xr.Dataset( + data_vars=dict( + channel=( + "single_target", + channel, + ), + ping_time=( + "single_target", + ping_time, + ), + range_sample=( + "single_target", + p, + ), + frequency_nominal=( + "single_target", + frequency_nominal, + ), + ping_index=( + "single_target", + it, + ), + iinf=( + "single_target", + feats["iinf"].values.astype(np.int64), + ), + isup=( + "single_target", + feats["isup"].values.astype(np.int64), + ), + pulse_len_samples=( + "single_target", + feats["pulse_len_samples"].values.astype(np.int64), + ), + norm_pulse_len=( + "single_target", + feats["norm_pulse_len"].values.astype(np.float64), + ), + single_target_range=( + "single_target", + feats["single_target_range"].values.astype(np.float64), + ), + single_target_alongship_angle=( + "single_target", + single_target_alongship_angle.astype(np.float64), + ), + single_target_athwartship_angle=( + "single_target", + single_target_athwartship_angle.astype(np.float64), + ), + single_target_athwartship_angle_sd=( + "single_target", + feats["single_target_athwartship_angle_sd"].values.astype(np.float64), + ), + single_target_alongship_angle_sd=( + "single_target", + feats["single_target_alongship_angle_sd"].values.astype(np.float64), + ), + beam_comp_db=( + "single_target", + feats["beam_comp_db"].values.astype(np.float64), + ), + plike_peak=( + "single_target", + feats["plike_peak"].values.astype(np.float64), + ), + ), + coords={ + "single_target": np.arange( + n_targets, + dtype=np.int64, + ) + }, + attrs=dict(method="from_Sp"), + ) + + +def detect_from_Sp(ds_sp: xr.Dataset, params: dict) -> xr.Dataset: + """ + Detect single-target candidate locations from point scattering strength Sp. + + This follows the detection part of Echoview split-beam Method 2, but stops + before target-strength calculation. TS and TS(f) should be computed later + from the returned target locations. + """ + params = _validate_params(params) + ds_sp = _validate_from_Sp_dataset(ds_sp) + + # no need for copy? + ds_sp = ds_sp.copy() + + deg2rad = np.pi / 180.0 + + ds_sp["angle_alongship"] = ds_sp["angle_alongship"] * deg2rad + ds_sp["angle_athwartship"] = ds_sp["angle_athwartship"] * deg2rad + + if params["beam_comp_model"] == "simrad_lobe": + for v in [ + "beamwidth_alongship", + "beamwidth_athwartship", + "angle_offset_alongship", + "angle_offset_athwartship", + ]: + if v not in ds_sp: + raise ValueError(f"ds_sp missing required variable for beam compensation: {v}") + ds_sp[v] = ds_sp[v] * deg2rad + + beam_comp_db = _beam_comp_db(ds_sp, params) + + feats = _phase1_simple(ds_sp, params, beam_comp_db) + feats = _reject_overlaps_per_ping(feats) + + return _pack_targets(feats, ds_sp) diff --git a/echopype/tests/calibrate/test_calibrate.py b/echopype/tests/calibrate/test_calibrate.py index c9628b67c..2ae7adaa1 100644 --- a/echopype/tests/calibrate/test_calibrate.py +++ b/echopype/tests/calibrate/test_calibrate.py @@ -101,7 +101,7 @@ def test_compute_Sv_ek60_matlab(ek60_path): # Calibrate to get Sv ds_Sv = ep.calibrate.compute_Sv(echodata) - ds_TS = ep.calibrate.compute_TS(echodata) + ds_Sp = ep.calibrate.compute_Sp(echodata) # Load matlab outputs and test @@ -117,10 +117,10 @@ def check_output(da_cmp, cal_type): assert np.allclose(pyel_vals, ep_vals) # Check Sv - check_output(ds_Sv['Sv'], 'Sv') + check_output(ds_Sv["Sv"], "Sv") - # Check TS - check_output(ds_TS['TS'], 'Sp') + # Check Sp + check_output(ds_Sp["Sp"], "Sp") @pytest.mark.integration @@ -137,10 +137,10 @@ def test_compute_Sv_ek60_duplicated_freq(ek60_path): # Calibrate to get Sv ds_Sv = ep.calibrate.compute_Sv(echodata) - ds_TS = ep.calibrate.compute_TS(echodata) + ds_Sp = ep.calibrate.compute_Sp(echodata) assert isinstance(ds_Sv, xr.Dataset) - assert isinstance(ds_TS, xr.Dataset) + assert isinstance(ds_Sp, xr.Dataset) @pytest.mark.integration @@ -446,15 +446,27 @@ def test_check_echodata_backscatter_size( @pytest.mark.integration def test_fm_equals_bb(ek80_path): - """Check that waveform_mode='BB' and waveform_mode='FM' result in the same Sv/TS.""" - # Open Raw and Compute both Sv and both TS - ed = ep.open_raw(ek80_path / "D20170912-T234910.raw", sonar_model = "EK80") - ds_Sv_bb = ep.calibrate.compute_Sv(ed, waveform_mode="BB", encode_mode="complex") - ds_Sv_fm = ep.calibrate.compute_Sv(ed, waveform_mode="FM", encode_mode="complex") - ds_TS_bb = ep.calibrate.compute_TS(ed, waveform_mode="BB", encode_mode="complex") - ds_TS_fm = ep.calibrate.compute_TS(ed, waveform_mode="FM", encode_mode="complex") - - # Check that they are equal + """Check that waveform_mode='BB' and waveform_mode='FM' produce identical Sv and TS.""" + ed = ep.open_raw(ek80_path / "D20170912-T234910.raw", sonar_model="EK80") + + with pytest.deprecated_call(match="'BB' is deprecated"): + ds_Sv_bb = ep.calibrate.compute_Sv( + ed, waveform_mode="BB", encode_mode="complex" + ) + + ds_Sv_fm = ep.calibrate.compute_Sv( + ed, waveform_mode="FM", encode_mode="complex" + ) + + with pytest.deprecated_call(match="'BB' is deprecated"): + ds_TS_bb = ep.calibrate.compute_TS( + ed, waveform_mode="BB", encode_mode="complex" + ) + + ds_TS_fm = ep.calibrate.compute_TS( + ed, waveform_mode="FM", encode_mode="complex" + ) + assert ds_Sv_bb.equals(ds_Sv_fm) assert ds_TS_bb.equals(ds_TS_fm) diff --git a/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py new file mode 100644 index 000000000..507db8cce --- /dev/null +++ b/echopype/tests/calibrate/test_calibrate_ek80_broadband_crimac.py @@ -0,0 +1,219 @@ +import numpy as np +import pytest +import xarray as xr + +import echopype as ep + +pytestmark = pytest.mark.integration + +CRIMAC_CHANNEL = "WBT 747022-15 ES120-7CD_ES" +CRIMAC_PING_INDEX = 509 + + +@pytest.fixture(scope="module") +def ts_spectrum_example_path(test_path): + return test_path["TS_SPECTRUM_EXAMPLE"] + + +@pytest.fixture(scope="module") +def ts_raw_path(ts_spectrum_example_path): + return ts_spectrum_example_path / "IMR-D20211215-T143432-TSf.raw" + + +@pytest.fixture(scope="module") +def ts_ref(ts_spectrum_example_path): + return np.load( + ts_spectrum_example_path / "crimac_tsf_reference_outputs.npz", + allow_pickle=True, + ) + + +@pytest.fixture(scope="module") +def ts_echodata(ts_raw_path): + return ep.open_raw(ts_raw_path, sonar_model="EK80") + + +def _target_locations_from_crimac(ed, ref, channel, ping_index): + ping_time = ed["Sonar/Beam_group1"]["ping_time"].isel(ping_time=ping_index).values + + return xr.Dataset( + data_vars={ + "target_range": ("target_id", [float(ref["r_t"])]), + "angle_alongship": ("target_id", [float(ref["theta_t"])]), + "angle_athwartship": ("target_id", [float(ref["phi_t"])]), + "target_range_min": ("target_id", [float(ref["dum_r"][0])]), + "target_range_max": ("target_id", [float(ref["dum_r"][-1])]), + }, + coords={ + "target_id": [0], + "ping_time": ("target_id", [ping_time]), + "channel": ("target_id", [channel]), + }, + ) + + +def test_compute_sp_fm_complex_runs(ts_echodata): + ds = ep.calibrate.compute_Sp( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + ) + + assert "Sp" in ds + assert set(("channel", "ping_time", "range_sample")).issubset(ds["Sp"].dims) + assert np.isfinite(ds["Sp"]).any() + + +def test_frequency_dependent_absorption_matches_crimac(ts_echodata, ts_ref): + cal_obj = ep.calibrate.calibrate_ek.CalibrateEK80( + echodata=ts_echodata, + waveform_mode="BB", #TODO change to FM after deprecation of BB + encode_mode="complex", + env_params=None, + cal_params=None, + ) + + sound_speed = float(cal_obj.env_params["sound_speed"]) + + absorption_f = cal_obj._compute_absorption_f( + frequency=ts_ref["f_m"], + channel=CRIMAC_CHANNEL, + ping_idx=CRIMAC_PING_INDEX, + sound_speed=sound_speed, + ) + + np.testing.assert_allclose( + absorption_f, + ts_ref["alpha_m"], + atol=1e-8, + rtol=0.0, + ) + + +def test_compute_ts_spectrum_matches_crimac(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + ds = ep.calibrate.compute_TS_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ts = ( + ds["TS_spectrum"] + .sel(channel=CRIMAC_CHANNEL) + .isel(target_id=0) + .values + ) + + assert ts.shape == ts_ref["TS_m"].shape + np.testing.assert_allclose(ts, ts_ref["TS_m"], atol=0.5, rtol=0.0) + + +def test_compute_Sv_spectrum_not_implemented(ts_echodata): + """ + Test that compute_Sv_spectrum raises NotImplementedError + for EK80 broadband complex data, since this is not currently implemented. + """ + with pytest.raises(NotImplementedError): + ep.calibrate.compute_Sv_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + ) + +@pytest.mark.parametrize("window", [None, "boxcar", "hann", "hamming", ("tukey", 0.25)]) +def test_compute_ts_spectrum_accepts_scipy_windows(ts_echodata, ts_ref, window): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + ds = ep.calibrate.compute_TS_spectrum( + ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + window=window, + ) + + assert "TS_spectrum" in ds + assert np.isfinite(ds["TS_spectrum"]).any() + +def test_compute_ts_spectrum_none_window_matches_boxcar(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_none = ep.calibrate.compute_TS_spectrum(**kwargs, window=None) + ds_boxcar = ep.calibrate.compute_TS_spectrum(**kwargs, window="boxcar") + + xr.testing.assert_allclose(ds_none["TS_spectrum"], ds_boxcar["TS_spectrum"]) + + +def test_compute_ts_spectrum_explicit_range_ignores_split_front(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_025 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.25) + ds_075 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.75) + + xr.testing.assert_allclose(ds_025["TS_spectrum"], ds_075["TS_spectrum"]) + +def test_compute_ts_spectrum_target_range_only_uses_split_front(ts_echodata, ts_ref): + point_locations = _target_locations_from_crimac( + ed=ts_echodata, + ref=ts_ref, + channel=CRIMAC_CHANNEL, + ping_index=CRIMAC_PING_INDEX, + ).drop_vars(["target_range_min", "target_range_max"]) + + kwargs = dict( + echodata=ts_echodata, + waveform_mode="FM", + encode_mode="complex", + point_locations=point_locations, + n_f_points=ts_ref["f_m"].size, + ) + + ds_025 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.25) + ds_075 = ep.calibrate.compute_TS_spectrum(**kwargs, split_front=0.75) + + assert not np.allclose( + ds_025["TS_spectrum"].values, + ds_075["TS_spectrum"].values, + equal_nan=True, + ) \ No newline at end of file diff --git a/echopype/tests/calibrate/test_ek80_complex.py b/echopype/tests/calibrate/test_ek80_complex.py index 595800496..e5bf54853 100644 --- a/echopype/tests/calibrate/test_ek80_complex.py +++ b/echopype/tests/calibrate/test_ek80_complex.py @@ -2,7 +2,15 @@ import numpy as np import xarray as xr -from echopype.calibrate.ek80_complex import get_vend_filter_EK80 +from echopype.calibrate.ek80_complex import ( + get_vend_filter_EK80, + _get_average_signal, + _compute_power_from_complex_signal, + _compute_ts_spectrum_power, + _align_autocorrelation, + _compute_ts_spectrum, + _compute_ts_spectrum_calibrated, +) pytestmark = pytest.mark.unit @@ -73,3 +81,103 @@ def test_get_vend_filter_EK80(ch_num, filter_len, has_nan): assert sel_vend[var_df].values == get_vend_filter_EK80( vend, channel_id=ch, filter_name=filter_name, param_type="decimation" ) + +def test_get_average_signal(): + signal = xr.DataArray( + np.array([[1 + 1j, 3 + 3j], [5 + 5j, 7 + 7j]]), + dims=("range_sample", "beam"), + coords={"beam": [0, 1]}, + ) + + out = _get_average_signal(signal) + expected = signal.mean(dim="beam") + + xr.testing.assert_allclose(out, expected) + assert out.name == "average_signal" + + +def test_compute_power_from_complex_signal(): + signal = xr.DataArray( + np.array([[1 + 1j, 3 + 3j]]), + dims=("range_sample", "beam"), + coords={"beam": [0, 1]}, + ) + + z_et = 75.0 + z_er = 5400.0 + n_beams = signal["beam"].size + avg = signal.mean(dim="beam") + + expected = ( + n_beams + * np.abs(avg) ** 2 + / (2 * np.sqrt(2)) ** 2 + * (np.abs(z_er + z_et) / np.abs(z_er)) ** 2 + / np.abs(z_et) + ) + + out = _compute_power_from_complex_signal(signal, z_et=z_et, z_er=z_er) + + xr.testing.assert_allclose(out, expected) + assert out.name == "received_power" + + +def test_align_autocorrelation(): + mf_auto = np.array([0, 0, 1, 0.5, 0.25, 0.1]) + pc_target = np.array([0.2, 1.0, 0.3]) + + out = _align_autocorrelation(mf_auto=mf_auto, pc_target=pc_target) + + np.testing.assert_allclose(out, np.array([0, 1, 0.5])) + + +def test_compute_ts_spectrum(): + pc_target = np.array([1.0, 2.0, 1.0]) + mf_auto_red = np.array([1.0, 1.0, 1.0]) + frequency = np.array([0.0, 1.0, 2.0]) + fs_dec = 8.0 + + y_pc, y_mf, y_norm = _compute_ts_spectrum( + pc_target=pc_target, + mf_auto_red=mf_auto_red, + NFFT=8, + frequency=frequency, + fs_dec=fs_dec, + ) + + assert y_pc.shape == frequency.shape + assert y_mf.shape == frequency.shape + assert y_norm.shape == frequency.shape + np.testing.assert_allclose(y_norm, y_pc / y_mf) + + +def test_compute_ts_spectrum_calibrated(): + power_spectrum = np.array([1e-12, 2e-12]) + frequency = np.array([90000.0, 100000.0]) + sound_speed = 1500.0 + target_range = 10.0 + absorption_f = np.array([0.03, 0.04]) + transmit_power = 1000.0 + gain_f = np.array([100.0, 120.0]) + + out = _compute_ts_spectrum_calibrated( + power_spectrum=power_spectrum, + target_range=target_range, + frequency=frequency, + sound_speed=sound_speed, + absorption_f=absorption_f, + transmit_power=transmit_power, + gain_f=gain_f, + ) + + wavelength = sound_speed / frequency + expected = ( + 10 * np.log10(power_spectrum) + + 40 * np.log10(target_range) + + 2 * absorption_f * target_range + - 10 * np.log10( + transmit_power * wavelength**2 * gain_f**2 / (16 * np.pi**2) + ) + ) + + np.testing.assert_allclose(out, expected) \ No newline at end of file diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 734d9836c..f1219b5e8 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -98,7 +98,9 @@ def _unpack(fname, action, pooch_instance): time.sleep(1) with ZipFile(z, "r") as f: - f.extractall(out) + for member in f.infolist(): + member.filename = member.filename.replace("\\", "/") + f.extract(member, out) # flatten single nested dir if needed try: diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index 631a4fe75..3a5d0f8d0 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -4,6 +4,7 @@ import os import pandas as pd # noqa: F401 +from echopype.calibrate import compute_TS import xarray as xr import numpy as np import dask.array @@ -22,6 +23,10 @@ # for schoals from echopype.mask import detect_shoal + +# for single targets +from echopype.mask import detect_single_targets + from scipy import ndimage as ndi from typing import List, Union, Optional # noqa: F811 @@ -2037,3 +2042,815 @@ def test_echoview_mincan_no_linking(): # Label the mask and confirm they are separate components (not connected) _, nlab = ndi.label(mask.values, structure=np.ones((3, 3), dtype=bool)) assert nlab == 2 + +# test for single target detections + +# Helpers: base coords +def _coords_ping_range(n_ping=5, n_range=10): + return { + "ping_time": np.arange(n_ping), + "range_sample": np.arange(n_range), + } + +# Stub: from-Sv detector input (3D with channel) +def _make_ds_from_Sv_minimal( + n_ping=5, + n_range=10, + channels=("chan1",), +) -> xr.Dataset: + coords = _coords_ping_range(n_ping, n_range) + ping_time = coords["ping_time"] + range_sample = coords["range_sample"] + channel = np.array(list(channels)) + + # Core 3D fields + Sv = xr.DataArray( + np.full((len(channel), n_ping, n_range), -90.0, dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="Sv", + ) + + al = xr.DataArray( + np.zeros((len(channel), n_ping, n_range), dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="angle_alongship", + ) + + ath = xr.DataArray( + np.zeros((len(channel), n_ping, n_range), dtype=float), + dims=("channel", "ping_time", "range_sample"), + coords={"channel": channel, "ping_time": ping_time, "range_sample": range_sample}, + name="angle_athwartship", + ) + + # Range axis (2D accepted) + echo_range = xr.DataArray( + np.tile(np.linspace(1.0, float(n_range), n_range)[None, :], (n_ping, 1)), + dims=("ping_time", "range_sample"), + coords={"ping_time": ping_time, "range_sample": range_sample}, + name="echo_range", + ) + + # SHIFT alignment requirement + start_depth_m = xr.DataArray( + np.zeros((n_ping, len(channel)), dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="start_depth_m", + ) + + # scalar-ish + sound_speed = xr.DataArray(1500.0, name="sound_speed") + equivalent_beam_angle = xr.DataArray(-20.0, name="equivalent_beam_angle") + sa_correction = xr.DataArray(0.0, name="sa_correction") + sample_interval = xr.DataArray(4e-5, name="sample_interval") # 40 us + + # tau_effective: channel-only + tau_effective = xr.DataArray( + np.full((len(channel),), 8e-4, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="tau_effective", + ) + + transmit_duration_nominal = xr.DataArray( + np.full((n_ping, len(channel)), 1e-3, dtype=float), # 1 ms + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="transmit_duration_nominal", + ) + + # ping_time x channel required + sound_absorption = xr.DataArray( + np.full((n_ping, len(channel)), 0.003, dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="sound_absorption", + ) + transducer_depth = xr.DataArray( + np.full((n_ping, len(channel)), 2.0, dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="transducer_depth", + ) + heave_compensation = xr.DataArray( + np.zeros((n_ping, len(channel)), dtype=float), + dims=("ping_time", "channel"), + coords={"ping_time": ping_time, "channel": channel}, + name="heave_compensation", + ) + + # channel-only beam geometry + beamwidth_alongship = xr.DataArray( + np.full((len(channel),), 7.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="beamwidth_alongship", + ) + beamwidth_athwartship = xr.DataArray( + np.full((len(channel),), 7.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="beamwidth_athwartship", + ) + angle_offset_alongship = xr.DataArray( + np.zeros((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_offset_alongship", + ) + angle_offset_athwartship = xr.DataArray( + np.zeros((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_offset_athwartship", + ) + angle_sensitivity_alongship = xr.DataArray( + np.ones((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_sensitivity_alongship", + ) + angle_sensitivity_athwartship = xr.DataArray( + np.ones((len(channel),), dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="angle_sensitivity_athwartship", + ) + + beam_type = xr.DataArray( + np.full((len(channel),), 1, dtype=np.int16), + dims=("channel",), + coords={"channel": channel}, + name="beam_type", + ) + + frequency_nominal = xr.DataArray( + np.full((len(channel),), 38000.0, dtype=float), + dims=("channel",), + coords={"channel": channel}, + name="frequency_nominal", + ) + + ds = xr.Dataset( + data_vars=dict( + Sv=Sv, + angle_alongship=al, + angle_athwartship=ath, + echo_range=echo_range, + start_depth_m=start_depth_m, + sound_speed=sound_speed, + transmit_duration_nominal=transmit_duration_nominal, + tau_effective=tau_effective, + sample_interval=sample_interval, + sound_absorption=sound_absorption, + equivalent_beam_angle=equivalent_beam_angle, + sa_correction=sa_correction, + beamwidth_alongship=beamwidth_alongship, + beamwidth_athwartship=beamwidth_athwartship, + angle_offset_alongship=angle_offset_alongship, + angle_offset_athwartship=angle_offset_athwartship, + angle_sensitivity_alongship=angle_sensitivity_alongship, + angle_sensitivity_athwartship=angle_sensitivity_athwartship, + transducer_depth=transducer_depth, + heave_compensation=heave_compensation, + beam_type=beam_type, + frequency_nominal=frequency_nominal, + ), + coords=dict(channel=channel, ping_time=ping_time, range_sample=range_sample), + ) + return ds + + +def _make_ds_from_Sv_wrong_dims_missing_range(channels=("chan1",), n_ping=5) -> xr.Dataset: + # Sv missing range_sample dim intentionally (and also missing range_sample coord/dim) + da = xr.DataArray( + np.full((len(channels), n_ping), -90.0, dtype=float), + dims=("channel", "ping_time"), + coords={"channel": list(channels), "ping_time": np.arange(n_ping)}, + name="Sv", + ) + ds = da.to_dataset() + ds = ds.assign( + beam_type=xr.DataArray( + np.ones((len(channels),), dtype=np.int16), + dims=("channel",), + coords={"channel": list(channels)}, + ) + ) + return ds + + +# Stub: from-Sp detector input (2D, no channel dim) +def _make_ds_from_Sp_minimal(n_ping=5, n_range=10) -> xr.Dataset: + coords = _coords_ping_range(n_ping, n_range) + + TS = xr.DataArray( + np.full((n_ping, n_range), -90.0, dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="TS", + ) + echo_range = xr.DataArray( + np.tile(np.linspace(1.0, float(n_range), n_range)[None, :], (n_ping, 1)), + dims=("ping_time", "range_sample"), + coords=coords, + name="echo_range", + ) + angle_al = xr.DataArray( + np.zeros((n_ping, n_range), dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_alongship", + ) + angle_ath = xr.DataArray( + np.zeros((n_ping, n_range), dtype=float), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_athwartship", + ) + + sound_absorption = xr.DataArray(0.003, name="sound_absorption") + sample_interval = xr.DataArray(4e-5, name="sample_interval") + tau_effective = xr.DataArray(8e-4, name="tau_effective") + + frequency_nominal = xr.DataArray(38000.0, name="frequency_nominal") + + # beam compensation inputs for simrad_lobe + beamwidth_al = xr.DataArray(7.0, name="beamwidth_alongship") + beamwidth_at = xr.DataArray(7.0, name="beamwidth_athwartship") + angle_off_al = xr.DataArray(0.0, name="angle_offset_alongship") + angle_off_at = xr.DataArray(0.0, name="angle_offset_athwartship") + + beam_type = xr.DataArray(np.int16(1), name="beam_type") + + ds = xr.Dataset( + data_vars=dict( + TS=TS, + echo_range=echo_range, + angle_alongship=angle_al, + angle_athwartship=angle_ath, + sound_absorption=sound_absorption, + sample_interval=sample_interval, + tau_effective=tau_effective, + frequency_nominal=frequency_nominal, + beamwidth_alongship=beamwidth_al, + beamwidth_athwartship=beamwidth_at, + angle_offset_alongship=angle_off_al, + angle_offset_athwartship=angle_off_at, + beam_type=beam_type, + ), + coords=coords, + ) + return ds + +# --------------------------------------------------------------------- +# Single-target detection tests +# --------------------------------------------------------------------- + +FROM_SP_REQ = { + "pldl_db": 6.0, + "min_norm_pulse": 0.7, + "max_norm_pulse": 1.5, + "beam_comp_model": "simrad_lobe", + "max_beam_comp_db": 4.0, + "max_sd_minor_deg": 0.6, + "max_sd_major_deg": 0.6, +} + + +def _make_ds_from_Sp_minimal( + n_ping: int = 5, + n_range: int = 30, + background_db: float = -90.0, +) -> xr.Dataset: + """Create a minimal single-channel Sp dataset for detector tests.""" + ping_time = np.arange(n_ping) + range_sample = np.arange(n_range) + + coords = { + "ping_time": ping_time, + "range_sample": range_sample, + } + + sp = xr.DataArray( + np.full( + (n_ping, n_range), + background_db, + dtype=np.float64, + ), + dims=("ping_time", "range_sample"), + coords=coords, + name="Sp", + ) + + echo_range = xr.DataArray( + np.tile( + np.linspace(1.0, float(n_range), n_range)[None, :], + (n_ping, 1), + ), + dims=("ping_time", "range_sample"), + coords=coords, + name="echo_range", + ) + + angle_alongship = xr.DataArray( + np.zeros((n_ping, n_range), dtype=np.float64), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_alongship", + ) + + angle_athwartship = xr.DataArray( + np.zeros((n_ping, n_range), dtype=np.float64), + dims=("ping_time", "range_sample"), + coords=coords, + name="angle_athwartship", + ) + + return xr.Dataset( + data_vars={ + "Sp": sp, + "echo_range": echo_range, + "angle_alongship": angle_alongship, + "angle_athwartship": angle_athwartship, + "sound_absorption": xr.DataArray( + 0.003, + name="sound_absorption", + ), + "sample_interval": xr.DataArray( + 4e-5, + name="sample_interval", + ), + "tau_effective": xr.DataArray( + 8e-4, + name="tau_effective", + ), + "frequency_nominal": xr.DataArray( + 38000.0, + name="frequency_nominal", + ), + "beamwidth_alongship": xr.DataArray( + 7.0, + name="beamwidth_alongship", + ), + "beamwidth_athwartship": xr.DataArray( + 7.0, + name="beamwidth_athwartship", + ), + "angle_offset_alongship": xr.DataArray( + 0.0, + name="angle_offset_alongship", + ), + "angle_offset_athwartship": xr.DataArray( + 0.0, + name="angle_offset_athwartship", + ), + "beam_type": xr.DataArray( + np.int16(1), + name="beam_type", + ), + "channel": xr.DataArray( + "chan1", + name="channel", + ), + }, + coords=coords, + ) + + +# --------------------------------------------------------------------- +# Dispatcher validation +# --------------------------------------------------------------------- + + +@pytest.mark.unit +def test_detect_single_targets_unknown_method_raises(): + ds = _make_ds_from_Sp_minimal() + + with pytest.raises( + ValueError, + match="Unsupported single-target method", + ): + detect_single_targets( + ds, + method="__bad__", + params=dict(FROM_SP_REQ), + ) + + +@pytest.mark.unit +def test_detect_single_targets_fm_not_implemented_raises(): + ds = _make_ds_from_Sp_minimal() + + with pytest.raises( + NotImplementedError, + match="FM single-target detection", + ): + detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + waveform_mode="FM", + ) + + +@pytest.mark.unit +def test_detect_single_targets_invalid_waveform_mode_raises(): + ds = _make_ds_from_Sp_minimal() + + with pytest.raises( + ValueError, + match="waveform_mode must be 'CW' or 'FM'", + ): + detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + waveform_mode="BB", + ) + + +@pytest.mark.unit +def test_detect_single_targets_params_required_raises(): + ds = _make_ds_from_Sp_minimal() + + with pytest.raises( + ValueError, + match="No parameters given", + ): + detect_single_targets( + ds, + method="from_Sp", + params=None, + ) + + +@pytest.mark.unit +def test_detect_single_targets_beam_type_missing_raises(): + ds = _make_ds_from_Sp_minimal().drop_vars("beam_type") + + with pytest.raises( + ValueError, + match="beam_type variable is missing", + ): + detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + +@pytest.mark.unit +def test_detect_single_targets_beam_type_not_split_raises(): + ds = _make_ds_from_Sp_minimal() + ds["beam_type"] = xr.DataArray( + np.int16(2), + name="beam_type", + ) + + with pytest.raises( + ValueError, + match="Only split-beam data supported", + ): + detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + +# --------------------------------------------------------------------- +# Detector input validation +# --------------------------------------------------------------------- + + +@pytest.mark.unit +def test_detect_single_targets_from_Sp_missing_required_param_raises(): + ds = _make_ds_from_Sp_minimal() + params = dict(FROM_SP_REQ) + params.pop("pldl_db") + + with pytest.raises( + ValueError, + match="Missing required parameters", + ): + detect_single_targets( + ds, + method="from_Sp", + params=params, + ) + + +@pytest.mark.unit +def test_detect_single_targets_from_Sp_unknown_param_raises(): + ds = _make_ds_from_Sp_minimal() + params = dict(FROM_SP_REQ) + params["unknown_parameter"] = 42 + + with pytest.raises( + ValueError, + match="Unknown parameters", + ): + detect_single_targets( + ds, + method="from_Sp", + params=params, + ) + + +@pytest.mark.unit +@pytest.mark.parametrize( + "missing_variable", + [ + "Sp", + "echo_range", + "angle_alongship", + "angle_athwartship", + "tau_effective", + "sample_interval", + "sound_absorption", + "frequency_nominal", + ], +) +def test_detect_single_targets_from_Sp_missing_required_ds_var_raises( + missing_variable, +): + ds = _make_ds_from_Sp_minimal().drop_vars(missing_variable) + + with pytest.raises( + ValueError, + match="missing required variable", + ): + detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + +# --------------------------------------------------------------------- +# Empty output and schema +# --------------------------------------------------------------------- + + +@pytest.mark.unit +def test_detect_single_targets_from_Sp_returns_empty_schema(): + ds = _make_ds_from_Sp_minimal() + + out = detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + assert isinstance(out, xr.Dataset) + assert "single_target" in out.dims + assert out.sizes["single_target"] == 0 + + required_variables = ( + "channel", + "ping_time", + "range_sample", + "frequency_nominal", + "ping_index", + "iinf", + "isup", + "pulse_len_samples", + "norm_pulse_len", + "single_target_range", + "single_target_alongship_angle", + "single_target_athwartship_angle", + "single_target_alongship_angle_sd", + "single_target_athwartship_angle_sd", + "beam_comp_db", + "plike_peak", + ) + + for variable in required_variables: + assert variable in out + assert out[variable].dims == ("single_target",) + assert out[variable].sizes["single_target"] == 0 + + +@pytest.mark.unit +def test_detect_single_targets_concat_empty_is_stable(): + ds = _make_ds_from_Sp_minimal() + + out1 = detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + out2 = detect_single_targets( + ds, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + combined = xr.concat( + [out1, out2], + dim="single_target", + ) + + assert combined.sizes["single_target"] == 0 + + for variable in out1.data_vars: + assert combined[variable].dims == ("single_target",) + + +# --------------------------------------------------------------------- +# Echoview reference comparison +# --------------------------------------------------------------------- + +@pytest.mark.skip(reason="Echoview reference data will be added from the notebook") +@pytest.mark.integration +def test_detect_single_targets_from_Sp_matches_echoview_reference( + echoview_single_target_reference_path, +): + """ + Compare echopype target locations against a fixed Echoview export. + + The fixture should provide: + - source_Sp.nc + - echoview_single_targets.csv + + Required CSV columns: + ping_index, range_sample + """ + source_path = ( + echoview_single_target_reference_path / "source_Sp.nc" + ) + reference_path = ( + echoview_single_target_reference_path + / "echoview_single_targets.csv" + ) + + if not source_path.exists() or not reference_path.exists(): + pytest.skip( + "Echoview single-target reference files are unavailable." + ) + + ds_sp = xr.open_dataset(source_path) + reference = pd.read_csv(reference_path) + + out = detect_single_targets( + ds_sp, + method="from_Sp", + params=dict(FROM_SP_REQ), + ) + + actual_locations = { + (int(ping), int(sample)) + for ping, sample in zip( + out["ping_index"].values, + out["range_sample"].values, + ) + } + + reference_locations = { + (int(row.ping_index), int(row.range_sample)) + for row in reference.itertuples() + } + + assert actual_locations == reference_locations + +@pytest.mark.unit +def test_compute_TS_from_Sp_returns_expected_values(): + ds_sp = xr.Dataset( + data_vars={ + "Sp": ( + ("channel", "ping_time", "range_sample"), + np.array( + [ + [ + [-90.0, -80.0, -70.0], + [-60.0, -50.0, -40.0], + ] + ] + ), + ) + }, + coords={ + "channel": ["chan1"], + "ping_time": [0, 1], + "range_sample": [0, 1, 2], + }, + ) + + point_locations = xr.Dataset( + data_vars={ + "channel": ( + "single_target", + np.array(["chan1", "chan1"], dtype=object), + ), + "ping_index": ( + "single_target", + np.array([0, 1], dtype=np.int64), + ), + "range_sample": ( + "single_target", + np.array([2, 1], dtype=np.int64), + ), + "beam_comp_db": ( + "single_target", + np.array([1.5, 2.0], dtype=np.float64), + ), + }, + coords={ + "single_target": np.arange(2), + }, + ) + + out = compute_TS( + ds_sp, + point_locations=point_locations, + ) + + np.testing.assert_allclose( + out["uncompensated_TS"].values, + [-70.0, -50.0], + ) + + np.testing.assert_allclose( + out["compensated_TS"].values, + [-68.5, -48.0], + ) + + assert out["uncompensated_TS"].dims == ("single_target",) + assert out["compensated_TS"].dims == ("single_target",) + + +@pytest.mark.unit +def test_compute_TS_from_Sp_requires_point_locations(): + ds_sp = xr.Dataset( + { + "Sp": ( + ("channel", "ping_time", "range_sample"), + np.zeros((1, 2, 3)), + ) + }, + coords={ + "channel": ["chan1"], + "ping_time": [0, 1], + "range_sample": [0, 1, 2], + }, + ) + + with pytest.raises( + ValueError, + match="point_locations must be provided", + ): + compute_TS(ds_sp) + + +@pytest.mark.unit +def test_compute_TS_from_Sp_empty_targets(): + ds_sp = xr.Dataset( + { + "Sp": ( + ("channel", "ping_time", "range_sample"), + np.zeros((1, 2, 3)), + ) + }, + coords={ + "channel": ["chan1"], + "ping_time": [0, 1], + "range_sample": [0, 1, 2], + }, + ) + + point_locations = xr.Dataset( + data_vars={ + "channel": ( + "single_target", + np.array([], dtype=object), + ), + "ping_index": ( + "single_target", + np.array([], dtype=np.int64), + ), + "range_sample": ( + "single_target", + np.array([], dtype=np.int64), + ), + "beam_comp_db": ( + "single_target", + np.array([], dtype=np.float64), + ), + }, + coords={ + "single_target": np.array([], dtype=np.int64), + }, + ) + + out = compute_TS( + ds_sp, + point_locations=point_locations, + ) + + assert out.sizes["single_target"] == 0 + assert "uncompensated_TS" in out + assert "compensated_TS" in out \ No newline at end of file