From ad487687d6062e552bfd6a04a69c6c411b667fed Mon Sep 17 00:00:00 2001 From: Wu-Jung Lee Date: Thu, 22 Dec 2022 09:00:20 -0800 Subject: [PATCH] Test merge split beam angle (#15) * create general structure of add_splitbeam_angle function * create simrad.py and start working on a check of the waveform and encode modes * rearange and improve logic for checking the encode and waveform modes * finish documenting functions in simrad and finalizing their outputs, incorporate new waveform/encode mode check into calibrate_ek and consolidate api * allow for 2 beam groups for waveform_mode=BB and add check_waveform_encode_mode to compute_Sv and compute_TS * move core split-beam angle functions to split_beam_angle.py, finish implementation of _get_splitbeam_angle_power_CW, _add_splitbeam_angle_to_ds, add an empty test suite for echodata.simrad, and start working on test_add_splitbeam_angle * make check_waveform_encode_mode return only the echodata group name corresponding to the encode mode * change test_add_splitbeam_angle so that it takes a list of echoview mat paths, construct a function that creates a numpy array from mat paths, and select echodata beam_group data in add_splitbeam_angle based on ds channel dims * add pulse_compression input to add_splitbeam_angle, implement (CW, complex) and (BB, complex, no pulse compression) options, and create tests for them * modify test_add_splitbeam_angle pytest parameters to account for data location, finish docstring for add_splitbeam_angle, drop beam wherever necessary * add comment at the top of test_consolidate.py specifying where the split-beam angle tests data exists * creat apply_mask function without implementing it * create general structure implementation of apply_mask * create first simple test for apply_mask with mock Dataset creation * create apply_mask paramter test for when mask is a list of DataArrays * change the function name validate_source_ds to validate_source_ds_da and allow it to handle a DataArray input, modify existing tests for validate_source_ds * create routines that check the input of apply_mask and form the inputs for the core method in apply_mask * allow get_mock_source_ds_apply_mask to produced delayed variables and add associated test in test_apply_mask * add code and parameters to test_apply_mask to test for a mask that is provided as a path * add tests for different fill_value types and accout for type(fill_value)=DataArray conflict with xr.where * correct docstring in _check_mode_input_without_data Co-authored-by: Wu-Jung Lee * raise ValueError instead of RuntimeError in _check_mode_input_without_data * remove TODO in _check_mode_input_with_data_EK60 * improve grammar in error output in _check_mode_input_with_data_EK80 Co-authored-by: Wu-Jung Lee * remove check for EK80-like sensors in _check_mode_input_with_data_EK80 as it is not necessary * improve docstring function description for add_splitbeam_angle Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove States in pulse_compression docstring description Co-authored-by: Wu-Jung Lee * improve logic associated with split-beam angle function calls in add_splitbeam_angle Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove remnant else statement caused by GH commit in add_splitbeam_angle * rename _check_mode_input_without_data to _check_input_args_combination * add reminder to add a test for _check_input_args_combination * change _check_mode_input_with_data_EK60/EK80 to _retrieve_correct_beam_group_EK60/EK80 * rename check_waveform_encode_mode to retrieve_correct_beam_group * finish constructing test_check_input_args_combination * remove reference to MVBS for ds input in add_splitbeam_angle function Co-authored-by: Wu-Jung Lee * rename ds to source_Sv in add_splitbeam_angle * allow source_Sv to be a Dataset, str, or pathlib.Path * add return_dataset argument to add_splitb_ang and modify _add_splitbeam_angle_to_ds so that it will write the data to a zarr or netcdf, if necessary * remove auto chunking when opening a Dataset since netcdf files have issues with this and free resources linked to ds in _add_splitbeam_angle_to_ds as netcdf also has an issue if the file remains open * modify test_add_splitbeam_angle to allow writing ds_Sv to a file and add associated test parameters * raise not implemented error if source_Sv corresponds to MVBS * implement split-beam angle calculation for waveform_mode=BB, encode_mode=complex, and pulse_compression=False, and add corresponding test * Improve documentation by clarifying statements Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * remove whitespace caused by GH commit * replace RuntimeError with ValueError * replace chunks=auto with chunks={} * add default value to apply_mask input var_name and modify documentation accordingly * Fix meta_source_filenames bug and enable (meta)source_filenames appending of path and list (#908) * Expand prov.source_files_vars to support path sequences that mix a str/path and another sequence * For EK80, ES70, ES80, EA640, conversion was inserting an unnecessary, empty meta_source_filenames variable * Add more comprehensive and readable prov source-files type hints; rename _source_files * Fix prov type hint bug with Py 3.8 * Add unit test for prov _sanitize_source_files; plus small fixes to np.ndarray type references in prov * fix validate_source_ds to validate_source_ds_da * improve Notes section in docstring for add_splitbeam_angle Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * clarify when the split-beam angle data potentially exist Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fix long lines * Add clarification in NotImplemented error statement Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * split up line in NotImplemented error statement * rename function compute_split_beam_beamtype1 to _e2f * put radian to degree conversion on one line Co-authored-by: Wu-Jung Lee * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * change sin to arcsin Co-authored-by: Wu-Jung Lee * remove unnecessary statement Co-authored-by: Wu-Jung Lee * replace RuntimeError with ValueError * change theta_fc --> theta and phi_fc --> phi * allow splitbeam calculation for power and CW options when beam_type is not equal to zero * replace core functions in _get_splitbeam_angle_complex_CW with get_splitbeam and then remove them * shorten function names, minor docstring tweaks * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add TODO to test * fix bug in logic for at least 1 channel being split-beam Co-authored-by: b-reyes Co-authored-by: b-reyes <53541061+b-reyes@users.noreply.github.com> Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Emilio Mayorga --- echopype/calibrate/calibrate_ek.py | 105 +--- echopype/consolidate/__init__.py | 4 +- echopype/consolidate/api.py | 177 ++++++- echopype/consolidate/split_beam_angle.py | 501 ++++++++++++++++++ echopype/echodata/simrad.py | 229 ++++++++ .../tests/consolidate/test_consolidate.py | 158 +++++- .../tests/echodata/test_echodata_simrad.py | 69 +++ 7 files changed, 1159 insertions(+), 84 deletions(-) create mode 100644 echopype/consolidate/split_beam_angle.py create mode 100644 echopype/echodata/simrad.py create mode 100644 echopype/tests/echodata/test_echodata_simrad.py diff --git a/echopype/calibrate/calibrate_ek.py b/echopype/calibrate/calibrate_ek.py index 06bdb0b2ac..11b03d7e20 100644 --- a/echopype/calibrate/calibrate_ek.py +++ b/echopype/calibrate/calibrate_ek.py @@ -3,6 +3,7 @@ from scipy import signal from ..echodata import EchoData +from ..echodata.simrad import retrieve_correct_beam_group from ..utils import uwa from ..utils.log import _init_logger from .calibrate_base import CAL_PARAMS, CalibrateBase @@ -134,19 +135,16 @@ def get_cal_params(self, cal_params, waveform_mode, encode_mode): else beam["equivalent_beam_angle"] ) - def _cal_power(self, cal_type, use_beam_power=False) -> xr.Dataset: + def _cal_power(self, cal_type: str, power_ed_group: str = None) -> xr.Dataset: """Calibrate power data from EK60 and EK80. Parameters ---------- - cal_type : str + cal_type: str 'Sv' for calculating volume backscattering strength, or 'TS' for calculating target strength - use_beam_power : bool - whether to use beam_power. - If ``True`` use ``echodata["Sonar/Beam_group2"]``; - if ``False`` use ``echodata["Sonar/Beam_group1"]``. - Note ``echodata["Sonar/Beam_group2"]`` could only exist for EK80 data. + power_ed_group: + The ``EchoData`` beam group path containing the power data Returns ------- @@ -154,10 +152,7 @@ def _cal_power(self, cal_type, use_beam_power=False) -> xr.Dataset: The calibrated dataset containing Sv or TS """ # Select source of backscatter data - if use_beam_power: - beam = self.echodata["Sonar/Beam_group2"] - else: - beam = self.echodata["Sonar/Beam_group1"] + beam = self.echodata[power_ed_group] # Harmonize time coordinate between Beam_groupX data and env_params for p in self.env_params.keys(): @@ -284,10 +279,16 @@ def get_env_params(self, **kwargs): ) def compute_Sv(self, **kwargs): - return self._cal_power(cal_type="Sv") + power_ed_group = retrieve_correct_beam_group( + echodata=self.echodata, waveform_mode="CW", encode_mode="power", pulse_compression=False + ) + return self._cal_power(cal_type="Sv", power_ed_group=power_ed_group) def compute_TS(self, **kwargs): - return self._cal_power(cal_type="TS") + power_ed_group = retrieve_correct_beam_group( + echodata=self.echodata, waveform_mode="CW", encode_mode="power", pulse_compression=False + ) + return self._cal_power(cal_type="TS", power_ed_group=power_ed_group) class CalibrateEK80(CalibrateEK): @@ -883,78 +884,22 @@ def _compute_cal(self, cal_type, waveform_mode, encode_mode) -> xr.Dataset: xr.Dataset An xarray Dataset containing either Sv or TS. """ - # Raise error for wrong inputs - if waveform_mode not in ("BB", "CW"): - raise ValueError( - "Input waveform_mode not recognized! " - "waveform_mode must be either 'BB' or 'CW' for EK80 data." - ) - if encode_mode not in ("complex", "power"): - raise ValueError( - "Input encode_mode not recognized! " - "encode_mode must be either 'complex' or 'power' for EK80 data." - ) + + power_ed_group = retrieve_correct_beam_group( + echodata=self.echodata, + waveform_mode=waveform_mode, + encode_mode=encode_mode, + pulse_compression=False, + ) # Set flag_complex # - True: complex cal # - False: power cal - # BB: complex only, CW: complex or power + flag_complex = False if waveform_mode == "BB": - if encode_mode == "power": # BB waveform forces to collect complex samples - raise ValueError("encode_mode='power' not allowed when waveform_mode='BB'!") flag_complex = True - else: # waveform_mode="CW" - if encode_mode == "complex": - flag_complex = True - else: - flag_complex = False - - # Raise error when waveform_mode and actual recording mode do not match - # This simple check is only possible for BB-only data, - # since for data with both BB and CW complex samples, - # frequency_start will exist in echodata["Sonar/Beam_group1"] for the BB channels - if waveform_mode == "BB" and "frequency_start" not in self.echodata["Sonar/Beam_group1"]: - raise ValueError("waveform_mode='BB' but broadband data not found!") - - # Set use_beam_power - # - True: use self.echodata["Sonar/Beam_group2"] for cal - # - False: use self.echodata["Sonar/Beam_group1"] for cal - use_beam_power = False - - # Warn user about additional data in the raw file if another type exists - # When both power and complex samples exist: - # complex samples will be stored in echodata["Sonar/Beam_group1"] - # power samples will be stored in echodata["Sonar/Beam_group2"] - # When only one type of samples exist, - # all samples with be stored in echodata["Sonar/Beam_group1"] - if self.echodata["Sonar/Beam_group2"] is not None: # both power and complex samples exist - # If both beam and beam_power groups exist, - # this means that CW data are encoded as power samples and in beam_power group - if waveform_mode == "CW" and encode_mode == "complex": - raise ValueError("File does not contain CW complex samples") - - if encode_mode == "power": - use_beam_power = True # switch source of backscatter data - logger.info( - "Only power samples are calibrated, but complex samples also exist in the raw data file!" # noqa - ) - else: - logger.info( - "Only complex samples are calibrated, but power samples also exist in the raw data file!" # noqa - ) - else: # only power OR complex samples exist - if ( - "backscatter_i" in self.echodata["Sonar/Beam_group1"].variables - ): # data contain only complex samples - if encode_mode == "power": - raise TypeError( - "File does not contain power samples! Use encode_mode='complex'" - ) # user selects the wrong encode_mode - else: # data contain only power samples - if encode_mode == "complex": - raise TypeError( - "File does not contain complex samples! Use encode_mode='power'" - ) # user selects the wrong encode_mode + elif encode_mode == "complex": + flag_complex = True # Compute Sv if flag_complex: @@ -964,7 +909,7 @@ def _compute_cal(self, cal_type, waveform_mode, encode_mode) -> xr.Dataset: else: # Power samples only make sense for CW mode data self.compute_range_meter(waveform_mode="CW", encode_mode=encode_mode) - ds_cal = self._cal_power(cal_type=cal_type, use_beam_power=use_beam_power) + ds_cal = self._cal_power(cal_type=cal_type, power_ed_group=power_ed_group) return ds_cal diff --git a/echopype/consolidate/__init__.py b/echopype/consolidate/__init__.py index d26813ab89..acca09fd80 100644 --- a/echopype/consolidate/__init__.py +++ b/echopype/consolidate/__init__.py @@ -1,3 +1,3 @@ -from .api import add_depth, add_location, swap_dims_channel_frequency +from .api import add_depth, add_location, add_splitbeam_angle, swap_dims_channel_frequency -__all__ = ["swap_dims_channel_frequency", "add_depth", "add_location"] +__all__ = ["swap_dims_channel_frequency", "add_depth", "add_location", "add_splitbeam_angle"] diff --git a/echopype/consolidate/api.py b/echopype/consolidate/api.py index de81481bc0..d80d0e3f32 100644 --- a/echopype/consolidate/api.py +++ b/echopype/consolidate/api.py @@ -1,10 +1,20 @@ import datetime -from typing import Optional +import pathlib +from typing import Optional, Union import numpy as np import xarray as xr from ..echodata import EchoData +from ..echodata.simrad import retrieve_correct_beam_group +from ..utils.io import validate_source_ds_da +from .split_beam_angle import ( + add_angle_to_ds, + get_angle_complex_BB_nopc, + get_angle_complex_BB_pc, + get_angle_complex_CW, + get_angle_power_CW, +) def swap_dims_channel_frequency(ds: xr.Dataset) -> xr.Dataset: @@ -174,3 +184,168 @@ def sel_interp(var): interp_ds["longitude"] = interp_ds["longitude"].assign_attrs({"history": history}) return interp_ds.drop_vars("time1") + + +def add_splitbeam_angle( + source_Sv: Union[xr.Dataset, str, pathlib.Path], + echodata: EchoData, + waveform_mode: str, + encode_mode: str, + pulse_compression: bool = False, + storage_options: dict = {}, + return_dataset: bool = True, +) -> xr.Dataset: + """ + Add split-beam (alongship/athwartship) angles into the Sv dataset. + This function calculates the alongship/athwartship angle using data stored + in the beam groups of ``echodata``. In cases when angle data does not already exist + or cannot be computed from the data, an error is issued and no angle variables are + added to the dataset. + + Parameters + ---------- + source_Sv: xr.Dataset or str or pathlib.Path + The Sv Dataset or path to a file containing the Sv Dataset, which will have the + split-beam angle data added to it + echodata: EchoData + An ``EchoData`` object holding the raw data + waveform_mode : {"CW", "BB"} + Type of transmit waveform + + - ``"CW"`` for narrowband transmission, + returned echoes recorded either as complex or power/angle samples + - ``"BB"`` for broadband transmission, + returned echoes recorded as complex samples + + encode_mode : {"complex", "power"} + Type of encoded return echo data + + - ``"complex"`` for complex samples + - ``"power"`` for power/angle samples, only allowed when + the echosounder is configured for narrowband transmission + pulse_compression: bool, False + Whether pulse compression should be used (only valid for + ``waveform_mode="BB"`` and ``encode_mode="complex"``) + storage_options: dict, default={} + Any additional parameters for the storage backend, corresponding to the + path provided for ``source_Sv`` + return_dataset: bool, default=True + If True, ``source_Sv`` with the split-beam angle data added to it + will be returned, else it will not be returned. A value of ``False`` + is useful in the situation where ``source_Sv`` is a path and the user + only wants to write the split-beam angle data to the path provided. + + Returns + ------- + xr.Dataset or None + If ``return_dataset=False``, nothing will be returned. If ``return_dataset=True`` + either the input dataset ``source_Sv`` or a lazy-loaded Dataset (obtained from + the path provided by ``source_Sv``) with the split-beam angle data added + will be returned. + + Raises + ------ + ValueError + If ``echodata`` has a sonar model that is not analogous to either EK60 or EK80 + ValueError + If the input ``source_Sv`` does not have a ``channel`` dimension + ValueError + If ``source_Sv`` does not have appropriate dimension lengths in + comparison to ``echodata`` data + ValueError + If the provided ``waveform_mode``, ``encode_mode``, and ``pulse_compression`` are not valid + NotImplementedError + If an unknown ``beam_type`` is encountered during the split-beam calculation + + Notes + ----- + Split-beam angle data potentially exist for the following echosounders depending on + the instrument configuration and recording setting: + + - Simrad EK60 echosounder paired with split-beam transducers and + configured to store angle data + - Simrad EK80 echosounder paired with split-beam transducers and + configured to store angle data + + In most cases where the type of samples collected by the echosounder (power/angle + samples or complex samples) and the transmit waveform (broadband or narrowband) + are identical across all channels, the channels existing in ``source_Sv`` and ` + `echodata`` will be identical. If this is not the case, only angle data corresponding + to channels existing in ``source_Sv`` will be added. + + For EK80 broadband data, the split-beam angles can be estimated from the complex data. + The current implementation generates angles estimated *without* applying pulse compression. + Estimating the angle with pulse compression will be added in the near future. + """ + + # ensure that echodata was produced by EK60 or EK80-like sensors + if echodata.sonar_model not in ["EK60", "ES70", "EK80", "ES80", "EA640"]: + raise ValueError( + "The sonar model that produced echodata does not have split-beam " + "transducers, split-beam angles cannot be added to source_Sv!" + ) + + # validate the source_Sv type or path (if it is provided) + source_Sv, file_type = validate_source_ds_da(source_Sv, storage_options) + + # initialize source_Sv_path + source_Sv_path = None + + if isinstance(source_Sv, str): + + # store source_Sv path so we can use it to write to later + source_Sv_path = source_Sv + + # TODO: In the future we can improve this by obtaining the variable names, channels, + # and dimension lengths directly from source_Sv using zarr or netcdf4. This would + # prevent the unnecessary loading in of the coordinates, which the below statement does. + # open up Dataset using source_Sv path + source_Sv = xr.open_dataset(source_Sv, engine=file_type, chunks={}, **storage_options) + + # raise not implemented error if source_Sv corresponds to MVBS + if source_Sv.attrs["processing_function"] == "preprocess.compute_MVBS": + raise NotImplementedError("Adding split-beam data to MVBS has not been implemented!") + + # check that the appropriate waveform and encode mode have been given + # and obtain the echodata group path corresponding to encode_mode + encode_mode_ed_group = retrieve_correct_beam_group( + echodata, waveform_mode, encode_mode, pulse_compression + ) + + # check that source_Sv at least has a channel dimension + if "channel" not in source_Sv.variables: + raise ValueError("The input source_Sv Dataset must have a channel dimension!") + + # set ds_beam, select the same channels that are in source_Sv + ds_beam = echodata[encode_mode_ed_group].sel(channel=source_Sv.channel.values) + + # fail if source_Sv and ds_beam do not have the same lengths + # for ping_time, range_sample, and channel + same_dim_lens = [ + ds_beam.dims[dim] == source_Sv.dims[dim] for dim in ["channel", "ping_time", "range_sample"] + ] + if not same_dim_lens: + raise ValueError( + "Input source_Sv does not have the same dimension lengths as all dimensions in ds_beam!" + ) + + # obtain split-beam angles from + # CW mode data + if waveform_mode == "CW": + if encode_mode == "power": # power data + theta, phi = get_angle_power_CW(ds_beam=ds_beam) + else: # complex data + theta, phi = get_angle_complex_CW(ds_beam=ds_beam) + # BB mode data + else: + if pulse_compression: # with pulse compression + theta, phi = get_angle_complex_BB_pc(ds_beam=ds_beam) + else: # without pulse compression + theta, phi = get_angle_complex_BB_nopc(ds_beam=ds_beam, ed=echodata) + + # add theta and phi to source_Sv input + source_Sv = add_angle_to_ds( + theta, phi, source_Sv, return_dataset, source_Sv_path, file_type, storage_options + ) + + return source_Sv diff --git a/echopype/consolidate/split_beam_angle.py b/echopype/consolidate/split_beam_angle.py new file mode 100644 index 0000000000..96796c8519 --- /dev/null +++ b/echopype/consolidate/split_beam_angle.py @@ -0,0 +1,501 @@ +""" +Contains functions necessary to compute the split-beam (alongship/athwartship) +angles and add them to a Dataset. +""" +from typing import List, Optional, Tuple + +import numpy as np +import xarray as xr + +from ..echodata import EchoData + + +def get_angle_power_CW(ds_beam: xr.Dataset) -> Tuple[xr.Dataset, xr.Dataset]: + """ + Obtains the split-beam angle data from power encoded data with CW waveform. + + Parameters + ---------- + ds_beam: xr.Dataset + An ``EchoData`` beam group containing angle information needed for + split-beam angle calculation + + Returns + ------- + theta: xr.Dataset + The calculated split-beam alongship angle + phi: xr.Dataset + The calculated split-beam athwartship angle + + Raises + ------ + NotImplementedError + If all ``beam_type`` values are not equal to 1 + + Notes + ----- + Can be used on both EK60 and EK80 data + + Computation done for ``beam_type=1``: + ``physical_angle = ((raw_angle * 180 / 128) / sensitivity) - offset`` + """ + + # raw_angle scaling constant + conversion_const = 180.0 / 128.0 + + def _e2f(angle_type: str) -> xr.Dataset: + """Convert electric angle to physical angle for split-beam data""" + return ( + conversion_const + * ds_beam[f"angle_{angle_type}"] + / ds_beam[f"angle_sensitivity_{angle_type}"] + - ds_beam[f"angle_offset_{angle_type}"] + ) + + # add split-beam angle if at least one channel is split-beam + # in the case when some channels are split-beam and some single-beam + # the single-beam channels will be all NaNs and _e2f would run through and output NaNs + if not np.all(ds_beam["beam_type"].data == 0): + + # obtain split-beam alongship angle + theta = _e2f(angle_type="alongship") + + # obtain split-beam athwartship angle + phi = _e2f(angle_type="athwartship") + + else: + raise ValueError( + "Computing physical split-beam angle is only available for data " + "from split-beam transducers!" + ) + + # drop the beam dimension in theta and phi, if it exists + if "beam" in theta.dims: + theta = theta.drop("beam").squeeze(dim="beam") + phi = phi.drop("beam").squeeze(dim="beam") + + return theta, phi + + +def get_angle_complex_CW(ds_beam: xr.Dataset) -> Tuple[xr.DataArray, xr.DataArray]: + """ + Obtains the split-beam angle data from complex encoded data with CW waveform. + + Parameters + ---------- + ds_beam: xr.Dataset + An ``EchoData`` beam group containing angle information needed for + split-beam angle calculation + + Returns + ------- + theta: xr.Dataset + The calculated split-beam alongship angle + phi: xr.Dataset + The calculated split-beam athwartship angle + """ + + # ensure that the beam_type is appropriate for calculation + if np.all(ds_beam["beam_type"].data == 1): + + # get complex representation of backscatter + backscatter = ds_beam["backscatter_r"] + 1j * ds_beam["backscatter_i"] + + # get angle sensitivity alongship and athwartship + angle_sensitivity_alongship = ds_beam["angle_sensitivity_alongship"].isel( + ping_time=0, beam=0 + ) + angle_sensitivity_athwartship = ds_beam["angle_sensitivity_athwartship"].isel( + ping_time=0, beam=0 + ) + + # get angle offset alongship and athwartship + angle_offset_alongship = ds_beam["angle_offset_alongship"].isel(ping_time=0, beam=0) + angle_offset_athwartship = ds_beam["angle_offset_athwartship"].isel(ping_time=0, beam=0) + + # obtain the split-beam angle data + theta, phi = _compute_angle_from_complex( + bs=backscatter, + beam_type=1, + sens=[angle_sensitivity_alongship, angle_sensitivity_athwartship], + offset=[angle_offset_alongship, angle_offset_athwartship], + ) + + else: + raise NotImplementedError("Computing split-beam angle is only available for beam_type=1!") + + # drop the beam dimension in theta and phi, if it exists + if "beam" in theta.coords: + theta = theta.drop_vars("beam") + phi = phi.drop("beam") + + return theta, phi + + +def _get_interp_offset( + param: str, chan_id: str, freq_center: xr.DataArray, ed: EchoData +) -> np.ndarray: + """ + Obtains an angle offset by first interpolating the + ``angle_offset_alongship`` or ``angle_offset_athwartship`` + data found in the ``Vendor_specific`` group and then + selecting the offset corresponding to the center frequency + value for ``channel=chan_id``. + + Parameters + ---------- + param: {"angle_offset_alongship", "angle_offset_athwartship"} + The angle offset data to select in the ``Vendor_specific`` group + chan_id: str + The channel used to select the center frequency value + freq_center: xr.DataArray + A DataArray filled with center frequency values with coordinate ``channel`` + ed: EchoData + An ``EchoData`` object holding the raw data + + Returns + ------- + np.ndarray + Array filled with the requested angle offset values + """ + + freq_wanted = freq_center.sel(channel=chan_id) + return ( + ed["Vendor_specific"][param].sel(cal_channel_id=chan_id).interp(cal_frequency=freq_wanted) + ).values + + +def _get_offset( + ds_beam: xr.Dataset, fc: xr.DataArray, freq_nominal: xr.DataArray, ed: EchoData +) -> Tuple[xr.DataArray, xr.DataArray]: + """ + Obtains the alongship and athwartship angle offsets. + + Parameters + ---------- + ds_beam: xr.Dataset + The dataset corresponding to a beam group + fc: xr.DataArray + Array corresponding to the center frequency + freq_nominal: xr.DataArray + Array of frequency nominal values + ed: EchoData + An ``EchoData`` object holding the raw data + + Returns + ------- + offset_along: xr.DataArray + Array corresponding to the angle alongship offset + offset_athwart: xr.DataArray + Array corresponding to the angle athwartship offset + """ + + # initialize lists that will hold offsets + offset_along = [] + offset_athwart = [] + + # obtain the offsets for each channel + for ch in fc["channel"].values: + if ch in ed["Vendor_specific"]["cal_channel_id"]: + # calculate offsets using Vendor_specific values + offset_along.append( + _get_interp_offset( + param="angle_offset_alongship", chan_id=ch, freq_center=fc, ed=ed + ) + ) + offset_athwart.append( + _get_interp_offset( + param="angle_offset_athwartship", chan_id=ch, freq_center=fc, ed=ed + ) + ) + else: + # calculate offsets using data in ds_beam + offset_along.append( + ds_beam["angle_offset_alongship"].sel(channel=ch).isel(ping_time=0, beam=0) + * fc.sel(channel=ch) + / freq_nominal.sel(channel=ch) + ) + offset_athwart.append( + ds_beam["angle_offset_athwartship"].sel(channel=ch).isel(ping_time=0, beam=0) + * fc.sel(channel=ch) + / freq_nominal.sel(channel=ch) + ) + + # construct offset DataArrays from lists + offset_along = xr.DataArray( + offset_along, coords={"channel": fc["channel"], "ping_time": fc["ping_time"]} + ) + offset_athwart = xr.DataArray( + offset_athwart, coords={"channel": fc["channel"], "ping_time": fc["ping_time"]} + ) + return offset_along, offset_athwart + + +def _compute_angle_from_complex( + bs: xr.Dataset, beam_type: int, sens: List[xr.DataArray], offset: List[xr.DataArray] +): + """ + Obtains the split-beam angle data alongship and athwartship + using data from a single channel. + + Parameters + ---------- + bs: xr.Dataset + Complex representation of backscatter + beam_type: int + The type of beam being considered + sens: list of xr.DataArray + A list of length two where the first element corresponds to the + angle sensitivity alongship and the second corresponds to the + angle sensitivity athwartship + offset: list of xr.DataArray + A list of length two where the first element corresponds to the + angle offset alongship and the second corresponds to the + angle offset athwartship + + Returns + ------- + theta: xr.Dataset + The calculated split-beam alongship angle for a specific channel + phi: xr.Dataset + The calculated split-beam athwartship angle for a specific channel + + Notes + ----- + This function should only be used for data with complex backscatter. + """ + + # 4-sector transducer + if beam_type == 1: + + bs_fore = (bs.isel(beam=2) + bs.isel(beam=3)) / 2 # forward + bs_aft = (bs.isel(beam=0) + bs.isel(beam=1)) / 2 # aft + bs_star = (bs.isel(beam=0) + bs.isel(beam=3)) / 2 # starboard + bs_port = (bs.isel(beam=1) + bs.isel(beam=2)) / 2 # port + + bs_theta = bs_fore * np.conj(bs_aft) + bs_phi = bs_star * np.conj(bs_port) + theta = np.arctan2(np.imag(bs_theta), np.real(bs_theta)) / np.pi * 180 + phi = np.arctan2(np.imag(bs_phi), np.real(bs_phi)) / np.pi * 180 + + # 3-sector transducer with or without center element + elif beam_type in [17, 49, 65, 81]: + # 3-sector + if beam_type == 17: + bs_star = bs.isel(beam=0) + bs_port = bs.isel(beam=1) + bs_fore = bs.isel(beam=2) + else: + # 3-sector + 1 center element + bs_star = (bs.isel(beam=0) + bs.isel(beam=3)) / 2 + bs_port = (bs.isel(beam=1) + bs.isel(beam=3)) / 2 + bs_fore = (bs.isel(beam=2) + bs.isel(beam=3)) / 2 + + bs_fac1 = bs_fore * np.conj(bs_star) + bs_fac2 = bs_fore * np.conj(bs_port) + fac1 = np.arctan2(np.imag(bs_fac1), np.real(bs_fac1)) / np.pi * 180 + fac2 = np.arctan2(np.imag(bs_fac2), np.real(bs_fac2)) / np.pi * 180 + + theta = (fac1 + fac2) / np.sqrt(3) + phi = fac2 - fac1 + + # EC150–3C + elif beam_type == 97: + raise NotImplementedError + + else: + raise ValueError("beam_type not recognized!") + + theta = theta / sens[0] - offset[0] + phi = phi / sens[1] - offset[1] + + return theta, phi + + +def get_angle_complex_BB_nopc( + ds_beam: xr.Dataset, ed: EchoData +) -> Tuple[xr.DataArray, xr.DataArray]: + """ + Obtains the split-beam angle data from complex samples from broadband transmit signals + without pulse compression. + + Parameters + ---------- + ds_beam: xr.Dataset + An ``EchoData`` beam group containing angle information needed for + split-beam angle calculation + ed: EchoData + An ``EchoData`` object holding the raw data + + Returns + ------- + theta: xr.Dataset + The calculated split-beam alongship angle + phi: xr.Dataset + The calculated split-beam athwartship angle + """ + + # nominal frequency [Hz] + freq_nominal = ds_beam["frequency_nominal"] + + # calculate center frequency + freq_center = (ds_beam["frequency_start"] + ds_beam["frequency_end"]).isel(beam=0) / 2 + + # obtain the angle alongship and athwartship offsets + offset_along, offset_athwart = _get_offset( + ds_beam=ds_beam, fc=freq_center, freq_nominal=freq_nominal, ed=ed + ) + + # obtain the angle sensitivity values alongship and athwartship + sens_along = ds_beam["angle_sensitivity_alongship"].isel(beam=0) * freq_center / freq_nominal + sens_athwart = ( + ds_beam["angle_sensitivity_athwartship"].isel(beam=0) * freq_center / freq_nominal + ) + + # get complex representation of backscatter + backscatter = ds_beam["backscatter_r"] + 1j * ds_beam["backscatter_i"] + + # initialize list that will hold split-beam angle data for each channel + theta_channels = [] + phi_channels = [] + + # obtain the split-beam angle data for each channel + for chan_id in backscatter.channel.values: + theta, phi = _compute_angle_from_complex( + bs=backscatter.sel(channel=chan_id), + beam_type=int(ds_beam["beam_type"].sel(channel=chan_id).isel(ping_time=0)), + sens=[sens_along.sel(channel=chan_id), sens_athwart.sel(channel=chan_id)], + offset=[offset_along.sel(channel=chan_id), offset_athwart.sel(channel=chan_id)], + ) + + theta_channels.append(theta) + phi_channels.append(phi) + + # collect and construct final DataArrays for split-beam angle data + theta = xr.DataArray( + data=theta_channels, + coords={ + "channel": backscatter.channel, + "ping_time": theta_channels[0].ping_time, + "range_sample": theta_channels[0].range_sample, + }, + ) + + phi = xr.DataArray( + data=phi_channels, + coords={ + "channel": backscatter.channel, + "ping_time": phi_channels[0].ping_time, + "range_sample": phi_channels[0].range_sample, + }, + ) + + return theta, phi + + +def get_angle_complex_BB_pc(ds_beam: xr.Dataset) -> Tuple[xr.DataArray, xr.DataArray]: + """ + Obtains the split-beam angle data from complex samples from broadband transmit signals + after pulse compression. + + Parameters + ---------- + ds_beam: xr.Dataset + An ``EchoData`` beam group containing angle information needed for + split-beam angle calculation + + Returns + ------- + theta: xr.Dataset + The calculated split-beam alongship angle + phi: xr.Dataset + The calculated split-beam athwartship angle + """ + + # TODO: make sure to check that the appropriate beam_type is being used + raise NotImplementedError( + "Obtaining the split-beam angle data using pulse compressed " + "backscatter has not been implemented!" + ) + + return xr.DataArray(), xr.DataArray() + + +def add_angle_to_ds( + theta: xr.Dataset, + phi: xr.Dataset, + ds: xr.Dataset, + return_dataset: bool, + source_ds_path: Optional[str] = None, + file_type: Optional[str] = None, + storage_options: dict = {}, +) -> Optional[xr.Dataset]: + """ + Adds the split-beam angle data to the provided input ``ds``. + + Parameters + ---------- + theta: xr.Dataset + The calculated split-beam alongship angle + phi: xr.Dataset + The calculated split-beam athwartship angle + ds: xr.Dataset + The Dataset that ``theta`` and ``phi`` will be added to + return_dataset: bool + Whether a dataset will be returned or not + source_ds_path: str, optional + The path to the file corresponding to ``ds``, if it exists + file_type: {"netcdf4", "zarr"}, optional + The file type corresponding to ``source_ds_path`` + storage_options: dict, default={} + Any additional parameters for the storage backend, corresponding to the + path ``source_ds_path`` + + Returns + ------- + xr.Dataset or None + If ``return_dataset=False``, nothing will be returned. If ``return_dataset=True`` + either the input dataset ``ds`` or a lazy-loaded Dataset (obtained from + the path provided by ``source_ds_path``) with the split-beam angle data added + will be returned. + """ + + # TODO: do we want to add anymore attributes to these variables? + # add appropriate attributes to theta and phi + theta.attrs["long_name"] = "split-beam alongship angle" + phi.attrs["long_name"] = "split-beam athwartship angle" + + if source_ds_path is not None: + + # put the variables into a Dataset, so they can be written at the same time + # add ds attributes to splitb_ds since they will be overwritten by to_netcdf/zarr + splitb_ds = xr.Dataset( + data_vars={"angle_alongship": theta, "angle_athwartship": phi}, + coords=theta.coords, + attrs=ds.attrs, + ) + + # release any resources linked to ds (necessary for to_netcdf) + ds.close() + + # write the split-beam angle data to the provided path + if file_type == "netcdf4": + splitb_ds.to_netcdf(path=source_ds_path, mode="a", **storage_options) + else: + splitb_ds.to_zarr(store=source_ds_path, mode="a", **storage_options) + + else: + + # add the split-beam angles to the provided Dataset + ds["angle_alongship"] = theta + ds["angle_athwartship"] = phi + + if return_dataset and (source_ds_path is not None): + + # open up and return Dataset in source_ds_path + return xr.open_dataset(source_ds_path, engine=file_type, chunks={}, **storage_options) + + elif return_dataset: + + # return input dataset with split-beam angle data + return ds diff --git a/echopype/echodata/simrad.py b/echopype/echodata/simrad.py new file mode 100644 index 0000000000..195f68e1f5 --- /dev/null +++ b/echopype/echodata/simrad.py @@ -0,0 +1,229 @@ +""" +Contains functions that are specific to Simrad echo sounders +""" +from typing import Optional, Tuple + +from .echodata import EchoData + + +def _check_input_args_combination( + waveform_mode: str, encode_mode: str, pulse_compression: bool +) -> None: + """ + Checks that the ``waveform_mode`` and ``encode_mode`` have + the correct values and that the combination of input arguments are valid, without + considering the actual data. + + Parameters + ---------- + waveform_mode: str + Type of transmit waveform + encode_mode: str + Type of encoded return echo data + pulse_compression: bool + States whether pulse compression should be used + """ + + if waveform_mode not in ["CW", "BB"]: + raise ValueError("The input waveform_mode must be either 'CW' or 'BB'!") + + if encode_mode not in ["complex", "power"]: + raise ValueError("The input encode_mode must be either 'complex' or 'power'!") + + # BB has complex data only, but CW can have complex or power data + if (waveform_mode == "BB") and (encode_mode == "power"): + raise ValueError("encode_mode='power' not allowed when waveform_mode='BB'!") + + # make sure that we have BB and complex inputs, if pulse compression is selected + if pulse_compression and ((waveform_mode != "BB") or (encode_mode != "complex")): + raise ValueError( + "Pulse compression can only be used with " + "waveform_mode='BB' and encode_mode='complex'" + ) + + +def _retrieve_correct_beam_group_EK60( + echodata: EchoData, waveform_mode: str, encode_mode: str +) -> Optional[str]: + """ + Ensures that the provided ``waveform_mode`` and ``encode_mode`` are consistent + with the EK60-like data supplied by ``echodata``. Additionally, select the + appropriate beam group corresponding to this input. + + Parameters + ---------- + echodata: EchoData + An ``EchoData`` object holding the data + waveform_mode : {"CW", "BB"} + Type of transmit waveform + encode_mode : {"complex", "power"} + Type of encoded return echo data + + Returns + ------- + power_ed_group: str, optional + The ``EchoData`` beam group path containing the power data + """ + + # initialize power EchoData group value + power_ed_group = None + + # EK60-like sensors must have 'power' and 'CW' modes only + if waveform_mode != "CW": + raise RuntimeError("Incorrect waveform_mode input provided!") + if encode_mode != "power": + raise RuntimeError("Incorrect encode_mode input provided!") + + # ensure that no complex data exists (this should never be triggered) + if "backscatter_i" in echodata["Sonar/Beam_group1"].variables: + raise RuntimeError( + "Provided echodata object does not correspond to an EK60-like " + "sensor, but is labeled as data from an EK60-like sensor!" + ) + else: + power_ed_group = "Sonar/Beam_group1" + + return power_ed_group + + +def _retrieve_correct_beam_group_EK80( + echodata: EchoData, waveform_mode: str, encode_mode: str +) -> Tuple[Optional[str], Optional[str]]: + """ + Ensures that the provided ``waveform_mode`` and ``encode_mode`` are consistent + with the EK80-like data supplied by ``echodata``. Additionally, select the + appropriate beam group corresponding to this input. + + Parameters + ---------- + echodata: EchoData + An ``EchoData`` object holding the data + waveform_mode : {"CW", "BB"} + Type of transmit waveform + encode_mode : {"complex", "power"} + Type of encoded return echo data + + Returns + ------- + power_ed_group: str, optional + The ``EchoData`` beam group path containing the power data + complex_ed_group: str, optional + The ``EchoData`` beam group path containing the complex data + """ + + # initialize power and complex EchoData group values + power_ed_group = None + complex_ed_group = None + + if waveform_mode == "BB": + + # check BB waveform_mode, BB must always have complex data, can have 2 beam groups + # when echodata contains CW power and BB complex samples, and frequency_start + # variable in Beam_group1 + if waveform_mode == "BB" and "frequency_start" not in echodata["Sonar/Beam_group1"]: + raise RuntimeError("waveform_mode='BB', but broadband data not found!") + elif "backscatter_i" not in echodata["Sonar/Beam_group1"].variables: + raise RuntimeError("waveform_mode='BB', but complex data does not exist!") + elif echodata["Sonar/Beam_group2"] is not None: + power_ed_group = "Sonar/Beam_group2" + complex_ed_group = "Sonar/Beam_group1" + else: + complex_ed_group = "Sonar/Beam_group1" + + else: + + # CW can have complex or power data, so we just need to make sure that + # 1) complex samples always exist in Sonar/Beam_group1 + # 2) power samples are in Sonar/Beam_group1 if only one beam group exists + # 3) power samples are in Sonar/Beam_group2 if two beam groups exist + if echodata["Sonar/Beam_group2"] is None: + + if encode_mode == "power": + # power samples must be in Sonar/Beam_group1 (thus no complex samples) + if "backscatter_i" in echodata["Sonar/Beam_group1"].variables: + raise RuntimeError("Data provided does not correspond to encode_mode='power'!") + else: + power_ed_group = "Sonar/Beam_group1" + elif encode_mode == "complex": + # complex samples must be in Sonar/Beam_group1 + if "backscatter_i" not in echodata["Sonar/Beam_group1"].variables: + raise RuntimeError( + "Data provided does not correspond to encode_mode='complex'!" + ) + else: + complex_ed_group = "Sonar/Beam_group1" + + else: + + # complex should be in Sonar/Beam_group1 and power in Sonar/Beam_group2 + # the RuntimeErrors below should never be triggered + if "backscatter_i" not in echodata["Sonar/Beam_group1"].variables: + raise RuntimeError( + "Complex data does not exist in Sonar/Beam_group1, " + "input echodata object must have been incorrectly constructed!" + ) + elif "backscatter_r" not in echodata["Sonar/Beam_group2"].variables: + raise RuntimeError( + "Power data does not exist in Sonar/Beam_group2, " + "input echodata object must have been incorrectly constructed!" + ) + else: + complex_ed_group = "Sonar/Beam_group1" + power_ed_group = "Sonar/Beam_group2" + + return power_ed_group, complex_ed_group + + +def retrieve_correct_beam_group( + echodata: EchoData, waveform_mode: str, encode_mode: str, pulse_compression: bool +) -> str: + """ + A function to make sure that the user has provided the correct + ``waveform_mode`` and ``encode_mode`` inputs based off of the + supplied ``echodata`` object. Additionally, determine the + ``EchoData`` beam group corresponding to ``encode_mode``. + + Parameters + ---------- + echodata: EchoData + An ``EchoData`` object holding the data corresponding to the + waveform and encode modes + waveform_mode : {"CW", "BB"} + Type of transmit waveform + encode_mode : {"complex", "power"} + Type of encoded return echo data + pulse_compression: bool + States whether pulse compression should be used + + Returns + ------- + str + The ``EchoData`` beam group path corresponding to the ``encode_mode`` input + """ + + # checks input and logic of modes without referencing data + _check_input_args_combination(waveform_mode, encode_mode, pulse_compression) + + if echodata.sonar_model in ["EK60", "ES70"]: + + # initialize complex_data_location (needed only for EK60) + complex_ed_group = None + + # check modes against data for EK60 and get power EchoData group + power_ed_group = _retrieve_correct_beam_group_EK60(echodata, waveform_mode, encode_mode) + + elif echodata.sonar_model in ["EK80", "ES80", "EA640"]: + + # check modes against data for EK80 and get power/complex EchoData groups + power_ed_group, complex_ed_group = _retrieve_correct_beam_group_EK80( + echodata, waveform_mode, encode_mode + ) + + else: + # raise error for unknown or unaccounted for sonar model + raise RuntimeError("EchoData was produced by a non-Simrad or unknown Simrad echo sounder!") + + if encode_mode == "complex": + return complex_ed_group + else: + return power_ed_group diff --git a/echopype/tests/consolidate/test_consolidate.py b/echopype/tests/consolidate/test_consolidate.py index a7d9162a7e..41070a5746 100644 --- a/echopype/tests/consolidate/test_consolidate.py +++ b/echopype/tests/consolidate/test_consolidate.py @@ -1,10 +1,25 @@ +import pathlib + import pytest import numpy as np import pandas as pd import xarray as xr - +import scipy.io as io import echopype as ep +from typing import List +import tempfile +import os + +""" +For future reference: + +For ``test_add_splitbeam_angle`` the test data is in the following locations: +- the EK60 raw file is in `test_data/ek60/DY1801_EK60-D20180211-T164025.raw` and the +associated echoview split-beam data is in `test_data/ek60/splitbeam`. +- the EK80 raw file is in `test_data/ek80_bb_with_calibration/2018115-D20181213-T094600.raw` and +the associated echoview split-beam data is in `test_data/ek80_bb_with_calibration/splitbeam` +""" @pytest.fixture( @@ -181,3 +196,144 @@ def _check_var(ds_test): ds_sel = ep.consolidate.add_location(ds=ds, echodata=ed, nmea_sentence="GGA") _check_var(ds_sel) + + +def _create_array_list_from_echoview_mats(paths_to_echoview_mat: List[pathlib.Path]) -> List[np.ndarray]: + """ + Opens each mat file in ``paths_to_echoview_mat``, selects the first ``ping_time``, + and then stores the array in a list. + + Parameters + ---------- + paths_to_echoview_mat: list of pathlib.Path + A list of paths corresponding to mat files, where each mat file contains the + echoview generated angle alongship and athwartship data for a channel + + Returns + ------- + list of np.ndarray + A list of numpy arrays generated by choosing the appropriate data from the mat files. + This list will have the same length as ``paths_to_echoview_mat`` + """ + + list_of_mat_arrays = [] + for mat_file in paths_to_echoview_mat: + + # open mat file and grab appropriate data + list_of_mat_arrays.append(io.loadmat(file_name=mat_file)["P0"]["Data_values"][0][0]) + + return list_of_mat_arrays + + +@pytest.mark.parametrize( + ("sonar_model", "test_path_key", "raw_file_name", "paths_to_echoview_mat", + "waveform_mode", "encode_mode", "pulse_compression", "write_Sv_to_file"), + [ + ( + "EK60", "EK60", "DY1801_EK60-D20180211-T164025.raw", + [ + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T1.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T2.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T3.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T4.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T5.mat' + ], + "CW", "power", False, False + ), + ( + "EK60", "EK60", "DY1801_EK60-D20180211-T164025.raw", + [ + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T1.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T2.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T3.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T4.mat', + 'splitbeam/DY1801_EK60-D20180211-T164025_angles_T5.mat' + ], + "CW", "power", False, True + ), + ( + "EK80", "EK80_CAL", "2018115-D20181213-T094600.raw", + [ + 'splitbeam/2018115-D20181213-T094600_angles_T1.mat', + 'splitbeam/2018115-D20181213-T094600_angles_T4.mat', + 'splitbeam/2018115-D20181213-T094600_angles_T6.mat', + 'splitbeam/2018115-D20181213-T094600_angles_T5.mat' + ], + "CW", "complex", False, False + ), + ( + "EK80", "EK80_CAL", "2018115-D20181213-T094600.raw", + [ + 'splitbeam/2018115-D20181213-T094600_angles_T3_nopc.mat', + 'splitbeam/2018115-D20181213-T094600_angles_T2_nopc.mat', + ], + "BB", "complex", False, False, + ), + ], + ids=["ek60_CW_power", "ek60_CW_power_Sv_path", "ek80_CW_complex", "ek80_BB_complex_no_pulse"] +) +def test_add_splitbeam_angle(sonar_model, test_path_key, raw_file_name, test_path, + paths_to_echoview_mat, waveform_mode, encode_mode, + pulse_compression, write_Sv_to_file): + + # obtain the EchoData object with the data needed for the calculation + ed = ep.open_raw(test_path[test_path_key] / raw_file_name, sonar_model=sonar_model) + + # compute Sv as it is required for the split-beam angle calculation + ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode=waveform_mode, encode_mode=encode_mode) + + # initialize temporary directory object + temp_dir = None + + # allows us to test for the case when source_Sv is a path + if write_Sv_to_file: + + # create temporary directory for mask_file + temp_dir = tempfile.TemporaryDirectory() + + # write DataArray to temporary directory + zarr_path = os.path.join(temp_dir.name, "Sv_data.zarr") + ds_Sv.to_zarr(zarr_path) + + # assign input to a path + ds_Sv = zarr_path + + # add the split-beam angles to an empty Dataset + ds_Sv = ep.consolidate.add_splitbeam_angle(source_Sv=ds_Sv, echodata=ed, + waveform_mode=waveform_mode, + encode_mode=encode_mode, + pulse_compression=pulse_compression) + + # obtain corresponding echoview output + full_echoview_path = [test_path[test_path_key] / path for path in paths_to_echoview_mat] + echoview_arr_list = _create_array_list_from_echoview_mats(full_echoview_path) + + # compare echoview output against computed output for all channels + for chan_ind in range(len(echoview_arr_list)): + + # grabs the appropriate ds data to compare against + reduced_angle_alongship = ds_Sv.isel(channel=chan_ind, ping_time=0).angle_alongship.dropna("range_sample") + reduced_angle_athwartship = ds_Sv.isel(channel=chan_ind, ping_time=0).angle_athwartship.dropna("range_sample") + + # TODO: make "start" below a parameter in the input so that this is not ad-hoc but something known + # for some files the echoview data is shifted by one index, here we account for that + if reduced_angle_alongship.shape == (echoview_arr_list[chan_ind].shape[1], ): + start = 0 + else: + start = 1 + + # check the computed angle_alongship values against the echoview output + assert np.allclose(reduced_angle_alongship.values[start:], + echoview_arr_list[chan_ind][0, :], rtol=1e-1, atol=1e-2) + + # check the computed angle_alongship values against the echoview output + assert np.allclose(reduced_angle_athwartship.values[start:], + echoview_arr_list[chan_ind][1, :], rtol=1e-1, atol=1e-2) + + if temp_dir: + # remove the temporary directory, if it was created + temp_dir.cleanup() + + +# TODO: need a test for power/angle data, with mock EchoData object +# containing some channels with single-beam data and some channels with split-beam data \ No newline at end of file diff --git a/echopype/tests/echodata/test_echodata_simrad.py b/echopype/tests/echodata/test_echodata_simrad.py new file mode 100644 index 0000000000..e8b3184727 --- /dev/null +++ b/echopype/tests/echodata/test_echodata_simrad.py @@ -0,0 +1,69 @@ +""" +Tests functions contained within echodata/simrad.py +""" +import pytest +from echopype.echodata.simrad import retrieve_correct_beam_group, _check_input_args_combination + + +@pytest.mark.parametrize( + ("waveform_mode", "encode_mode", "pulse_compression"), + [ + pytest.param("CW", "comp_power", None, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since comp_power ' + 'is not an acceptable choice for encode_mode.')), + pytest.param("CB", None, None, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since CB is not an ' + 'acceptable choice for waveform_mode.')), + pytest.param("BB", "power", None, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since BB and power is ' + 'not an acceptable combination.')), + pytest.param("BB", "power", True, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since BB and complex ' + 'must be used if pulse_compression is True.')), + pytest.param("CW", "complex", True, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since BB and complex ' + 'must be used if pulse_compression is True.')), + pytest.param("CW", "power", True, + marks=pytest.mark.xfail(strict=True, + reason='This test should fail since BB and complex ' + 'must be used if pulse_compression is True.')), + ("CW", "complex", False), + ("CW", "power", False), + ("BB", "complex", False), + ("BB", "complex", True), + + ], + ids=["incorrect_encode_mode", "incorrect_waveform_mode", "BB_power_combo", + "BB_power_pc_True", "CW_complex_pc_True", "CW_power_pc_True", "CW_complex_pc_False", + "CW_power_pc_False", "BB_complex_pc_False", "BB_complex_pc_True"] +) +def test_check_input_args_combination(waveform_mode: str, encode_mode: str, + pulse_compression: bool): + """ + Ensures that ``check_input_args_combination`` functions correctly when + provided various combinations of the input parameters. + + Parameters + ---------- + waveform_mode: str + Type of transmit waveform + encode_mode: str + Type of encoded return echo data + pulse_compression: bool + States whether pulse compression should be used + """ + + _check_input_args_combination(waveform_mode, encode_mode, pulse_compression) + + +def test_retrieve_correct_beam_group(): + + # TODO: create this test once we are happy with the form of retrieve_correct_beam_group + + pytest.skip("We need to add tests for retrieve_correct_beam_group!") +