diff --git a/echopype/mask/__init__.py b/echopype/mask/__init__.py index c2233c112..5867b1656 100644 --- a/echopype/mask/__init__.py +++ b/echopype/mask/__init__.py @@ -1,4 +1,3 @@ from .api import apply_mask, frequency_differencing, get_seabed_mask, shoal_weill __all__ = ["frequency_differencing", "apply_mask", "get_seabed_mask", "shoal_weill"] - diff --git a/echopype/mask/api.py b/echopype/mask/api.py index 90ac0c6ca..26d2c2fc5 100644 --- a/echopype/mask/api.py +++ b/echopype/mask/api.py @@ -7,17 +7,10 @@ import xarray as xr from pandas import Index -from echopype.mask.seabed import ( - _get_seabed_mask_ariza, - _get_seabed_mask_blackwell, - _get_seabed_mask_blackwell_mod, - _get_seabed_mask_deltaSv, - _get_seabed_mask_experimental, - _get_seabed_mask_maxSv, -) - -from ..utils.io import validate_source_ds_da +from ..utils.io import get_dataset, validate_source_ds_da +from ..utils.misc import frequency_nominal_to_channel from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level +from . import seabed from .freq_diff import _check_freq_diff_source_Sv, _parse_freq_diff_eq from .shoal import _weill as shoal_weill @@ -534,72 +527,41 @@ def frequency_differencing( return da -def get_shoal_mask( - source_Sv: Union[xr.Dataset, str, pathlib.Path], - desired_channel: str, - mask_type: str = "will", - **kwargs, -): +def create_multichannel_mask(masks: [xr.Dataset], channels: [str]) -> xr.Dataset: """ - Wrapper function for (future) multiple shoal masking algorithms - (currently, only MOVIES-B (Will) is implemented) + Given a set of single-channel masks and a list of channels, + creates a multichannel mask - Args: - source_Sv: xr.Dataset or str or pathlib.Path - If a Dataset this value contains the Sv data to create a mask for, - else it specifies the path to a zarr or netcdf file containing - a Dataset. This input must correspond to a Dataset that has the - coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. - desired_channel: str specifying the channel to generate the mask on - mask_type: string specifying the algorithm to use - currently, 'weill' is the only one implemented + Parameters + ========== + masks(xr.Dataset): a list of single-channel masks + channels(str): a list of channel names Returns - ------- - mask: xr.DataArray - A DataArray containing the mask for the Sv data. Regions satisfying the thresholding - criteria are filled with ``True``, else the regions are filled with ``False``. - mask_: xr.DataArray - A DataArray containing the mask for areas in which shoals were searched. - Edge regions are filled with 'False', whereas the portion - in which shoals could be detected is 'True' - - - Raises - ------ - ValueError - If 'weill' is not given + mask: a multi-channel mask + ====== """ - assert mask_type in ["will"] - if mask_type == "will": - # Define a list of the keyword arguments your function can handle - valid_args = {"thr", "maxvgap", "maxhgap", "minvlen", "minhlen"} - # Filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask, mask_ = shoal_weill(source_Sv, desired_channel, **filtered_kwargs) - else: - raise ValueError("The provided mask type must be Will") - return_mask = xr.DataArray( - mask, - dims=("ping_time", "range_sample"), - coords={"ping_time": source_Sv.ping_time, "range_sample": source_Sv.range_sample}, - ) - return_mask_ = xr.DataArray( - mask_, - dims=("ping_time", "range_sample"), - coords={"ping_time": source_Sv.ping_time, "range_sample": source_Sv.range_sample}, + if len(masks) != len(channels): + raise ValueError("number of masks and of channels provided should be the same") + for i in range(0, len(masks)): + mask = masks[i] + if "channel" in mask.coords: + masks[i] = mask.isel(channel=0) + result = xr.concat( + masks, Index(channels, name="channel"), data_vars="all", coords="all", join="exact" ) - return return_mask, return_mask_ + return result def get_seabed_mask( source_Sv: Union[xr.Dataset, str, pathlib.Path], - desired_channel: str, - mask_type: str = "ariza", - **kwargs, + parameters: dict, + desired_channel: str = None, + desired_frequency: int = None, + method: str = "ariza", ) -> xr.DataArray: """ - Create a mask based on the identified signal attenuations of Sv values. + Create a mask based on the identified signal attenuations of Sv values at 38KHz. Parameters ---------- source_Sv: xr.Dataset or str or pathlib.Path @@ -608,9 +570,10 @@ def get_seabed_mask( a Dataset. This input must correspond to a Dataset that has the coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. desired_channel: str - channel to generate the mask for - mask_type: str with either "ariza", "experimental", "blackwell_mod", + desired_freuency: int - desired frequency, in case the channel isn't directly specified + method: str with either "ariza", "experimental", "blackwell_mod", "blackwell", "deltaSv", "maxSv" - based on the preferred method for signal attenuation mask generation + based on the preferred method for seabed mask generation Returns ------- xr.DataArray @@ -631,117 +594,83 @@ def get_seabed_mask( -------- """ - assert mask_type in [ - "ariza", - "experimental", - "blackwell_mod", - "blackwell", - "deltaSv", - "maxSv", - ], "mask_type must be either 'ariza', 'experimental', 'blackwell', 'maxSv', 'deltaSv'" - - channel_Sv = source_Sv.sel(channel=desired_channel) - Sv = channel_Sv["Sv"].values.T - r = source_Sv["echo_range"].values[0, 0] - - if mask_type == "ariza": - # Define a list of the keyword arguments your function can handle - valid_args = {"r0", "r1", "roff", "thr", "ec", "ek", "dc", "dk"} - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask = _get_seabed_mask_ariza(Sv, r, **filtered_kwargs) - elif mask_type == "experimental": - # Define a list of the keyword arguments your function can handle - valid_args = {"r0", "r1", "roff", "thr", "ns", "n_dil"} - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask = _get_seabed_mask_experimental(Sv, r, **filtered_kwargs) - elif mask_type == "blackwell": - # Define a list of the keyword arguments your function can handle - # valid_args = {"theta", "phi", "r0", "r1", "tSv", "ttheta", "tphi", "wtheta", "wphi"} - valid_args = {"r0", "r1", "tSv", "ttheta", "tphi", "wtheta", "wphi"} - theta = channel_Sv["angle_alongship"].values.T - phi = channel_Sv["angle_athwartship"].values.T - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask = _get_seabed_mask_blackwell(Sv, r, theta=theta, phi=phi, **filtered_kwargs) - elif mask_type == "blackwell_mod": - # Define a list of the keyword arguments your function can handle - valid_args = { - # "theta", - # "phi", - "r0", - "r1", - "tSv", - "ttheta", - "tphi", - "wtheta", - "wphi", - "rlog", - "tpi", - "freq", - "rank", - } - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - theta = channel_Sv["angle_alongship"].values.T - phi = channel_Sv["angle_athwartship"].values.T - mask = _get_seabed_mask_blackwell_mod(Sv, r, theta=theta, phi=phi, **filtered_kwargs) - elif mask_type == "deltaSv": - # Define a list of the keyword arguments your function can handle - valid_args = {"r0", "r1", "roff", "thr"} - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask = _get_seabed_mask_deltaSv(Sv, r, **filtered_kwargs) - elif mask_type == "maxSv": - # Define a list of the keyword arguments your function can handle - valid_args = {"r0", "r1", "roff", "thr"} - # Use dictionary comprehension to filter out any kwargs not in your list - filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} - mask = _get_seabed_mask_maxSv(Sv, r, **filtered_kwargs) - else: - raise ValueError( - "The provided mask_type must be 'ariza', " - + "'experimental', 'blackwell', 'maxSv' or 'deltaSv'!" - ) + source_Sv = get_dataset(source_Sv) + mask_map = { + "ariza": seabed._ariza, + "experimental": seabed._experimental, + "blackwell": seabed._blackwell, + "blackwell_mod": seabed._blackwell_mod, + "delta_Sv": seabed._deltaSv, + "max_Sv": seabed._maxSv, + } - mask = np.logical_not(mask.T) - return_mask = xr.DataArray( - mask, - dims=("ping_time", "range_sample"), - coords={"ping_time": source_Sv.ping_time, "range_sample": source_Sv.range_sample}, - ) - return return_mask + if method not in mask_map.keys(): + raise ValueError(f"Unsupported method: {method}") + if desired_channel is None: + if desired_frequency is None: + raise ValueError("Must specify either desired channel or desired frequency") + else: + desired_channel = frequency_nominal_to_channel(source_Sv, desired_frequency) + mask = mask_map[method](source_Sv, desired_channel, parameters) + return mask -def create_multichannel_mask(masks: [xr.Dataset], channels: [str]) -> xr.Dataset: - """ - Given a set of single-channel masks and a list of channels, - creates a multichannel mask +def get_seabed_mask_multichannel( + source_Sv: Union[xr.Dataset, str, pathlib.Path], + parameters: dict, + method: str = "ariza", +) -> xr.DataArray: + """ + Create a mask based on the identified signal attenuations of Sv values at 38KHz. Parameters - ========== - masks(xr.Dataset): a list of single-channel masks - channels(str): a list of channel names - + ---------- + source_Sv: xr.Dataset or str or pathlib.Path + If a Dataset this value contains the Sv data to create a mask for, + else it specifies the path to a zarr or netcdf file containing + a Dataset. This input must correspond to a Dataset that has the + coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. + method: str with either "ariza", "experimental", "blackwell_mod", + "blackwell", "deltaSv", "maxSv" + based on the preferred method for seabed mask generation Returns - mask: a multi-channel mask - ====== + ------- + xr.DataArray + A DataArray containing the mask for the Sv data. Regions satisfying the thresholding + criteria are filled with ``True``, else the regions are filled with ``False``. + + Raises + ------ + ValueError + If neither "ariza", "experimental", "blackwell_mod", + "blackwell", "deltaSv", "maxSv" are given + + Notes + ----- + + + Examples + -------- + """ - if len(masks) != len(channels): - raise ValueError("number of masks and of channels provided should be the same") - for i in range(0, len(masks)): - mask = masks[i] - if "channel" in mask.coords: - masks[i] = mask.isel(channel=0) - result = xr.concat( - masks, Index(channels, name="channel"), data_vars="all", coords="all", join="exact" - ) - return result + source_Sv = get_dataset(source_Sv) + channel_list = source_Sv["channel"].values + mask_list = [] + for channel in channel_list: + mask = get_seabed_mask( + source_Sv, + parameters=parameters, + desired_channel=channel, + method=method, + ) + mask_list.append(mask) + mask = create_multichannel_mask(mask_list, channel_list) + return mask -def get_shoal_mask_multichannel( +def get_shoal_mask( source_Sv: Union[xr.Dataset, str, pathlib.Path], + desired_channel: str, mask_type: str = "will", **kwargs, ): @@ -755,17 +684,17 @@ def get_shoal_mask_multichannel( else it specifies the path to a zarr or netcdf file containing a Dataset. This input must correspond to a Dataset that has the coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. + desired_channel: str specifying the channel to generate the mask on mask_type: string specifying the algorithm to use currently, 'weill' is the only one implemented Returns ------- mask: xr.DataArray - A DataArray containing the multichannel mask for the Sv data. - Regions satisfying the thresholding criteria are filled with ``True``, - else the regions are filled with ``False``. + A DataArray containing the mask for the Sv data. Regions satisfying the thresholding + criteria are filled with ``True``, else the regions are filled with ``False``. mask_: xr.DataArray - A DataArray containing the multichannel mask for areas in which shoals were searched. + A DataArray containing the mask for areas in which shoals were searched. Edge regions are filled with 'False', whereas the portion in which shoals could be detected is 'True' @@ -775,53 +704,70 @@ def get_shoal_mask_multichannel( ValueError If 'weill' is not given """ - channel_list = source_Sv["channel"].values - mask_list = [] - _mask_list = [] - for channel in channel_list: - mask, _mask = get_shoal_mask(source_Sv, channel, mask_type, **kwargs) - mask_list.append(mask) - _mask_list.append(_mask) - mask = create_multichannel_mask(mask_list, channel_list) - _mask = create_multichannel_mask(_mask_list, channel_list) - return mask, _mask + assert mask_type in ["will"] + if mask_type == "will": + # Define a list of the keyword arguments your function can handle + valid_args = {"thr", "maxvgap", "maxhgap", "minvlen", "minhlen"} + # Filter out any kwargs not in your list + filtered_kwargs = {k: v for k, v in kwargs.items() if k in valid_args} + mask, mask_ = shoal_weill(source_Sv, desired_channel, **filtered_kwargs) + else: + raise ValueError("The provided mask type must be Will") + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": source_Sv.ping_time, "range_sample": source_Sv.range_sample}, + ) + return_mask_ = xr.DataArray( + mask_, + dims=("ping_time", "range_sample"), + coords={"ping_time": source_Sv.ping_time, "range_sample": source_Sv.range_sample}, + ) + return return_mask, return_mask_ -def get_seabed_mask_multichannel( +def get_shoal_mask_multichannel( source_Sv: Union[xr.Dataset, str, pathlib.Path], - mask_type: str = "ariza", + mask_type: str = "will", **kwargs, -) -> xr.DataArray: +): """ - Create a mask based on the identified signal attenuations of Sv values. + Wrapper function for (future) multiple shoal masking algorithms + (currently, only MOVIES-B (Will) is implemented) + + Args: + source_Sv: xr.Dataset or str or pathlib.Path + If a Dataset this value contains the Sv data to create a mask for, + else it specifies the path to a zarr or netcdf file containing + a Dataset. This input must correspond to a Dataset that has the + coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. + mask_type: string specifying the algorithm to use + currently, 'weill' is the only one implemented - Parameters - ---------- - source_Sv: xr.Dataset or str or pathlib.Path - If a Dataset this value contains the Sv data to create a mask for, - else it specifies the path to a zarr or netcdf file containing - a Dataset. This input must correspond to a Dataset that has the - coordinate ``channel`` and variables ``frequency_nominal`` and ``Sv``. - mask_type: str with either "ariza", "experimental", "blackwell_mod", - "blackwell", "deltaSv", "maxSv" - based on the preferred method for signal attenuation mask generation Returns ------- - xr.DataSet - A DataSet containing the multichannel mask for the Sv data. - Regions satisfying the thresholding criteria are filled with ``True``, - else the regions are filled with ``False``. + mask: xr.DataArray + A DataArray containing the multichannel mask for the Sv data. + Regions satisfying the thresholding criteria are filled with ``True``, + else the regions are filled with ``False``. + mask_: xr.DataArray + A DataArray containing the multichannel mask for areas in which shoals were searched. + Edge regions are filled with 'False', whereas the portion + in which shoals could be detected is 'True' + Raises ------ ValueError - If neither "ariza", "experimental", "blackwell_mod", - "blackwell", "deltaSv", "maxSv" are given + If 'weill' is not given """ channel_list = source_Sv["channel"].values mask_list = [] + _mask_list = [] for channel in channel_list: - mask = get_seabed_mask(source_Sv, channel, mask_type, **kwargs) + mask, _mask = get_shoal_mask(source_Sv, channel, mask_type, **kwargs) mask_list.append(mask) + _mask_list.append(_mask) mask = create_multichannel_mask(mask_list, channel_list) - return mask + _mask = create_multichannel_mask(_mask_list, channel_list) + return mask, _mask diff --git a/echopype/mask/seabed.py b/echopype/mask/seabed.py index d940925e7..9a9b401cb 100644 --- a/echopype/mask/seabed.py +++ b/echopype/mask/seabed.py @@ -33,14 +33,62 @@ import numpy as np import scipy.ndimage as nd_img +import xarray as xr from scipy.signal import convolve2d from skimage.measure import label from skimage.morphology import dilation, erosion, remove_small_objects, square from ..utils.mask_transformation import lin, log - -def _get_seabed_mask_maxSv(Sv, r, r0=10, r1=1000, roff=0, thr=(-40, -60)): +MAX_SV_DEFAULT_PARAMS = {"r0": 10, "r1": 1000, "roff": 0, "thr": (-40, -60)} +DELTA_SV_DEFAULT_PARAMS = {"r0": 10, "r1": 1000, "roff": 0, "thr": 20} +BLACKWELL_DEFAULT_PARAMS = { + "theta": None, + "phi": None, + "r0": 10, + "r1": 1000, + "tSv": -75, + "ttheta": 702, + "tphi": 282, + "wtheta": 28, + "wphi": 52, +} +BLACKWELL_MOD_DEFAULT_PARAMS = { + "theta": None, + "phi": None, + "r0": 10, + "r1": 1000, + "tSv": -75, + "ttheta": 702, + "tphi": 282, + "wtheta": 28, + "wphi": 52, + "rlog": None, + "tpi": None, + "freq": None, + "rank": 50, +} +EXPERIMENTAL_DEFAULT_PARAMS = { + "r0": 10, + "r1": 1000, + "roff": 0, + "thr": (-30, -70), + "ns": 150, + "n_dil": 3, +} +ARIZA_DEFAULT_PARAMS = { + "r0": 10, + "r1": 1000, + "roff": 0, + "thr": -40, + "ec": 1, + "ek": (1, 3), + "dc": 10, + "dk": (3, 7), +} + + +def _maxSv(Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS): """ Initially detects the seabed as the ping sample with the strongest Sv value, as long as it exceeds a dB threshold. Then it searches up along the ping @@ -48,16 +96,34 @@ def _get_seabed_mask_maxSv(Sv, r, r0=10, r1=1000, roff=0, thr=(-40, -60)): seabed is set. Args: - Sv (float): 2D Sv array (dB). - r (float): 1D range array (m). - r0 (int): minimum range below which the search will be performed (m). - r1 (int): maximum range above which the search will be performed (m). - roff (int): seabed range offset (m). - thr (tuple): 2 integers with 1st and 2nd Sv threshold (dB). + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m). + r1 (int): maximum range above which the search will be performed (m). + roff (int): seabed range offset (m). + thr (tuple): 2 integers with 1st and 2nd Sv threshold (dB). Returns: - bool: 2D array with seabed mask. + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = ["r0", "r1", "roff", "thr"] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + roff = parameters["roff"] + thr = parameters["thr"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] # get offset and range indexes roff = np.nanargmin(abs(r - roff)) @@ -93,27 +159,51 @@ def _get_seabed_mask_maxSv(Sv, r, r0=10, r1=1000, roff=0, thr=(-40, -60)): i = 0 mask[i:, j] = True - return mask + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask -def _get_seabed_mask_deltaSv(Sv, r, r0=10, r1=1000, roff=0, thr=20): +def _deltaSv(Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS): """ Examines the difference in Sv over a 2-samples moving window along every ping, and returns the range of the first value that exceeded a user-defined dB threshold (likely, the seabed). Args: - Sv (float): 2D Sv array (dB). - r (float): 1D range array (m). - r0 (int): minimum range below which the search will be performed (m). - r1 (int): maximum range above which the search will be performed (m). - roff (int): seabed range offset (m). - thr (int): threshold value (dB). - start (int): ping index to start processing. + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m). + r1 (int): maximum range above which the search will be performed (m). + roff (int): seabed range offset (m). + thr (int): threshold value (dB). Returns: - bool: 2D array with seabed mask. + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = ["r0", "r1", "roff", "thr"] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + roff = parameters["roff"] + thr = parameters["thr"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] + # get offset as number of samples roff = np.nanargmin(abs(r - roff)) @@ -140,43 +230,58 @@ def _get_seabed_mask_deltaSv(Sv, r, r0=10, r1=1000, roff=0, thr=20): i = 0 mask[i:, j] = True - return mask - - -def _get_seabed_mask_blackwell( - Sv, - r, - theta=None, - phi=None, - r0=10, - r1=1000, - tSv=-75, - ttheta=702, # on raw EK60 data - tphi=282, # on raw EK60 data - wtheta=28, - wphi=52, -): + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask + + +def _blackwell(Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS): """ Detects and mask seabed using the split-beam angle and Sv, based in "Blackwell et al (2019), Aliased seabed detection in fisheries acoustic data". Complete article here: https://arxiv.org/abs/1904.10736 Args: - Sv (float): 2D numpy array with Sv data (dB) - theta (float): 2D numpy array with the along-ship angle (degrees) - phi (float): 2D numpy array with the athwart-ship angle (degrees) - r (float): 1D range array (m) - r0 (int): minimum range below which the search will be performed (m) - r1 (int): maximum range above which the search will be performed (m) - tSv (float): Sv threshold above which seabed is pre-selected (dB) - ttheta (int): Theta threshold above which seabed is pre-selected (dB) - tphi (int): Phi threshold above which seabed is pre-selected (dB) - wtheta (int): window's size for mean square operation in Theta field - wphi (int): window's size for mean square operation in Phi field + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m) + r1 (int): maximum range above which the search will be performed (m) + tSv (float): Sv threshold above which seabed is pre-selected (dB) + ttheta (int): Theta threshold above which seabed is pre-selected (dB) + tphi (int): Phi threshold above which seabed is pre-selected (dB) + wtheta (int): window's size for mean square operation in Theta field + wphi (int): window's size for mean square operation in Phi field Returns: - bool: 2D array with seabed mask + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = ["r0", "r1", "tSv", "ttheta", "tphi", "wtheta", "wphi"] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + tSv = parameters["tSv"] + ttheta = parameters["ttheta"] + tphi = parameters["tphi"] + wtheta = parameters["wtheta"] + wphi = parameters["wphi"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] + theta = channel_Sv["angle_alongship"].values.T + phi = channel_Sv["angle_athwartship"].values.T # apply reverse correction on theta & phi to match Blackwell's constants theta = theta * 22 * 128 / 180 phi = phi * 22 * 128 / 180 @@ -228,25 +333,17 @@ def _get_seabed_mask_blackwell( else: mask = np.zeros_like(Sv, dtype=bool) - return mask - - -def _get_seabed_mask_blackwell_mod( - Sv, - r, - theta=None, - phi=None, - r0=10, - r1=1000, - tSv=-75, - ttheta=702, # on raw EK60 data rather than physical angle - tphi=282, # on raw EK60 data rather than physical angle - wtheta=28, - wphi=52, - rlog=None, - tpi=None, - freq=None, - rank=50, + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask + + +def _blackwell_mod( + Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS ): """ Detects and mask seabed using the split-beam angle and Sv, based in @@ -259,25 +356,62 @@ def _get_seabed_mask_blackwell_mod( rank. Args: - Sv (float): 2D numpy array with Sv data (dB) - theta (float): 2D numpy array with the along-ship angle (degrees) - phi (float): 2D numpy array with the athwart-ship angle (degrees) - r (float): 1D range array (m) - r0 (int): minimum range below which the search will be performed (m) - r1 (int): maximum range above which the search will be performed (m) - tSv (float): Sv threshold above which seabed is pre-selected (dB) - ttheta (int): Theta threshold above which seabed is pre-selected (dB) - tphi (int): Phi threshold above which seabed is pre-selected (dB) - wtheta (int): window's size for mean square operation in Theta field - wphi (int): window's size for mean square operation in Phi field - rlog (float): Maximum logging range of the echosounder (m) - tpi (float): Transmit pulse interval, or ping rate (s) - freq (int): frequecy (kHz) - rank (int): Rank for percentile operation: [0, 100] + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m) + r1 (int): maximum range above which the search will be performed (m) + tSv (float): Sv threshold above which seabed is pre-selected (dB) + ttheta (int): Theta threshold above which seabed is pre-selected (dB) + tphi (int): Phi threshold above which seabed is pre-selected (dB) + wtheta (int): window's size for mean square operation in Theta field + wphi (int): window's size for mean square operation in Phi field + rlog (float): Maximum logging range of the echosounder (m) + tpi (float): Transmit pulse interval, or ping rate (s) + freq (int): frequecy (kHz) + rank (int): Rank for percentile operation: [0, 100] Returns: - bool: 2D array with seabed mask + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = [ + "r0", + "r1", + "tSv", + "ttheta", + "tphi", + "wtheta", + "wphi", + "rlog", + "tpi", + "freq", + "rank", + ] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + tSv = parameters["tSv"] + ttheta = parameters["ttheta"] + tphi = parameters["tphi"] + wtheta = parameters["wtheta"] + wphi = parameters["wphi"] + rlog = parameters["rlog"] + tpi = parameters["tpi"] + freq = parameters["freq"] + rank = parameters["rank"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] + theta = channel_Sv["angle_alongship"].values.T + phi = channel_Sv["angle_athwartship"].values.T # apply reverse correction on theta & phi to match Blackwell's constants theta = theta * 22 * 128 / 180 phi = phi * 22 * 128 / 180 @@ -351,7 +485,13 @@ def _get_seabed_mask_blackwell_mod( else: mask = np.zeros_like(Sv, dtype=bool) - return mask + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask def _aliased2seabed( @@ -418,7 +558,9 @@ def _seabed2aliased( return aliased -def _get_seabed_mask_experimental(Sv, r, r0=10, r1=1000, roff=0, thr=(-30, -70), ns=150, n_dil=3): +def _experimental( + Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS +): """ Mask Sv above a threshold to get a potential seabed mask. Then, the mask is dilated to fill seabed breaches, and small objects are removed to prevent @@ -427,18 +569,38 @@ def _get_seabed_mask_experimental(Sv, r, r0=10, r1=1000, roff=0, thr=(-30, -70), threshold, Finally, the mask is extended all the way down. Args: - Sv (float): 2D Sv array (dB). - r (float): 1D range array (m). - r0 (int): minimum range below which the search will be performed (m). - r1 (int): maximum range above which the search will be performed (m). - roff (int): seabed range offset (m). - thr (tuple): 2 integers with 1st and 2nd Sv threshold (dB). - ns (int): maximum number of samples for an object to be removed. - n_dil (int): number of dilations performed to the seabed mask. + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m). + r1 (int): maximum range above which the search will be performed (m). + roff (int): seabed range offset (m). + thr (tuple): 2 integers with 1st and 2nd Sv threshold (dB). + ns (int): maximum number of samples for an object to be removed. + n_dil (int): number of dilations performed to the seabed mask. Returns: - bool: 2D array with seabed mask. + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = ["r0", "r1", "roff", "thr", "ns", "n_dil"] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + roff = parameters["roff"] + thr = parameters["thr"] + ns = parameters["ns"] + n_dil = parameters["n_dil"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] # get indexes for range offset and range limits roff = np.nanargmin(abs(r - roff)) @@ -479,12 +641,16 @@ def _get_seabed_mask_experimental(Sv, r, r0=10, r1=1000, roff=0, thr=(-30, -70), # mask = cv2.dilate(np.uint8(mask), kernel, iterations = 2) # mask = np.array(mask, dtype = 'bool') - return mask + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask -def _get_seabed_mask_ariza( - Sv, r, r0=10, r1=1000, roff=0, thr=-40, ec=1, ek=(1, 3), dc=10, dk=(3, 7) -): +def _ariza(Sv_ds: xr.DataArray, desired_channel: str, parameters: dict = MAX_SV_DEFAULT_PARAMS): """ Mask Sv above a threshold to get potential seabed features. These features are eroded first to get rid of fake seabeds (spikes, schools, etc.) and @@ -494,22 +660,44 @@ def _get_seabed_mask_ariza( reconmended for non-supervised processing. Args: - Sv (float): 2D Sv array (dB). - r (float): 1D range array (m). - r0 (int): minimum range below which the search will be performed (m). - r1 (int): maximum range above which the search will be performed (m). - roff (int): seabed range offset (m). - thr (int): Sv threshold above which seabed might occur (dB). - ec (int): number of erosion cycles. - ek (int): 2-elements tuple with vertical and horizontal dimensions - of the erosion kernel. - dc (int): number of dilation cycles. - dk (int): 2-elements tuple with vertical and horizontal dimensions - of the dilation kernel. + Sv_ds (xr.DataArray): xr.DataArray with Sv data for multiple channels (dB) + desired_channel(str): Name of the desired frequency channel + parameters: parameter dict, should contain: + r0 (int): minimum range below which the search will be performed (m). + r1 (int): maximum range above which the search will be performed (m). + roff (int): seabed range offset (m). + thr (int): Sv threshold above which seabed might occur (dB). + ec (int): number of erosion cycles. + ek (int): 2-elements tuple with vertical and horizontal dimensions + of the erosion kernel. + dc (int): number of dilation cycles. + dk (int): 2-elements tuple with vertical and horizontal dimensions + of the dilation kernel. Returns: - bool: 2D array with seabed mask. + xr.DataArray: A DataArray containing the mask for the Sv data. + Regions satisfying the thresholding criteria are True, others are False """ + parameter_names = ["r0", "r1", "roff", "thr", "ec", "ek", "dc", "dk"] + if not all(name in parameters.keys() for name in parameter_names): + raise ValueError( + "Missing parameters - should be: " + + str(parameter_names) + + ", are: " + + str(parameters.keys()) + ) + r0 = parameters["r0"] + r1 = parameters["r1"] + roff = parameters["roff"] + thr = parameters["thr"] + ec = parameters["ec"] + ek = parameters["ek"] + dc = parameters["dc"] + dk = parameters["dk"] + + channel_Sv = Sv_ds.sel(channel=desired_channel) + Sv = channel_Sv["Sv"].values.T + r = Sv_ds["echo_range"].values[0, 0] # raise errors if wrong arguments if r0 > r1: @@ -562,4 +750,11 @@ def _get_seabed_mask_ariza( i = 0 mask[i:, j] = True - return mask + mask = np.logical_not(mask.T) + return_mask = xr.DataArray( + mask, + dims=("ping_time", "range_sample"), + coords={"ping_time": channel_Sv.ping_time, "range_sample": channel_Sv.range_sample}, + ) + return return_mask + \ No newline at end of file diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 0e22135f6..740448a20 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -130,4 +130,3 @@ def raw_dataset_jr179(setup_test_data_jr179): ed = _get_raw_dataset(setup_test_data_jr179) return ed - diff --git a/echopype/tests/mask/test_mask.py b/echopype/tests/mask/test_mask.py index 1152d5729..8c5552e39 100644 --- a/echopype/tests/mask/test_mask.py +++ b/echopype/tests/mask/test_mask.py @@ -1024,5 +1024,7 @@ def test_shoal_mask_all(sv_dataset_jr161): def test_seabed_mask_all(complete_dataset_jr179): source_Sv = complete_dataset_jr179 - ml = echopype.mask.api.get_seabed_mask_multichannel(source_Sv) + ml = ep.mask.api.get_seabed_mask_multichannel( + source_Sv, method="ariza", parameters=ep.mask.seabed.ARIZA_DEFAULT_PARAMS + ) assert np.all(ml["channel"] == source_Sv["channel"]) diff --git a/echopype/tests/mask/test_mask_seabed.py b/echopype/tests/mask/test_mask_seabed.py index f3e87f6ca..3aeb027b3 100644 --- a/echopype/tests/mask/test_mask_seabed.py +++ b/echopype/tests/mask/test_mask_seabed.py @@ -1,26 +1,35 @@ import numpy as np import pytest from echopype.mask.api import get_seabed_mask +from echopype.mask.seabed import ( + ARIZA_DEFAULT_PARAMS, + EXPERIMENTAL_DEFAULT_PARAMS, + BLACKWELL_DEFAULT_PARAMS, + BLACKWELL_MOD_DEFAULT_PARAMS, +) DESIRED_CHANNEL = "GPT 38 kHz 009072033fa5 1 ES38" @pytest.mark.parametrize( - "desired_channel,mask_type,expected_true_false_counts", + "desired_channel,method,parameters,expected_true_false_counts", [ - (DESIRED_CHANNEL, "ariza", (1430880, 736051)), - (DESIRED_CHANNEL, "experimental", (1514853, 652078)), - (DESIRED_CHANNEL, "blackwell", (1746551, 420380)), - (DESIRED_CHANNEL, "blackwell_mod", (1945202, 221729)), + (DESIRED_CHANNEL, "ariza", ARIZA_DEFAULT_PARAMS, (1430880, 736051)), + (DESIRED_CHANNEL, "experimental", EXPERIMENTAL_DEFAULT_PARAMS, (1514853, 652078)), + (DESIRED_CHANNEL, "blackwell", BLACKWELL_DEFAULT_PARAMS, (1746551, 420380)), + (DESIRED_CHANNEL, "blackwell_mod", BLACKWELL_MOD_DEFAULT_PARAMS, (1945202, 221729)), ], ) def test_mask_seabed( - complete_dataset_jr179, desired_channel, mask_type, expected_true_false_counts + complete_dataset_jr179, desired_channel, method, parameters, expected_true_false_counts ): source_Sv = complete_dataset_jr179 - mask = get_seabed_mask(source_Sv, desired_channel=desired_channel, mask_type=mask_type) + mask = get_seabed_mask( + source_Sv, desired_channel=desired_channel, method=method, parameters=parameters + ) count_true = np.count_nonzero(mask) count_false = mask.size - count_true true_false_counts = (count_true, count_false) assert true_false_counts == expected_true_false_counts + \ No newline at end of file diff --git a/requirements.txt b/requirements.txt index 42422871d..62a9315db 100644 --- a/requirements.txt +++ b/requirements.txt @@ -18,4 +18,4 @@ psutil>=5.9.1 geopy scikit-image flox>=0.7.2,<1.0.0 - +charset-normalizer<3.2