From 68d66b2ecd1ae62be43d31ba30fc9c4bad38eadc Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Tue, 20 Jan 2026 11:11:41 -0800 Subject: [PATCH 01/31] "Add regrid functionality based on match geometry" --- echopype/commongrid/utils.py | 207 +++++++++++++++++++++++++++++++++++ 1 file changed, 207 insertions(+) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index d831c39b1..eb30e6100 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -649,3 +649,210 @@ def assign_actual_range(ds_MVBS: xr.Dataset) -> xr.Dataset: round(float(ds_MVBS["Sv"].max().values), 2), ] return ds_MVBS.assign_attrs({"actual_range": actual_range}) + + +def get_valid_max_depth_ping(ds: xr.Dataset, channel_idx: int) -> int: + """ + Finds the integer indices of the last non-NaN value in a DataArray. + + Parameters + ---------- + ds : xr.Dataset + The xarray Dataset containing the variable of interest. + channel_idx : int + The index of the channel to access within the Dataset. + deepest_ping : int + The index of the ping (row) that contains the deepest valid (non-NaN) sample in the specified variable. + + Returns + ------- + deepest_ping: int + index of the "deepest" sample + """ + da = ds.isel(channel=channel_idx)["echo_range"].values + + deepest_ping = 0 + max_range_sample = -np.inf + + for i in range(da.shape[0]): + all_nan = np.isnan(da[i, :]) + row_indices = np.where(~all_nan)[0] + max_range_sample = max(max_range_sample, row_indices.max()) + if max_range_sample == row_indices.max(): + deepest_ping = i + return deepest_ping + + +def _weighted_mean_kernel(target_ranges, source_ranges, source_values): + """ + The Numba/Numpy kernel. + Resamples a single ping (1D array) from source geometry to target geometry. + + Parameters + ---------- + target_ranges: np.ndarray + 1D array of target bin centers (range values) to which the source data will be resampled. + source_ranges: np.ndarray + 1D array of source bin centers (range values) representing the original data geometry. + source_values: np.ndarray + 1D array of values corresponding to each source bin (e.g., power or intensity values). + + Returns + ------- + output: np.ndarray + 1D array of resampled values at the target bin centers. + """ + # 1. Quick Checks for Empty/Invalid Data + if np.all(np.isnan(target_ranges)) or np.all(np.isnan(source_ranges)): + return np.full_like(target_ranges, np.nan) + + valid_mask_source = ~np.isnan(source_ranges) & ~np.isnan(source_values) + valid_mask_target = ~np.isnan(target_ranges) + + if not np.any(valid_mask_source) or not np.any(valid_mask_target): + return np.full_like(target_ranges, np.nan) + + source_range = source_ranges[valid_mask_source] + source_value = source_values[valid_mask_source] + target_range_valid = target_ranges[valid_mask_target] + + # Define edges source + # Estimate bin edges from sample centers + if len(source_range) > 1: + source_mid = 0.5 * (source_range[:-1] + source_range[1:]) + source_edges = np.concatenate( + ( + [source_range[0] - (source_range[1] - source_range[0]) / 2], + source_mid, + [source_range[-1] + (source_range[-1] - source_range[-2]) / 2], + ) + ) + else: + source_edges = np.array([source_range[0] - 0.5, source_range[0] + 0.5]) + + # Define edges target + if len(target_range_valid) > 1: + target_mid = 0.5 * (target_range_valid[:-1] + target_range_valid[1:]) + target_edges = np.concatenate( + ( + [target_range_valid[0] - (target_range_valid[1] - target_range_valid[0]) / 2], + target_mid, + [target_range_valid[-1] + (target_range_valid[-1] - target_range_valid[-2]) / 2], + ) + ) + else: + target_edges = np.array([target_range_valid[0] - 0.5, target_range_valid[0] + 0.5]) + + # Weighted mean integration + # Calculate accumulated energy + source_thickness = np.diff(source_edges) + source_energies = source_value * source_thickness + source_cdf = np.concatenate(([0], np.cumsum(source_energies))) + + # Interpolate CDF onto Target Edges + target_cdf_interp = np.interp(target_edges, source_edges, source_cdf) + + # Differentiate to get Energy per Target Bin + target_energies = np.diff(target_cdf_interp) + target_thickness = np.diff(target_edges) + + with np.errstate(divide="ignore", invalid="ignore"): + target_mean_power = target_energies / target_thickness + + output = np.full_like(target_ranges, np.nan) + output[valid_mask_target] = target_mean_power + + return output + + +def regrid_all_channels(ds_Sv, target_channel_idx=0): + """ + Regrids all channels in the EchoData object to match the geometry + of the specified target channel index. + + Parameters + ---------- + ds_Sv : xr.Dataset + Input Dataset containing Sv data + target_channel_idx : int + The index of the channel to serve as the master grid. + + Returns + ------- + xr.Dataset + A new Dataset where all channels share the same `ping_time`, + `range_sample`, and `echo_range` as the target. + """ + + channels = ds_Sv.channel.values + + # Check bounds + if not (0 <= target_channel_idx < len(channels)): + raise IndexError( + f"Target index {target_channel_idx} out of range for {len(channels)} channels." + ) + + # Extract target ds + ds_target = ds_Sv.isel(channel=target_channel_idx).copy() + target_range_da = ds_target["echo_range"] + + # List to hold the aligned DataArrays + aligned_arrays = [] + + print(f"Master Grid: {channels[target_channel_idx]} (Index {target_channel_idx})") + + for i, channel in enumerate(channels): + + if i == target_channel_idx: + print(f" - Copying {channel} (Target)...") + aligned_arrays.append(ds_target["Sv"]) + continue + + print(f" - Matching {channel} (Index {i}) to Target...") + ds_source = ds_Sv.isel(channel=i) + + ds_source_aligned = ds_source.reindex(ping_time=ds_target.ping_time, method="nearest") + + # Linear domain for resampling + s_linear = 10 ** (ds_source_aligned["Sv"] / 10.0) + source_range_da = ds_source_aligned["echo_range"] + + # Apply weighted mean resapling as Ufunc + + result_linear = xr.apply_ufunc( + _weighted_mean_kernel, + target_range_da, + source_range_da, + s_linear, + input_core_dims=[["range_sample"], ["range_sample"], ["range_sample"]], + output_core_dims=[["range_sample"]], + vectorize=True, + dask="parallelized", + output_dtypes=[float], + ) + + # Convert back to log domain + result_linear = result_linear.where(result_linear > 0) + result_sv = 10 * np.log10(result_linear) + + result_sv.name = "Sv" + + result_sv = result_sv.assign_coords(channel=channel) + + aligned_arrays.append(result_sv) + + ds_combined = xr.concat(aligned_arrays, dim="channel") + + # Construct final Dataset based on original ds_Sv + deepest_ping_index = get_valid_max_depth_ping(ds_Sv, channel_idx=target_channel_idx) + valid_range_sample = np.argmax( + ds_Sv.isel(channel=target_channel_idx, ping_time=deepest_ping_index)["echo_range"].values + ) + + ds_final = ds_Sv.copy(deep=True) + + ds_final["echo_range"][:] = ds_Sv["echo_range"].isel(channel=target_channel_idx) + ds_final = ds_final.isel(range_sample=slice(0, valid_range_sample)) + ds_final["Sv"] = ds_combined.isel(range_sample=slice(0, valid_range_sample)) + + return ds_final From 1d1ff7643e3b840e0464d32a89924dc17ff9e548 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Tue, 20 Jan 2026 11:43:42 -0800 Subject: [PATCH 02/31] Fix format issue --- echopype/commongrid/utils.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index eb30e6100..bb8c4d6dc 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -662,7 +662,7 @@ def get_valid_max_depth_ping(ds: xr.Dataset, channel_idx: int) -> int: channel_idx : int The index of the channel to access within the Dataset. deepest_ping : int - The index of the ping (row) that contains the deepest valid (non-NaN) sample in the specified variable. + The index of the ping that contains the deepest valid sample in the specified variable. Returns ------- @@ -691,11 +691,11 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): Parameters ---------- target_ranges: np.ndarray - 1D array of target bin centers (range values) to which the source data will be resampled. + 1D array of target bin centers to which the source data will be resampled. source_ranges: np.ndarray - 1D array of source bin centers (range values) representing the original data geometry. + 1D array of source bin centers representing the original data geometry. source_values: np.ndarray - 1D array of values corresponding to each source bin (e.g., power or intensity values). + 1D array of values corresponding to each source bin. Returns ------- From af90b803910766daede4f0df3a15405f1010fde9 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Fri, 23 Jan 2026 12:59:41 -0800 Subject: [PATCH 03/31] Move regrid function, and add user specific grid --- echopype/commongrid/api.py | 103 ++++++++++++++++++++++++- echopype/commongrid/utils.py | 144 +++++++---------------------------- 2 files changed, 128 insertions(+), 119 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 54e3b71f9..c4a219773 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -21,6 +21,8 @@ compute_raw_MVBS, compute_raw_NASC, get_distance_from_latlon, + _weighted_mean_kernel, + get_valid_max_depth_ping, ) logger = logging.getLogger(__name__) @@ -441,5 +443,102 @@ def compute_NASC( return ds_NASC -def regrid(): - return 1 +def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) -> xr.Dataset: + """ + Regrids all channels in the EchoData object to match the geometry + of the specified target channel index. + + Parameters + ---------- + ds_Sv : xr.Dataset + Input Dataset containing Sv data + target_channel : string, optional + The label of the channel to serve as the master grid. + target_grid : xr.DataArray, optional + A data array containing specific 'echo_range' values specifying the target grid. + Data array must have dimension ('ping_time', 'range_sample'). + Returns + ------- + xr.Dataset + A new Dataset where all channels share the same `ping_time`, + `range_sample`, and `echo_range` as the target. + """ + if target_channel is not None and target_grid is not None: + raise ValueError("Provide only one of target_channel or target_grid, not both.") + + channels = ds_Sv.channel.values + ds_final = ds_Sv.copy(deep=True) + + # Check bounds + if target_channel and not target_channel in channels: + raise IndexError( + f"{target_channel} is not part of the channel names in : {channels}" + ) + + if target_grid is not None and target_grid.dims != ("ping_time", "range_sample"): + raise NotImplementedError("target_grid dimensions do not match expected dimensions.") + + # Target channel is given + if target_channel: + ds_target = ds_Sv.sel(channel=target_channel).copy() + target_range_da = ds_target["echo_range"] + + deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_channel=target_channel) + valid_range_sample = np.argmax( + ds_Sv["echo_range"].sel(channel=target_channel).isel(ping_time=deepest_ping_index).values + ) + ds_final["echo_range"][:] = ds_Sv["echo_range"].sel(channel=target_channel) + + # Target grid is given + else: + target_range_da = target_grid.copy() + + deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_grid=target_grid) + valid_range_sample = np.argmax( + target_grid.isel(ping_time=deepest_ping_index).values + ) + ds_final["echo_range"][:] = target_grid + + # List to hold the aligned DataArrays + aligned_arrays = [] + + for i, channel in enumerate(channels): + + ds_source = ds_Sv.sel(channel=channel) + + # Linear domain for resampling + source_linear = 10 ** (ds_source["Sv"] / 10.0) + source_range_da = ds_source["echo_range"] + + # Apply weighted mean resapling as Ufunc + + result_linear = xr.apply_ufunc( + _weighted_mean_kernel, + target_range_da, + source_range_da, + source_linear, + input_core_dims=[["range_sample"], ["range_sample"], ["range_sample"]], + output_core_dims=[["range_sample"]], + vectorize=True, + dask="parallelized", + output_dtypes=[np.float64], + ) + + # Convert back to log domain + result_linear = result_linear.where(result_linear > 0) + result_sv = 10 * np.log10(result_linear) + + result_sv.name = "Sv" + + result_sv = result_sv.assign_coords(channel=channel) + + aligned_arrays.append(result_sv) + + ds_combined = xr.concat(aligned_arrays, dim="channel") + + # Construct final Dataset based on original ds_Sv + + ds_final = ds_final.isel(range_sample=slice(0, valid_range_sample + 20)) + ds_final["Sv"] = ds_combined.isel(range_sample=slice(0, valid_range_sample + 20)) + + return ds_final diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index bb8c4d6dc..818248b2a 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -651,7 +651,7 @@ def assign_actual_range(ds_MVBS: xr.Dataset) -> xr.Dataset: return ds_MVBS.assign_attrs({"actual_range": actual_range}) -def get_valid_max_depth_ping(ds: xr.Dataset, channel_idx: int) -> int: +def get_valid_max_depth_ping(ds: xr.Dataset, target_channel: str = None, target_grid: xr.DataArray = None) -> int: """ Finds the integer indices of the last non-NaN value in a DataArray. @@ -659,8 +659,8 @@ def get_valid_max_depth_ping(ds: xr.Dataset, channel_idx: int) -> int: ---------- ds : xr.Dataset The xarray Dataset containing the variable of interest. - channel_idx : int - The index of the channel to access within the Dataset. + channel : str + The label used for the channel to access within the Dataset. deepest_ping : int The index of the ping that contains the deepest valid sample in the specified variable. @@ -669,7 +669,11 @@ def get_valid_max_depth_ping(ds: xr.Dataset, channel_idx: int) -> int: deepest_ping: int index of the "deepest" sample """ - da = ds.isel(channel=channel_idx)["echo_range"].values + da = xr.DataArray() + if target_channel is not None: + da = ds.sel(channel=target_channel)["echo_range"].values + if target_grid is not None: + da = target_grid.values deepest_ping = 0 max_range_sample = -np.inf @@ -702,25 +706,23 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): output: np.ndarray 1D array of resampled values at the target bin centers. """ - # 1. Quick Checks for Empty/Invalid Data - if np.all(np.isnan(target_ranges)) or np.all(np.isnan(source_ranges)): - return np.full_like(target_ranges, np.nan) - valid_mask_source = ~np.isnan(source_ranges) & ~np.isnan(source_values) - valid_mask_target = ~np.isnan(target_ranges) + valid_source_range_mask = ~np.isnan(source_ranges) + valid_target_range_mask = ~np.isnan(target_ranges) - if not np.any(valid_mask_source) or not np.any(valid_mask_target): + if not np.any(valid_source_range_mask) or not np.any(valid_target_range_mask): return np.full_like(target_ranges, np.nan) - source_range = source_ranges[valid_mask_source] - source_value = source_values[valid_mask_source] - target_range_valid = target_ranges[valid_mask_target] + source_range = source_ranges[valid_source_range_mask].astype(np.float64) + target_range_valid = target_ranges[valid_target_range_mask].astype(np.float64) + + source_value = source_values[valid_source_range_mask].astype(np.float64) # Define edges source # Estimate bin edges from sample centers if len(source_range) > 1: source_mid = 0.5 * (source_range[:-1] + source_range[1:]) - source_edges = np.concatenate( + source_edge = np.concatenate( ( [source_range[0] - (source_range[1] - source_range[0]) / 2], source_mid, @@ -728,7 +730,7 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): ) ) else: - source_edges = np.array([source_range[0] - 0.5, source_range[0] + 0.5]) + source_edge = np.array([source_range[0] - 0.5, source_range[0] + 0.5]) # Define edges target if len(target_range_valid) > 1: @@ -742,15 +744,16 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): ) else: target_edges = np.array([target_range_valid[0] - 0.5, target_range_valid[0] + 0.5]) + source_value = np.nan_to_num(source_value, nan=0.0) - # Weighted mean integration - # Calculate accumulated energy - source_thickness = np.diff(source_edges) - source_energies = source_value * source_thickness - source_cdf = np.concatenate(([0], np.cumsum(source_energies))) + # Interpolate CDF onto Target Edges + source_thickness = np.diff(source_edge) + + source_value = source_value * source_thickness + source_cdf = np.concatenate(([0], np.cumsum(source_value))) # Interpolate CDF onto Target Edges - target_cdf_interp = np.interp(target_edges, source_edges, source_cdf) + target_cdf_interp = np.interp(target_edges, source_edge, source_cdf) # Differentiate to get Energy per Target Bin target_energies = np.diff(target_cdf_interp) @@ -759,100 +762,7 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): with np.errstate(divide="ignore", invalid="ignore"): target_mean_power = target_energies / target_thickness - output = np.full_like(target_ranges, np.nan) - output[valid_mask_target] = target_mean_power - - return output - - -def regrid_all_channels(ds_Sv, target_channel_idx=0): - """ - Regrids all channels in the EchoData object to match the geometry - of the specified target channel index. - - Parameters - ---------- - ds_Sv : xr.Dataset - Input Dataset containing Sv data - target_channel_idx : int - The index of the channel to serve as the master grid. - - Returns - ------- - xr.Dataset - A new Dataset where all channels share the same `ping_time`, - `range_sample`, and `echo_range` as the target. - """ - - channels = ds_Sv.channel.values - - # Check bounds - if not (0 <= target_channel_idx < len(channels)): - raise IndexError( - f"Target index {target_channel_idx} out of range for {len(channels)} channels." - ) - - # Extract target ds - ds_target = ds_Sv.isel(channel=target_channel_idx).copy() - target_range_da = ds_target["echo_range"] - - # List to hold the aligned DataArrays - aligned_arrays = [] - - print(f"Master Grid: {channels[target_channel_idx]} (Index {target_channel_idx})") - - for i, channel in enumerate(channels): - - if i == target_channel_idx: - print(f" - Copying {channel} (Target)...") - aligned_arrays.append(ds_target["Sv"]) - continue - - print(f" - Matching {channel} (Index {i}) to Target...") - ds_source = ds_Sv.isel(channel=i) - - ds_source_aligned = ds_source.reindex(ping_time=ds_target.ping_time, method="nearest") - - # Linear domain for resampling - s_linear = 10 ** (ds_source_aligned["Sv"] / 10.0) - source_range_da = ds_source_aligned["echo_range"] - - # Apply weighted mean resapling as Ufunc - - result_linear = xr.apply_ufunc( - _weighted_mean_kernel, - target_range_da, - source_range_da, - s_linear, - input_core_dims=[["range_sample"], ["range_sample"], ["range_sample"]], - output_core_dims=[["range_sample"]], - vectorize=True, - dask="parallelized", - output_dtypes=[float], - ) - - # Convert back to log domain - result_linear = result_linear.where(result_linear > 0) - result_sv = 10 * np.log10(result_linear) - - result_sv.name = "Sv" - - result_sv = result_sv.assign_coords(channel=channel) - - aligned_arrays.append(result_sv) - - ds_combined = xr.concat(aligned_arrays, dim="channel") - - # Construct final Dataset based on original ds_Sv - deepest_ping_index = get_valid_max_depth_ping(ds_Sv, channel_idx=target_channel_idx) - valid_range_sample = np.argmax( - ds_Sv.isel(channel=target_channel_idx, ping_time=deepest_ping_index)["echo_range"].values - ) - - ds_final = ds_Sv.copy(deep=True) - - ds_final["echo_range"][:] = ds_Sv["echo_range"].isel(channel=target_channel_idx) - ds_final = ds_final.isel(range_sample=slice(0, valid_range_sample)) - ds_final["Sv"] = ds_combined.isel(range_sample=slice(0, valid_range_sample)) + output = np.full_like(target_ranges, np.nan, dtype=np.float64) + output[valid_target_range_mask] = target_mean_power - return ds_final + return output \ No newline at end of file From adcb2e8badbba0255cd6f1ee2285c06c5822f8e0 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 15:12:26 -0800 Subject: [PATCH 04/31] Add test, and have higher precision in Sv --- echopype/commongrid/utils.py | 75 ++++- .../tests/commongrid/test_commongrid_api.py | 283 ++++++++++++++++-- 2 files changed, 325 insertions(+), 33 deletions(-) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index 818248b2a..cf95538d7 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -687,10 +687,65 @@ def get_valid_max_depth_ping(ds: xr.Dataset, target_channel: str = None, target_ return deepest_ping +def _exact_integration_kernel(target_edges, source_edge, source_value): + """ + Helper Kernel: Performs exact interval intersection to determine energy distribution. + Replaces the standard CDF interpolation method to avoid floating point cancellation errors. + """ + n_target = len(target_edges) - 1 + n_source = len(source_edge) - 1 + + # Result array + output_values = np.empty(n_target, dtype=np.float64) + + # Sliding window pointer for source bins (O(N+M) complexity) + s_idx = 0 + + for t_i in range(n_target): + t_left = target_edges[t_i] + t_right = target_edges[t_i + 1] + t_width = t_right - t_left + + if t_width <= 0: + output_values[t_i] = np.nan + continue + + total_energy = 0.0 + + # 1. Fast-forward source pointer to the start of overlap + # Stop when the NEXT source bin starts after the current target starts + while s_idx < n_source and source_edge[s_idx + 1] <= t_left: + s_idx += 1 + + # 2. Integrate overlapping source bins + # Use a temp pointer so we don't lose the place for the next target bin + temp_s = s_idx + while temp_s < n_source and source_edge[temp_s] < t_right: + s_left = source_edge[temp_s] + s_right = source_edge[temp_s + 1] + + # Calculate geometric overlap (Intersection) + overlap_start = max(t_left, s_left) + overlap_end = min(t_right, s_right) + overlap_len = overlap_end - overlap_start + + if overlap_len > 0: + val = source_value[temp_s] + # Accumulate energy (Density * Overlap Width) + total_energy += val * overlap_len + + temp_s += 1 + + output_values[t_i] = total_energy / t_width + + return output_values + def _weighted_mean_kernel(target_ranges, source_ranges, source_values): """ The Numba/Numpy kernel. Resamples a single ping (1D array) from source geometry to target geometry. + + Updated: Uses exact interval integration for high precision. Parameters ---------- @@ -744,23 +799,15 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): ) else: target_edges = np.array([target_range_valid[0] - 0.5, target_range_valid[0] + 0.5]) + source_value = np.nan_to_num(source_value, nan=0.0) - # Interpolate CDF onto Target Edges - source_thickness = np.diff(source_edge) - - source_value = source_value * source_thickness - source_cdf = np.concatenate(([0], np.cumsum(source_value))) - - # Interpolate CDF onto Target Edges - target_cdf_interp = np.interp(target_edges, source_edge, source_cdf) - - # Differentiate to get Energy per Target Bin - target_energies = np.diff(target_cdf_interp) - target_thickness = np.diff(target_edges) + # --- CHANGED: High Precision Integration --- + # Replaces the CDF/Interp/Diff method with exact geometric overlap. + + target_mean_power = _exact_integration_kernel(target_edges, source_edge, source_value) - with np.errstate(divide="ignore", invalid="ignore"): - target_mean_power = target_energies / target_thickness + # ------------------------------------------- output = np.full_like(target_ranges, np.nan, dtype=np.float64) output[valid_target_range_mask] = target_mean_power diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 20dc6c36c..709daaca5 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -10,7 +10,8 @@ _parse_x_bin, _groupby_x_along_channels, get_distance_from_latlon, - compute_raw_NASC + compute_raw_NASC, + _weighted_mean_kernel, ) from echopype.tests.commongrid.conftest import get_NASC_echoview @@ -46,9 +47,7 @@ def test__parse_x_bin(x_bin, x_label, expected_result): @pytest.mark.unit -@pytest.mark.parametrize( - ["range_var", "lat_lon"], [("depth", False), ("echo_range", False)] -) +@pytest.mark.parametrize(["range_var", "lat_lon"], [("depth", False), ("echo_range", False)]) def test__groupby_x_along_channels(request, range_var, lat_lon): """Testing the underlying function of compute_MVBS and compute_NASC""" range_bin = 20 @@ -75,7 +74,7 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): .indexes["ping_time"] ) ping_interval = d_index.union([d_index[-1] + pd.Timedelta(ping_time_bin)]) - + sv_mean = _groupby_x_along_channels( ds_Sv, range_interval, @@ -83,12 +82,250 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): x_var="ping_time", range_var=range_var, method=method, - **flox_kwargs + **flox_kwargs, ) # Check that the range_var is in the dimension assert f"{range_var}_bins" in sv_mean.dims +#Regrid Tests +@pytest.mark.unit +@pytest.mark.parametrize( + ["target_ranges", "source_ranges", "source_values", "expected_linear"], + [ + # Simple Downsample (2 Source -> 1 Target) + # Expected: Average of 1.0 and 1.0 is 1.0 + (np.array([15.0]), np.array([10.0, 20.0]), np.array([0.0, 0.0]), np.array([1.0])), + # NaN Handling (Data Gap) + # Logic: NaN data is treated as Linear 0.0, but geometry is preserved. + # Expected: = 0.666 + ( + np.array([20.0]), + np.array([10.0, 20.0, 30.0]), + np.array([0.0, np.nan, 0.0]), + np.array([2 / 3]), + ), + # Case 4: Precision Check + # Ensure float64 and scaling handles tiny numbers without underflow + (np.array([10.0]), np.array([10.0]), np.array([-120.0]), np.array([10 ** (-12.0)])), + ], +) +def test__weighted_mean_kernel(target_ranges, source_ranges, source_values, expected_linear): + """Testing the underlying mathematical kernel of match_geometry""" + + # Run the kernel (which returns linear values) + result = _weighted_mean_kernel(target_ranges, source_ranges, source_values) + + # Assert + np.testing.assert_allclose( + result, expected_linear, rtol=1e-5, atol=1e-20, err_msg="Kernel calculation mismatch" + ) + +@pytest.mark.integration +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ("irregular"), + ], +) +def test_regrid_with_channel(request, er_type): + """Testing the regrid_all_channels_robust wrapper function""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + ds_regridded = ep.commongrid.regrid(ds_Sv, target_channel = "WBT 987763-15 ES38-7_ES") + + def calculate_total_energy(ds, channel): + """ + Calculates the total integrated energy (Linear Sv * Thickness) + per ping for a specific channel. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing 'Sv' and 'echo_range'. + channel_idx : int + Index of the channel to calculate. + + Returns + ------- + np.ndarray + 1D array of Total Energy per ping. + """ + # Select Data + ds_ch = ds.sel(channel=channel) + sv = ds_ch['Sv'].values + r = ds_ch['echo_range'].values # Shape: (ping_time, range_sample) + + # Calculate Thickness (Delta R) per ping. + n_pings = sv.shape[0] + total_energy = np.zeros(n_pings) + + # 2. Iterate per ping to handle ragged arrays (variable valid ranges) + for i in range(n_pings): + # Extract row + r_row = r[i, :] + sv_row = sv[i, :] + + # Mask + valid_r = ~np.isnan(r_row) + + if not np.any(valid_r): + total_energy[i] = 0.0 + continue + + # Get valid geometry + r_valid = r_row[valid_r] + sv_valid = sv_row[valid_r] + + # Calculate Thickness using Midpoints + if len(r_valid) > 1: + midpoints = 0.5 * (r_valid[:-1] + r_valid[1:]) + # Extrapolate edges for first and last bin + d_start = r_valid[1] - r_valid[0] + d_end = r_valid[-1] - r_valid[-2] + + edges = np.concatenate([ + [r_valid[0] - d_start/2], + midpoints, + [r_valid[-1] + d_end/2] + ]) + + thickness = np.diff(edges) + else: + # Single sample fallback (assume unit thickness or 1.0) + thickness = np.array([1.0]) + + # 4. Convert dB to Linear + linear_sv = 10 ** (sv_valid / 10.0) + + # Treat NaN Sv as 0.0 (Empty water) + linear_sv = np.nan_to_num(linear_sv, nan=0.0) + + # 5. Integrate + total_energy[i] = np.sum(linear_sv * thickness) + + return total_energy + + # Calculate total energy for original and regridded datasets + # Calculate total energy for all channels + channels = ds_Sv["channel"].values + total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] + total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] + + np.testing.assert_allclose( + total_energy_original, + total_energy_regridded, + atol=1e-5, + rtol=1e-5, + err_msg="Total energy was not conserved during regridding!" + ) + +@pytest.mark.integration +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ("irregular"), + ], +) +def test_regrid_with_grid(request, er_type): + """Testing the regrid_all_channels_robust wrapper function""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + ds_regridded = ep.commongrid.regrid(ds_Sv, target_grid = ds_Sv["echo_range"].sel(channel = "WBT 987763-15 ES38-7_ES")) + + def calculate_total_energy(ds, channel): + """ + Calculates the total integrated energy (Linear Sv * Thickness) + per ping for a specific channel. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing 'Sv' and 'echo_range'. + channel_idx : int + Index of the channel to calculate. + + Returns + ------- + np.ndarray + 1D array of Total Energy per ping. + """ + # Select Data + ds_ch = ds.sel(channel=channel) + sv = ds_ch['Sv'].values + r = ds_ch['echo_range'].values # Shape: (ping_time, range_sample) + + # Calculate Thickness (Delta R) per ping. + n_pings = sv.shape[0] + total_energy = np.zeros(n_pings) + + # 2. Iterate per ping to handle ragged arrays (variable valid ranges) + for i in range(n_pings): + # Extract row + r_row = r[i, :] + sv_row = sv[i, :] + + # Mask + valid_r = ~np.isnan(r_row) + + if not np.any(valid_r): + total_energy[i] = 0.0 + continue + + # Get valid geometry + r_valid = r_row[valid_r] + sv_valid = sv_row[valid_r] + + # Calculate Thickness using Midpoints + if len(r_valid) > 1: + midpoints = 0.5 * (r_valid[:-1] + r_valid[1:]) + # Extrapolate edges for first and last bin + d_start = r_valid[1] - r_valid[0] + d_end = r_valid[-1] - r_valid[-2] + + edges = np.concatenate([ + [r_valid[0] - d_start/2], + midpoints, + [r_valid[-1] + d_end/2] + ]) + + thickness = np.diff(edges) + else: + # Single sample fallback (assume unit thickness or 1.0) + thickness = np.array([1.0]) + + # 4. Convert dB to Linear + linear_sv = 10 ** (sv_valid / 10.0) + + # Treat NaN Sv as 0.0 (Empty water) + linear_sv = np.nan_to_num(linear_sv, nan=0.0) + + # 5. Integrate + total_energy[i] = np.sum(linear_sv * thickness) + + return total_energy + + # Calculate total energy for original and regridded datasets + # Calculate total energy for all channels + channels = ds_Sv["channel"].values + total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] + total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] + + np.testing.assert_allclose( + total_energy_original, + total_energy_regridded, + atol=1e-5, + rtol=1e-5, + err_msg="Total energy was not conserved during regridding!" + ) + + # NASC Tests @pytest.mark.integration @@ -256,6 +493,7 @@ def test_compute_MVBS_invalid_range_var(ds_Sv_echo_range_regular, range_var): else: pass + @pytest.mark.integration def test_compute_MVBS_w_dim_swapped(request): """ @@ -269,10 +507,14 @@ def test_compute_MVBS_w_dim_swapped(request): # Test to see if ping_time was resampled correctly expected_ping_time = ( - ds_Sv["ping_time"].resample(ping_time=ping_time_bin, skipna=True).asfreq().indexes["ping_time"] + ds_Sv["ping_time"] + .resample(ping_time=ping_time_bin, skipna=True) + .asfreq() + .indexes["ping_time"] ) assert np.array_equal(ds_MVBS.ping_time.data, expected_ping_time.values) + @pytest.mark.integration def test_compute_MVBS(test_data_samples): """ @@ -465,9 +707,10 @@ def test_compute_NASC_values(request, er_type): ds_NASC.NASC.values, expected_nasc.values, atol=1e-10, rtol=1e-10, equal_nan=True ) + @pytest.mark.integration @pytest.mark.parametrize( - ("operation","skipna", "range_var"), + ("operation", "skipna", "range_var"), [ ("MVBS", True, "depth"), ("MVBS", False, "depth"), @@ -490,7 +733,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( # Get fixture for irregular Sv ds_Sv = request.getfixturevalue("mock_Sv_dataset_irregular") # Already has 2 channels and 20 range samples, so subset for only ping time - subset_ds_Sv = ds_Sv.isel(ping_time=slice(0,2)) + subset_ds_Sv = ds_Sv.isel(ping_time=slice(0, 2)) # Compute MVBS / Compute NASC if operation == "MVBS": @@ -499,10 +742,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( ep.utils.log.verbose(override=False) da = ep.commongrid.compute_MVBS( - subset_ds_Sv, - range_var=range_var, - range_bin="2m", - skipna=skipna + subset_ds_Sv, range_var=range_var, range_bin="2m", skipna=skipna )["Sv"] if range_var == "echo_range": @@ -512,7 +752,9 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( "```Sv``` values that have corresponding NaN coordinate values. Consider handling " "these values before calling your intended commongrid function." ) - expected_warning = f"The ```echo_range``` coordinate array contain NaNs. {aggregation_msg}" + expected_warning = ( + f"The ```echo_range``` coordinate array contain NaNs. {aggregation_msg}" + ) assert any(expected_warning in record.message for record in caplog.records) # Turn off logger verbosity @@ -532,7 +774,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( # have any `NaNs` that are aggregated into them. expected_values = [ [[False, False, False, False, False]], - [[False, False, False, False, False]] + [[False, False, False, False, False]], ] assert np.array_equal(da_nan_mask, np.array(expected_values)) else: @@ -542,13 +784,13 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( if skipna: expected_values = [ [[True, False, False, False, False, False]], - [[True, False, False, False, False, False]] + [[True, False, False, False, False, False]], ] assert np.array_equal(da_nan_mask, np.array(expected_values)) else: expected_values = [ [[True, True, True, False, False, True]], - [[True, False, False, True, True, True]] + [[True, False, False, True, True, True]], ] assert np.array_equal(da_nan_mask, np.array(expected_values)) @@ -570,7 +812,7 @@ def test_assign_actual_range(request): [ round(float(ds_MVBS["Sv"].min().values), 2), round(float(ds_MVBS["Sv"].max().values), 2), - ] + ], ) @@ -600,5 +842,8 @@ def test_compute_reindex_non_NaN_not_map_reduce(request): # Compute MVBS with invalid parameters for method in ["blockwise", "cohorts"]: for reindex in [True, False]: - with pytest.raises(ValueError, match=f"Passing in reindex={reindex} is only allowed when method='map_reduce'."): + with pytest.raises( + ValueError, + match=f"Passing in reindex={reindex} is only allowed when method='map_reduce'.", + ): ep.commongrid.compute_MVBS(ds_Sv, method=method, reindex=reindex) From bfe0cf0cfc4d2611d3b73390eaf8fcaec7af72c8 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 15:26:39 -0800 Subject: [PATCH 05/31] Refactor names --- echopype/commongrid/utils.py | 49 +++++------ .../tests/commongrid/test_commongrid_api.py | 85 +++++++------------ 2 files changed, 52 insertions(+), 82 deletions(-) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index cf95538d7..9d664372e 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -698,45 +698,40 @@ def _exact_integration_kernel(target_edges, source_edge, source_value): # Result array output_values = np.empty(n_target, dtype=np.float64) - # Sliding window pointer for source bins (O(N+M) complexity) - s_idx = 0 + source_idx = 0 - for t_i in range(n_target): - t_left = target_edges[t_i] - t_right = target_edges[t_i + 1] - t_width = t_right - t_left + for target_i in range(n_target): + target_left = target_edges[target_i] + target_right = target_edges[target_i + 1] + target_width = target_right - target_left - if t_width <= 0: - output_values[t_i] = np.nan + if target_width <= 0: + output_values[target_i] = np.nan continue total_energy = 0.0 - # 1. Fast-forward source pointer to the start of overlap - # Stop when the NEXT source bin starts after the current target starts - while s_idx < n_source and source_edge[s_idx + 1] <= t_left: - s_idx += 1 + # Fast-forward source pointer to the start of overlap + while source_idx < n_source and source_edge[source_idx + 1] <= target_left: + source_idx += 1 - # 2. Integrate overlapping source bins - # Use a temp pointer so we don't lose the place for the next target bin - temp_s = s_idx - while temp_s < n_source and source_edge[temp_s] < t_right: - s_left = source_edge[temp_s] - s_right = source_edge[temp_s + 1] + # Integrate overlapping source bins + temp_source = source_idx + while temp_source < n_source and source_edge[temp_source] < target_right: + source_left = source_edge[temp_source] + source_right = source_edge[temp_source + 1] - # Calculate geometric overlap (Intersection) - overlap_start = max(t_left, s_left) - overlap_end = min(t_right, s_right) + overlap_start = max(target_left, source_left) + overlap_end = min(target_right, source_right) overlap_len = overlap_end - overlap_start if overlap_len > 0: - val = source_value[temp_s] - # Accumulate energy (Density * Overlap Width) + val = source_value[temp_source] total_energy += val * overlap_len - temp_s += 1 + temp_source += 1 - output_values[t_i] = total_energy / t_width + output_values[target_i] = total_energy / target_width return output_values @@ -801,13 +796,9 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): target_edges = np.array([target_range_valid[0] - 0.5, target_range_valid[0] + 0.5]) source_value = np.nan_to_num(source_value, nan=0.0) - - # --- CHANGED: High Precision Integration --- - # Replaces the CDF/Interp/Diff method with exact geometric overlap. target_mean_power = _exact_integration_kernel(target_edges, source_edge, source_value) - # ------------------------------------------- output = np.full_like(target_ranges, np.nan, dtype=np.float64) output[valid_target_range_mask] = target_mean_power diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 709daaca5..8c4beab1b 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -94,18 +94,14 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): ["target_ranges", "source_ranges", "source_values", "expected_linear"], [ # Simple Downsample (2 Source -> 1 Target) - # Expected: Average of 1.0 and 1.0 is 1.0 (np.array([15.0]), np.array([10.0, 20.0]), np.array([0.0, 0.0]), np.array([1.0])), # NaN Handling (Data Gap) - # Logic: NaN data is treated as Linear 0.0, but geometry is preserved. - # Expected: = 0.666 ( np.array([20.0]), np.array([10.0, 20.0, 30.0]), np.array([0.0, np.nan, 0.0]), np.array([2 / 3]), ), - # Case 4: Precision Check # Ensure float64 and scaling handles tiny numbers without underflow (np.array([10.0]), np.array([10.0]), np.array([-120.0]), np.array([10 ** (-12.0)])), ], @@ -113,10 +109,8 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): def test__weighted_mean_kernel(target_ranges, source_ranges, source_values, expected_linear): """Testing the underlying mathematical kernel of match_geometry""" - # Run the kernel (which returns linear values) result = _weighted_mean_kernel(target_ranges, source_ranges, source_values) - # Assert np.testing.assert_allclose( result, expected_linear, rtol=1e-5, atol=1e-20, err_msg="Kernel calculation mismatch" ) @@ -155,62 +149,55 @@ def calculate_total_energy(ds, channel): 1D array of Total Energy per ping. """ # Select Data - ds_ch = ds.sel(channel=channel) - sv = ds_ch['Sv'].values - r = ds_ch['echo_range'].values # Shape: (ping_time, range_sample) + ds_channel = ds.sel(channel=channel) + sv = ds_channel['Sv'].values + echo_range = ds_channel['echo_range'].values - # Calculate Thickness (Delta R) per ping. n_pings = sv.shape[0] total_energy = np.zeros(n_pings) - # 2. Iterate per ping to handle ragged arrays (variable valid ranges) - for i in range(n_pings): + for i in echo_range(n_pings): # Extract row - r_row = r[i, :] + range_row = echo_range[i, :] sv_row = sv[i, :] # Mask - valid_r = ~np.isnan(r_row) + valid_range = ~np.isnan(range_row) - if not np.any(valid_r): + if not np.any(valid_range): total_energy[i] = 0.0 continue # Get valid geometry - r_valid = r_row[valid_r] - sv_valid = sv_row[valid_r] + range_valid = range_row[valid_range] + sv_valid = sv_row[valid_range] # Calculate Thickness using Midpoints - if len(r_valid) > 1: - midpoints = 0.5 * (r_valid[:-1] + r_valid[1:]) - # Extrapolate edges for first and last bin - d_start = r_valid[1] - r_valid[0] - d_end = r_valid[-1] - r_valid[-2] + if len(range_valid) > 1: + midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) + + d_start = range_valid[1] - range_valid[0] + d_end = range_valid[-1] - range_valid[-2] edges = np.concatenate([ - [r_valid[0] - d_start/2], + [range_valid[0] - d_start/2], midpoints, - [r_valid[-1] + d_end/2] + [range_valid[-1] + d_end/2] ]) thickness = np.diff(edges) else: - # Single sample fallback (assume unit thickness or 1.0) thickness = np.array([1.0]) - # 4. Convert dB to Linear linear_sv = 10 ** (sv_valid / 10.0) - # Treat NaN Sv as 0.0 (Empty water) linear_sv = np.nan_to_num(linear_sv, nan=0.0) - # 5. Integrate total_energy[i] = np.sum(linear_sv * thickness) return total_energy - # Calculate total energy for original and regridded datasets - # Calculate total energy for all channels + channels = ds_Sv["channel"].values total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] @@ -257,62 +244,54 @@ def calculate_total_energy(ds, channel): 1D array of Total Energy per ping. """ # Select Data - ds_ch = ds.sel(channel=channel) - sv = ds_ch['Sv'].values - r = ds_ch['echo_range'].values # Shape: (ping_time, range_sample) + ds_channel = ds.sel(channel=channel) + sv = ds_channel['Sv'].values + range = ds_channel['echo_range'].values - # Calculate Thickness (Delta R) per ping. n_pings = sv.shape[0] total_energy = np.zeros(n_pings) - # 2. Iterate per ping to handle ragged arrays (variable valid ranges) for i in range(n_pings): # Extract row - r_row = r[i, :] + range_row = range[i, :] sv_row = sv[i, :] # Mask - valid_r = ~np.isnan(r_row) + valid_range = ~np.isnan(range_row) - if not np.any(valid_r): + if not np.any(valid_range): total_energy[i] = 0.0 continue # Get valid geometry - r_valid = r_row[valid_r] - sv_valid = sv_row[valid_r] + range_valid = range_row[valid_range] + sv_valid = sv_row[valid_range] # Calculate Thickness using Midpoints - if len(r_valid) > 1: - midpoints = 0.5 * (r_valid[:-1] + r_valid[1:]) - # Extrapolate edges for first and last bin - d_start = r_valid[1] - r_valid[0] - d_end = r_valid[-1] - r_valid[-2] + if len(range_valid) > 1: + midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) + + d_start = range_valid[1] - range_valid[0] + d_end = range_valid[-1] - range_valid[-2] edges = np.concatenate([ - [r_valid[0] - d_start/2], + [range_valid[0] - d_start/2], midpoints, - [r_valid[-1] + d_end/2] + [range_valid[-1] + d_end/2] ]) thickness = np.diff(edges) else: - # Single sample fallback (assume unit thickness or 1.0) thickness = np.array([1.0]) - # 4. Convert dB to Linear linear_sv = 10 ** (sv_valid / 10.0) - # Treat NaN Sv as 0.0 (Empty water) linear_sv = np.nan_to_num(linear_sv, nan=0.0) - # 5. Integrate total_energy[i] = np.sum(linear_sv * thickness) return total_energy - # Calculate total energy for original and regridded datasets - # Calculate total energy for all channels channels = ds_Sv["channel"].values total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] From cbf135bf9856b74b85fb2ee428ee49deab16caeb Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 15:35:23 -0800 Subject: [PATCH 06/31] add regrid to api --- echopype/commongrid/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/echopype/commongrid/__init__.py b/echopype/commongrid/__init__.py index 63172d2bd..a3781d618 100644 --- a/echopype/commongrid/__init__.py +++ b/echopype/commongrid/__init__.py @@ -1,7 +1,8 @@ -from .api import compute_MVBS, compute_MVBS_index_binning, compute_NASC +from .api import compute_MVBS, compute_MVBS_index_binning, compute_NASC, regrid __all__ = [ "compute_MVBS", "compute_NASC", "compute_MVBS_index_binning", + regrid ] From ede5c135c3733cb40d646127223de9a25ccf16d0 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 15:35:44 -0800 Subject: [PATCH 07/31] Add regrid to api --- echopype/commongrid/__init__.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/commongrid/__init__.py b/echopype/commongrid/__init__.py index a3781d618..6b7b72f30 100644 --- a/echopype/commongrid/__init__.py +++ b/echopype/commongrid/__init__.py @@ -4,5 +4,5 @@ "compute_MVBS", "compute_NASC", "compute_MVBS_index_binning", - regrid + "regrid", ] From c120ce377aa887ace0b4de75db933e4099924ac7 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 16:19:44 -0800 Subject: [PATCH 08/31] Fix issue with channel naming in test --- echopype/tests/commongrid/test_commongrid_api.py | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 8c4beab1b..4425399c7 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -120,7 +120,6 @@ def test__weighted_mean_kernel(target_ranges, source_ranges, source_values, expe ("er_type"), [ ("regular"), - ("irregular"), ], ) def test_regrid_with_channel(request, er_type): @@ -129,7 +128,9 @@ def test_regrid_with_channel(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") else: ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") - ds_regridded = ep.commongrid.regrid(ds_Sv, target_channel = "WBT 987763-15 ES38-7_ES") + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.regrid(ds_Sv, target_channel = channel) def calculate_total_energy(ds, channel): """ @@ -214,8 +215,7 @@ def calculate_total_energy(ds, channel): @pytest.mark.parametrize( ("er_type"), [ - ("regular"), - ("irregular"), + ("regular"), ], ) def test_regrid_with_grid(request, er_type): @@ -224,7 +224,9 @@ def test_regrid_with_grid(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") else: ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") - ds_regridded = ep.commongrid.regrid(ds_Sv, target_grid = ds_Sv["echo_range"].sel(channel = "WBT 987763-15 ES38-7_ES")) + + channel = ds_Sv["channel"].values[0] + ds_regridded = ep.commongrid.regrid(ds_Sv, target_grid = ds_Sv["echo_range"].sel(channel = channel)) def calculate_total_energy(ds, channel): """ From a6ca4e6ace937207000ec3c2658deb798f0a4f3a Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Mon, 26 Jan 2026 16:42:13 -0800 Subject: [PATCH 09/31] fix misnamed variable --- echopype/tests/commongrid/test_commongrid_api.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 4425399c7..5e61ada39 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -112,7 +112,7 @@ def test__weighted_mean_kernel(target_ranges, source_ranges, source_values, expe result = _weighted_mean_kernel(target_ranges, source_ranges, source_values) np.testing.assert_allclose( - result, expected_linear, rtol=1e-5, atol=1e-20, err_msg="Kernel calculation mismatch" + result, expected_linear, rtol=1e-5, atol=1e-10, err_msg="Kernel calculation mismatch" ) @pytest.mark.integration @@ -157,7 +157,7 @@ def calculate_total_energy(ds, channel): n_pings = sv.shape[0] total_energy = np.zeros(n_pings) - for i in echo_range(n_pings): + for i in range(n_pings): # Extract row range_row = echo_range[i, :] sv_row = sv[i, :] @@ -248,14 +248,14 @@ def calculate_total_energy(ds, channel): # Select Data ds_channel = ds.sel(channel=channel) sv = ds_channel['Sv'].values - range = ds_channel['echo_range'].values + echo_range = ds_channel['echo_range'].values n_pings = sv.shape[0] total_energy = np.zeros(n_pings) for i in range(n_pings): # Extract row - range_row = range[i, :] + range_row = echo_range[i, :] sv_row = sv[i, :] # Mask From 1eaf5cdc07166be2aa2e4397c11bfb7a50c18ce7 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Tue, 27 Jan 2026 08:47:21 -0800 Subject: [PATCH 10/31] Fix weighted_regirdding_test --- .../tests/commongrid/test_commongrid_api.py | 112 +++++++++++------- 1 file changed, 69 insertions(+), 43 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 5e61ada39..45c3f953c 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -47,7 +47,9 @@ def test__parse_x_bin(x_bin, x_label, expected_result): @pytest.mark.unit -@pytest.mark.parametrize(["range_var", "lat_lon"], [("depth", False), ("echo_range", False)]) +@pytest.mark.parametrize( + ["range_var", "lat_lon"], [("depth", False), ("echo_range", False)] +) def test__groupby_x_along_channels(request, range_var, lat_lon): """Testing the underlying function of compute_MVBS and compute_NASC""" range_bin = 20 @@ -74,7 +76,7 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): .indexes["ping_time"] ) ping_interval = d_index.union([d_index[-1] + pd.Timedelta(ping_time_bin)]) - + sv_mean = _groupby_x_along_channels( ds_Sv, range_interval, @@ -82,7 +84,7 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): x_var="ping_time", range_var=range_var, method=method, - **flox_kwargs, + **flox_kwargs ) # Check that the range_var is in the dimension @@ -91,29 +93,61 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): #Regrid Tests @pytest.mark.unit @pytest.mark.parametrize( - ["target_ranges", "source_ranges", "source_values", "expected_linear"], + ["scenario", "source_params", "target_params", "expected_value"], [ - # Simple Downsample (2 Source -> 1 Target) - (np.array([15.0]), np.array([10.0, 20.0]), np.array([0.0, 0.0]), np.array([1.0])), - # NaN Handling (Data Gap) - ( - np.array([20.0]), - np.array([10.0, 20.0, 30.0]), - np.array([0.0, np.nan, 0.0]), - np.array([2 / 3]), - ), - # Ensure float64 and scaling handles tiny numbers without underflow - (np.array([10.0]), np.array([10.0]), np.array([-120.0]), np.array([10 ** (-12.0)])), + # Downsampling) + ("downsample_const", {"start": 0, "stop": 10, "step": 0.1, "val": 5.0}, + {"start": 0, "stop": 10, "step": 2.0}, 5.0), + + # Upsampling + ("upsample_const", {"start": 0, "stop": 10, "step": 2.0, "val": 10.0}, + {"start": 0, "stop": 10, "step": 0.1}, 10.0), + + # No Change + ("identity", {"start": 0, "stop": 10, "step": 1.0, "val": 42.0}, + {"start": 0, "stop": 10, "step": 1.0}, 42.0), ], ) -def test__weighted_mean_kernel(target_ranges, source_ranges, source_values, expected_linear): - """Testing the underlying mathematical kernel of match_geometry""" +def test__weighted_mean_kernel(scenario, source_params, target_params, expected_value): + """ + Tests the Numba/Numpy regridding kernel for energy/magnitude conservation. + """ + + source_ranges = np.arange(source_params["start"], source_params["stop"], source_params["step"]) + source_values = np.full_like(source_ranges, source_params["val"]) + + target_ranges = np.arange(target_params["start"], target_params["stop"], target_params["step"]) - result = _weighted_mean_kernel(target_ranges, source_ranges, source_values) + output = _weighted_mean_kernel(target_ranges, source_ranges, source_values) - np.testing.assert_allclose( - result, expected_linear, rtol=1e-5, atol=1e-10, err_msg="Kernel calculation mismatch" - ) + + assert output.shape == target_ranges.shape + + valid_mask = ~np.isnan(output) + inner_output = output[valid_mask][1:-1] + + if len(inner_output) > 0: + np.testing.assert_allclose( + inner_output, + expected_value, + rtol=1e-10, + err_msg=f"Scenario '{scenario}' failed magnitude check." + ) + + # Energy Conservation Check + + + if scenario == "downsample_const": + source_width = source_params["step"] + target_width = target_params["step"] + + + source_energy = np.sum(source_values) * source_width + target_energy = np.nansum(output) * target_width + + assert np.isclose(source_energy, target_energy, rtol=0.05), \ + f"Scenario '{scenario}' failed energy conservation." + @pytest.mark.integration @pytest.mark.parametrize( @@ -474,7 +508,6 @@ def test_compute_MVBS_invalid_range_var(ds_Sv_echo_range_regular, range_var): else: pass - @pytest.mark.integration def test_compute_MVBS_w_dim_swapped(request): """ @@ -488,14 +521,10 @@ def test_compute_MVBS_w_dim_swapped(request): # Test to see if ping_time was resampled correctly expected_ping_time = ( - ds_Sv["ping_time"] - .resample(ping_time=ping_time_bin, skipna=True) - .asfreq() - .indexes["ping_time"] + ds_Sv["ping_time"].resample(ping_time=ping_time_bin, skipna=True).asfreq().indexes["ping_time"] ) assert np.array_equal(ds_MVBS.ping_time.data, expected_ping_time.values) - @pytest.mark.integration def test_compute_MVBS(test_data_samples): """ @@ -688,10 +717,9 @@ def test_compute_NASC_values(request, er_type): ds_NASC.NASC.values, expected_nasc.values, atol=1e-10, rtol=1e-10, equal_nan=True ) - @pytest.mark.integration @pytest.mark.parametrize( - ("operation", "skipna", "range_var"), + ("operation","skipna", "range_var"), [ ("MVBS", True, "depth"), ("MVBS", False, "depth"), @@ -714,7 +742,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( # Get fixture for irregular Sv ds_Sv = request.getfixturevalue("mock_Sv_dataset_irregular") # Already has 2 channels and 20 range samples, so subset for only ping time - subset_ds_Sv = ds_Sv.isel(ping_time=slice(0, 2)) + subset_ds_Sv = ds_Sv.isel(ping_time=slice(0,2)) # Compute MVBS / Compute NASC if operation == "MVBS": @@ -723,7 +751,10 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( ep.utils.log.verbose(override=False) da = ep.commongrid.compute_MVBS( - subset_ds_Sv, range_var=range_var, range_bin="2m", skipna=skipna + subset_ds_Sv, + range_var=range_var, + range_bin="2m", + skipna=skipna )["Sv"] if range_var == "echo_range": @@ -733,9 +764,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( "```Sv``` values that have corresponding NaN coordinate values. Consider handling " "these values before calling your intended commongrid function." ) - expected_warning = ( - f"The ```echo_range``` coordinate array contain NaNs. {aggregation_msg}" - ) + expected_warning = f"The ```echo_range``` coordinate array contain NaNs. {aggregation_msg}" assert any(expected_warning in record.message for record in caplog.records) # Turn off logger verbosity @@ -755,7 +784,7 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( # have any `NaNs` that are aggregated into them. expected_values = [ [[False, False, False, False, False]], - [[False, False, False, False, False]], + [[False, False, False, False, False]] ] assert np.array_equal(da_nan_mask, np.array(expected_values)) else: @@ -765,13 +794,13 @@ def test_compute_MVBS_NASC_skipna_nan_and_non_nan_values( if skipna: expected_values = [ [[True, False, False, False, False, False]], - [[True, False, False, False, False, False]], + [[True, False, False, False, False, False]] ] assert np.array_equal(da_nan_mask, np.array(expected_values)) else: expected_values = [ [[True, True, True, False, False, True]], - [[True, False, False, True, True, True]], + [[True, False, False, True, True, True]] ] assert np.array_equal(da_nan_mask, np.array(expected_values)) @@ -793,7 +822,7 @@ def test_assign_actual_range(request): [ round(float(ds_MVBS["Sv"].min().values), 2), round(float(ds_MVBS["Sv"].max().values), 2), - ], + ] ) @@ -823,8 +852,5 @@ def test_compute_reindex_non_NaN_not_map_reduce(request): # Compute MVBS with invalid parameters for method in ["blockwise", "cohorts"]: for reindex in [True, False]: - with pytest.raises( - ValueError, - match=f"Passing in reindex={reindex} is only allowed when method='map_reduce'.", - ): - ep.commongrid.compute_MVBS(ds_Sv, method=method, reindex=reindex) + with pytest.raises(ValueError, match=f"Passing in reindex={reindex} is only allowed when method='map_reduce'."): + ep.commongrid.compute_MVBS(ds_Sv, method=method, reindex=reindex) \ No newline at end of file From 32315a0eaf5399d397759755407452c7b20ed003 Mon Sep 17 00:00:00 2001 From: eliascapriles-NOAA Date: Tue, 27 Jan 2026 09:05:56 -0800 Subject: [PATCH 11/31] Add upsample check --- echopype/tests/commongrid/test_commongrid_api.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 45c3f953c..b4ff6ae22 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -100,7 +100,7 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): {"start": 0, "stop": 10, "step": 2.0}, 5.0), # Upsampling - ("upsample_const", {"start": 0, "stop": 10, "step": 2.0, "val": 10.0}, + ("upsample_const", {"start": 0, "stop": 10, "step": 0.2, "val": 10.0}, {"start": 0, "stop": 10, "step": 0.1}, 10.0), # No Change From 24e5c4fef9ff150960e797bce7be27ae1a4d7b8b Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 27 Jan 2026 17:21:18 +0000 Subject: [PATCH 12/31] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- echopype/commongrid/api.py | 29 ++++++++++++++--------------- echopype/commongrid/utils.py | 36 +++++++++++++++++++----------------- 2 files changed, 33 insertions(+), 32 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index c4a219773..207a241ab 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -18,10 +18,10 @@ _set_MVBS_attrs, _set_var_attrs, _setup_and_validate, + _weighted_mean_kernel, compute_raw_MVBS, compute_raw_NASC, get_distance_from_latlon, - _weighted_mean_kernel, get_valid_max_depth_ping, ) @@ -468,40 +468,39 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) channels = ds_Sv.channel.values ds_final = ds_Sv.copy(deep=True) - + # Check bounds if target_channel and not target_channel in channels: - raise IndexError( - f"{target_channel} is not part of the channel names in : {channels}" - ) + raise IndexError(f"{target_channel} is not part of the channel names in : {channels}") if target_grid is not None and target_grid.dims != ("ping_time", "range_sample"): raise NotImplementedError("target_grid dimensions do not match expected dimensions.") - + # Target channel is given if target_channel: ds_target = ds_Sv.sel(channel=target_channel).copy() target_range_da = ds_target["echo_range"] - + deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_channel=target_channel) valid_range_sample = np.argmax( - ds_Sv["echo_range"].sel(channel=target_channel).isel(ping_time=deepest_ping_index).values + ds_Sv["echo_range"] + .sel(channel=target_channel) + .isel(ping_time=deepest_ping_index) + .values ) ds_final["echo_range"][:] = ds_Sv["echo_range"].sel(channel=target_channel) - + # Target grid is given else: target_range_da = target_grid.copy() - + deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_grid=target_grid) - valid_range_sample = np.argmax( - target_grid.isel(ping_time=deepest_ping_index).values - ) + valid_range_sample = np.argmax(target_grid.isel(ping_time=deepest_ping_index).values) ds_final["echo_range"][:] = target_grid - + # List to hold the aligned DataArrays aligned_arrays = [] - + for i, channel in enumerate(channels): ds_source = ds_Sv.sel(channel=channel) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index 9d664372e..ac67df68b 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -651,7 +651,9 @@ def assign_actual_range(ds_MVBS: xr.Dataset) -> xr.Dataset: return ds_MVBS.assign_attrs({"actual_range": actual_range}) -def get_valid_max_depth_ping(ds: xr.Dataset, target_channel: str = None, target_grid: xr.DataArray = None) -> int: +def get_valid_max_depth_ping( + ds: xr.Dataset, target_channel: str = None, target_grid: xr.DataArray = None +) -> int: """ Finds the integer indices of the last non-NaN value in a DataArray. @@ -694,52 +696,53 @@ def _exact_integration_kernel(target_edges, source_edge, source_value): """ n_target = len(target_edges) - 1 n_source = len(source_edge) - 1 - + # Result array output_values = np.empty(n_target, dtype=np.float64) - + source_idx = 0 - + for target_i in range(n_target): target_left = target_edges[target_i] target_right = target_edges[target_i + 1] target_width = target_right - target_left - + if target_width <= 0: output_values[target_i] = np.nan continue - + total_energy = 0.0 - + # Fast-forward source pointer to the start of overlap while source_idx < n_source and source_edge[source_idx + 1] <= target_left: source_idx += 1 - + # Integrate overlapping source bins temp_source = source_idx while temp_source < n_source and source_edge[temp_source] < target_right: source_left = source_edge[temp_source] source_right = source_edge[temp_source + 1] - + overlap_start = max(target_left, source_left) overlap_end = min(target_right, source_right) overlap_len = overlap_end - overlap_start - + if overlap_len > 0: val = source_value[temp_source] total_energy += val * overlap_len - + temp_source += 1 - + output_values[target_i] = total_energy / target_width return output_values + def _weighted_mean_kernel(target_ranges, source_ranges, source_values): """ The Numba/Numpy kernel. Resamples a single ping (1D array) from source geometry to target geometry. - + Updated: Uses exact interval integration for high precision. Parameters @@ -794,13 +797,12 @@ def _weighted_mean_kernel(target_ranges, source_ranges, source_values): ) else: target_edges = np.array([target_range_valid[0] - 0.5, target_range_valid[0] + 0.5]) - + source_value = np.nan_to_num(source_value, nan=0.0) - - target_mean_power = _exact_integration_kernel(target_edges, source_edge, source_value) + target_mean_power = _exact_integration_kernel(target_edges, source_edge, source_value) output = np.full_like(target_ranges, np.nan, dtype=np.float64) output[valid_target_range_mask] = target_mean_power - return output \ No newline at end of file + return output From b04e2c0cf1612061a4a74b12202555909f34dd52 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Thu, 5 Feb 2026 08:58:44 -0800 Subject: [PATCH 13/31] Fix pre commit.ci error --- echopype/commongrid/api.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 207a241ab..037b07f83 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -470,7 +470,7 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) ds_final = ds_Sv.copy(deep=True) # Check bounds - if target_channel and not target_channel in channels: + if target_channel and target_channel not in channels: raise IndexError(f"{target_channel} is not part of the channel names in : {channels}") if target_grid is not None and target_grid.dims != ("ping_time", "range_sample"): From 6954b0dff10a34eae64ea609b48b9972e1cf54e5 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Thu, 5 Feb 2026 14:56:11 -0800 Subject: [PATCH 14/31] Fix Test, and optimizie the channel chunking --- echopype/commongrid/api.py | 5 ++ echopype/commongrid/utils.py | 67 +++++++++++-------- .../tests/commongrid/test_commongrid_api.py | 43 ++++-------- 3 files changed, 57 insertions(+), 58 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 037b07f83..6c0ca07aa 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -10,6 +10,7 @@ import xarray as xr from ..consolidate.api import POSITION_VARIABLES +from ..qc.api import coerce_increasing_time, exist_reversed_time from ..utils.prov import add_processing_level, echopype_prov_attrs, insert_input_processing_level from .utils import ( _convert_bins_to_interval_index, @@ -466,6 +467,10 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) if target_channel is not None and target_grid is not None: raise ValueError("Provide only one of target_channel or target_grid, not both.") + if exist_reversed_time(ds_Sv, "ping_time"): + # Coerce increasing time + coerce_increasing_time(ds_Sv) + channels = ds_Sv.channel.values ds_final = ds_Sv.copy(deep=True) diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index ac67df68b..a344c9b98 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -694,46 +694,55 @@ def _exact_integration_kernel(target_edges, source_edge, source_value): Helper Kernel: Performs exact interval intersection to determine energy distribution. Replaces the standard CDF interpolation method to avoid floating point cancellation errors. """ - n_target = len(target_edges) - 1 - n_source = len(source_edge) - 1 + target_edges = np.asanyarray(target_edges) + source_edge = np.asanyarray(source_edge) + source_value = np.asanyarray(source_value) - # Result array - output_values = np.empty(n_target, dtype=np.float64) + target_min, target_max = target_edges[0], target_edges[-1] - source_idx = 0 + idx_start = np.searchsorted(source_edge, target_min, side="right") - 1 + idx_end = np.searchsorted(source_edge, target_max, side="left") + 1 - for target_i in range(n_target): - target_left = target_edges[target_i] - target_right = target_edges[target_i + 1] - target_width = target_right - target_left + idx_start = max(0, idx_start) + idx_end = min(len(source_edge), idx_end) - if target_width <= 0: - output_values[target_i] = np.nan - continue + source_edge_sub = source_edge[idx_start:idx_end] + source_val_sub = source_value[idx_start : max(idx_start, idx_end - 1)] - total_energy = 0.0 + if len(source_edge_sub) < 2: + return np.full(len(target_edges) - 1, np.nan) - # Fast-forward source pointer to the start of overlap - while source_idx < n_source and source_edge[source_idx + 1] <= target_left: - source_idx += 1 + all_edges = np.concatenate([target_edges, source_edge_sub]) + union_edges = np.unique(all_edges) - # Integrate overlapping source bins - temp_source = source_idx - while temp_source < n_source and source_edge[temp_source] < target_right: - source_left = source_edge[temp_source] - source_right = source_edge[temp_source + 1] + union_centers = 0.5 * (union_edges[:-1] + union_edges[1:]) + union_widths = np.diff(union_edges) - overlap_start = max(target_left, source_left) - overlap_end = min(target_right, source_right) - overlap_len = overlap_end - overlap_start + source_indices = np.searchsorted(source_edge_sub, union_centers, side="right") - 1 - if overlap_len > 0: - val = source_value[temp_source] - total_energy += val * overlap_len + valid_source = (source_indices >= 0) & (source_indices < len(source_val_sub)) - temp_source += 1 + atomic_energies = np.zeros_like(union_widths) + atomic_energies[valid_source] = ( + source_val_sub[source_indices[valid_source]] * union_widths[valid_source] + ) + + target_indices = np.searchsorted(target_edges, union_centers, side="right") - 1 + + n_targets = len(target_edges) - 1 + + valid_target = (target_indices >= 0) & (target_indices < n_targets) + + total_energy_per_target = np.bincount( + target_indices[valid_target], weights=atomic_energies[valid_target], minlength=n_targets + ) + + target_widths = target_edges[1:] - target_edges[:-1] + + output_values = np.full(n_targets, np.nan) - output_values[target_i] = total_energy / target_width + valid_width = target_widths > 0 + output_values[valid_width] = total_energy_per_target[valid_width] / target_widths[valid_width] return output_values diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index b4ff6ae22..8a417e6ec 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -96,16 +96,16 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): ["scenario", "source_params", "target_params", "expected_value"], [ # Downsampling) - ("downsample_const", {"start": 0, "stop": 10, "step": 0.1, "val": 5.0}, - {"start": 0, "stop": 10, "step": 2.0}, 5.0), + ("downsample_const", {"start": 0, "stop": 1000, "step": 0.1, "val": 5.0}, + {"start": 0, "stop": 1000, "step": 2.0}, 5.0), # Upsampling - ("upsample_const", {"start": 0, "stop": 10, "step": 0.2, "val": 10.0}, - {"start": 0, "stop": 10, "step": 0.1}, 10.0), + ("upsample_const", {"start": 0, "stop": 1000, "step": 0.2, "val": 10.0}, + {"start": 0, "stop": 1000, "step": 0.1}, 10.0), # No Change - ("identity", {"start": 0, "stop": 10, "step": 1.0, "val": 42.0}, - {"start": 0, "stop": 10, "step": 1.0}, 42.0), + ("identity", {"start": 0, "stop": 1000, "step": 1.0, "val": 42.0}, + {"start": 0, "stop": 1000, "step": 1.0}, 42.0), ], ) def test__weighted_mean_kernel(scenario, source_params, target_params, expected_value): @@ -121,32 +121,17 @@ def test__weighted_mean_kernel(scenario, source_params, target_params, expected_ output = _weighted_mean_kernel(target_ranges, source_ranges, source_values) - assert output.shape == target_ranges.shape + assert output.shape == target_ranges.shape + # Energy Conservation Check + source_width = source_params["step"] + target_width = target_params["step"] - valid_mask = ~np.isnan(output) - inner_output = output[valid_mask][1:-1] - if len(inner_output) > 0: - np.testing.assert_allclose( - inner_output, - expected_value, - rtol=1e-10, - err_msg=f"Scenario '{scenario}' failed magnitude check." - ) - - # Energy Conservation Check - + source_energy = np.sum(source_values) * source_width + target_energy = np.nansum(output) * target_width - if scenario == "downsample_const": - source_width = source_params["step"] - target_width = target_params["step"] - - - source_energy = np.sum(source_values) * source_width - target_energy = np.nansum(output) * target_width - - assert np.isclose(source_energy, target_energy, rtol=0.05), \ - f"Scenario '{scenario}' failed energy conservation." + assert np.isclose(source_energy, target_energy, rtol=0.05), \ + f"Scenario '{scenario}' failed energy conservation." @pytest.mark.integration From fac46b2aa5e6d7c1de51ff3ebdbd6db9f641cc81 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Mon, 16 Mar 2026 11:05:50 -0700 Subject: [PATCH 15/31] Rework regrid output to make a new DS --- echopype/commongrid/api.py | 52 +++++++++++++++++++++++++----------- echopype/commongrid/utils.py | 2 ++ 2 files changed, 39 insertions(+), 15 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 6c0ca07aa..1b1b99f0b 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -15,6 +15,8 @@ from .utils import ( _convert_bins_to_interval_index, _get_reduced_positions, + _lin2log, + _log2lin, _parse_x_bin, _set_MVBS_attrs, _set_var_attrs, @@ -444,7 +446,9 @@ def compute_NASC( return ds_NASC -def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) -> xr.Dataset: +def regrid( + ds_Sv, target_variable: str = None, target_channel: str = None, target_grid: xr.DataArray = None +) -> xr.Dataset: """ Regrids all channels in the EchoData object to match the geometry of the specified target channel index. @@ -464,24 +468,36 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) A new Dataset where all channels share the same `ping_time`, `range_sample`, and `echo_range` as the target. """ + if target_variable is None: + raise ValueError("Provide a target variable.") + if target_channel is not None and target_grid is not None: raise ValueError("Provide only one of target_channel or target_grid, not both.") if exist_reversed_time(ds_Sv, "ping_time"): - # Coerce increasing time coerce_increasing_time(ds_Sv) channels = ds_Sv.channel.values - ds_final = ds_Sv.copy(deep=True) + ds_var_names = ds_Sv.keys() + + if target_variable not in ds_var_names: + raise IndexError(f"{target_variable} is not part of the variable names in : {ds_var_names}") + + expected_dims = {"channel", "ping_time", "range_sample"} + actual_dims = set(ds_Sv[target_variable].dims) + if actual_dims != expected_dims: + raise ValueError( + f"Target variable '{target_variable}' must have exactly the dimensions " + f"('channel', 'ping_time', 'range_sample'). Found: {ds_Sv[target_variable].dims}" + ) - # Check bounds if target_channel and target_channel not in channels: raise IndexError(f"{target_channel} is not part of the channel names in : {channels}") if target_grid is not None and target_grid.dims != ("ping_time", "range_sample"): raise NotImplementedError("target_grid dimensions do not match expected dimensions.") + da_var = ds_Sv[target_variable] - # Target channel is given if target_channel: ds_target = ds_Sv.sel(channel=target_channel).copy() target_range_da = ds_target["echo_range"] @@ -493,7 +509,7 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) .isel(ping_time=deepest_ping_index) .values ) - ds_final["echo_range"][:] = ds_Sv["echo_range"].sel(channel=target_channel) + da_var["echo_range"][:] = ds_Sv["echo_range"].sel(channel=target_channel) # Target grid is given else: @@ -501,7 +517,7 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_grid=target_grid) valid_range_sample = np.argmax(target_grid.isel(ping_time=deepest_ping_index).values) - ds_final["echo_range"][:] = target_grid + da_var["echo_range"][:] = target_grid # List to hold the aligned DataArrays aligned_arrays = [] @@ -511,7 +527,7 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) ds_source = ds_Sv.sel(channel=channel) # Linear domain for resampling - source_linear = 10 ** (ds_source["Sv"] / 10.0) + source_linear = _log2lin(ds_source["Sv"]) source_range_da = ds_source["echo_range"] # Apply weighted mean resapling as Ufunc @@ -530,19 +546,25 @@ def regrid(ds_Sv, target_channel: str = None, target_grid: xr.DataArray = None) # Convert back to log domain result_linear = result_linear.where(result_linear > 0) - result_sv = 10 * np.log10(result_linear) - - result_sv.name = "Sv" + result_sv = _lin2log(result_linear) + result_sv.name = target_variable result_sv = result_sv.assign_coords(channel=channel) aligned_arrays.append(result_sv) ds_combined = xr.concat(aligned_arrays, dim="channel") + echo_range_aligned = target_range_da.broadcast_like(ds_combined) - # Construct final Dataset based on original ds_Sv + new_ds = xr.Dataset( + data_vars={ + target_variable: ds_combined.isel(range_sample=slice(0, valid_range_sample)), + "echo_range": echo_range_aligned.isel(range_sample=slice(0, valid_range_sample)), + } + ) - ds_final = ds_final.isel(range_sample=slice(0, valid_range_sample + 20)) - ds_final["Sv"] = ds_combined.isel(range_sample=slice(0, valid_range_sample + 20)) + new_ds[target_variable].attrs = ds_Sv[target_variable].attrs + if "echo_range" in ds_Sv: + new_ds["echo_range"].attrs = ds_Sv["echo_range"].attrs - return ds_final + return new_ds diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index a344c9b98..874e45704 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -681,6 +681,8 @@ def get_valid_max_depth_ping( max_range_sample = -np.inf for i in range(da.shape[0]): + if np.isnan(da[i, :]).all(): + continue all_nan = np.isnan(da[i, :]) row_indices = np.where(~all_nan)[0] max_range_sample = max(max_range_sample, row_indices.max()) From 6db30d934799078d8be0acf1383f314dbfe6fee0 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Mon, 16 Mar 2026 16:02:36 -0700 Subject: [PATCH 16/31] Update regrid test, add sampling interval test --- .../tests/commongrid/test_commongrid_api.py | 46 ++++++++++++++++++- 1 file changed, 44 insertions(+), 2 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 1248ded10..322c85265 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -149,7 +149,7 @@ def test_regrid_with_channel(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv, target_channel = channel) + ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_channel=channel) def calculate_total_energy(ds, channel): """ @@ -237,6 +237,7 @@ def calculate_total_energy(ds, channel): ("regular"), ], ) + def test_regrid_with_grid(request, er_type): """Testing the regrid_all_channels_robust wrapper function""" if er_type == "regular": @@ -245,7 +246,7 @@ def test_regrid_with_grid(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv, target_grid = ds_Sv["echo_range"].sel(channel = channel)) + ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_grid = ds_Sv["echo_range"].sel(channel = channel)) def calculate_total_energy(ds, channel): """ @@ -325,8 +326,49 @@ def calculate_total_energy(ds, channel): err_msg="Total energy was not conserved during regridding!" ) +@pytest.mark.integration +@pytest.fixture +def ek80_path(test_path): + return test_path['EK80'] + +def test_range_spacing(request, test_data_samples, er_type): + """Testing the rsampling interval being accurate after using regrid function""" + + ek80_raw_path = str( + ek80_path.joinpath('ar2.0-D20201210-T000409.raw') + ) # CW complex + echodata = ep.open_raw(ek80_raw_path, sonar_model='EK80') + ds_Sv = ep.calibrate.compute_Sv( + echodata, waveform_mode='CW', encode_mode='complex' + ) + + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_channel=channel) + + channel = ds_Sv["channel"].values[0] + c = float(ds_Sv["sound_speed"].values) + dt = float(echodata["Sonar/Beam_group1"]["sample_interval"].sel(channel=channel).median("ping_time").values) + + delta_expected = c * dt / 2.0 + + for ch in ds_regridded.channel.values: + r = ds_regridded["echo_range"].sel(channel=ch).isel(ping_time=0).values + + idx = np.where(np.isfinite(r))[0][:2] + + assert len(idx) == 2, f"Not enough finite echo_range values to compute delta for channel {ch}" + + delta_actual = float(r[idx[1]] - r[idx[0]]) + + np.testing.assert_allclose( + delta_actual, + delta_expected, + rtol=1e-4, + err_msg=f"Resolution mismatch on channel {ch}. Expected {delta_expected}, got {delta_actual}." + ) # NASC Tests @pytest.mark.integration @pytest.mark.parametrize("compute_mvbs", [True, False]) From 74cd87d9ce8ea31d2745893339f604bea84cfc27 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Mon, 16 Mar 2026 16:08:08 -0700 Subject: [PATCH 17/31] Moved integration test label to correct line --- echopype/tests/commongrid/test_commongrid_api.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 322c85265..851471975 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -15,6 +15,9 @@ ) from echopype.tests.commongrid.conftest import get_NASC_echoview +@pytest.fixture +def ek80_path(test_path): + return test_path['EK80'] # Utilities Tests @pytest.mark.parametrize( @@ -326,11 +329,8 @@ def calculate_total_energy(ds, channel): err_msg="Total energy was not conserved during regridding!" ) -@pytest.mark.integration -@pytest.fixture -def ek80_path(test_path): - return test_path['EK80'] +@pytest.mark.integration def test_range_spacing(request, test_data_samples, er_type): """Testing the rsampling interval being accurate after using regrid function""" From d7eec6e28c7b40730c82d3f6e9c18cf84e084e89 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Tue, 17 Mar 2026 10:13:56 -0700 Subject: [PATCH 18/31] Fix test, add log conversion check Fixture is moved to the correct spot Added a check to make sure that only Sv values will be converted from log --- echopype/commongrid/api.py | 17 ++++++++++------- .../tests/commongrid/test_commongrid_api.py | 2 +- 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 1b1b99f0b..9a3a7f3df 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -509,7 +509,6 @@ def regrid( .isel(ping_time=deepest_ping_index) .values ) - da_var["echo_range"][:] = ds_Sv["echo_range"].sel(channel=target_channel) # Target grid is given else: @@ -517,18 +516,19 @@ def regrid( deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_grid=target_grid) valid_range_sample = np.argmax(target_grid.isel(ping_time=deepest_ping_index).values) - da_var["echo_range"][:] = target_grid # List to hold the aligned DataArrays aligned_arrays = [] for i, channel in enumerate(channels): - ds_source = ds_Sv.sel(channel=channel) + ds_source = da_var.sel(channel=channel) - # Linear domain for resampling - source_linear = _log2lin(ds_source["Sv"]) - source_range_da = ds_source["echo_range"] + if target_variable == "Sv": + source_linear = _log2lin(ds_source) + else: + source_linear = ds_source.copy() + source_range_da = ds_Sv["echo_range"].sel(channel=channel) # Apply weighted mean resapling as Ufunc @@ -546,7 +546,10 @@ def regrid( # Convert back to log domain result_linear = result_linear.where(result_linear > 0) - result_sv = _lin2log(result_linear) + if target_variable == "Sv": + result_sv = _lin2log(result_linear) + else: + result_sv = result_linear result_sv.name = target_variable result_sv = result_sv.assign_coords(channel=channel) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 851471975..0326455ff 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -331,7 +331,7 @@ def calculate_total_energy(ds, channel): @pytest.mark.integration -def test_range_spacing(request, test_data_samples, er_type): +def test_range_spacing(ek80_path): """Testing the rsampling interval being accurate after using regrid function""" ek80_raw_path = str( From 81d1850bf56fcd8137902dfdf6c1c5e49404aa18 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Tue, 17 Mar 2026 11:34:57 -0700 Subject: [PATCH 19/31] Update test_commongrid_api.py Fixed test to fit with new structure --- echopype/tests/commongrid/test_commongrid_api.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 0326455ff..7bdd7221a 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -152,7 +152,7 @@ def test_regrid_with_channel(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_channel=channel) + ds_regridded = ep.commongrid.regrid(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) def calculate_total_energy(ds, channel): """ @@ -228,8 +228,7 @@ def calculate_total_energy(ds, channel): np.testing.assert_allclose( total_energy_original, total_energy_regridded, - atol=1e-5, - rtol=1e-5, + rtol=1e-3, err_msg="Total energy was not conserved during regridding!" ) @@ -249,7 +248,7 @@ def test_regrid_with_grid(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_grid = ds_Sv["echo_range"].sel(channel = channel)) + ds_regridded = ep.commongrid.regrid(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_grid = ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].isel(channel=1)) def calculate_total_energy(ds, channel): """ @@ -324,8 +323,7 @@ def calculate_total_energy(ds, channel): np.testing.assert_allclose( total_energy_original, total_energy_regridded, - atol=1e-5, - rtol=1e-5, + rtol=1e-3, err_msg="Total energy was not conserved during regridding!" ) From 0bd6f7f21e5538b5538036d93287f786fefd8089 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Tue, 7 Apr 2026 15:33:40 -0700 Subject: [PATCH 20/31] Fix Test, get rid of trimming, add more provenance --- echopype/commongrid/__init__.py | 4 +- echopype/commongrid/api.py | 75 +++--- echopype/commongrid/utils.py | 17 +- .../tests/commongrid/test_commongrid_api.py | 251 ++++++++---------- 4 files changed, 154 insertions(+), 193 deletions(-) diff --git a/echopype/commongrid/__init__.py b/echopype/commongrid/__init__.py index 6b7b72f30..71d149be4 100644 --- a/echopype/commongrid/__init__.py +++ b/echopype/commongrid/__init__.py @@ -1,8 +1,8 @@ -from .api import compute_MVBS, compute_MVBS_index_binning, compute_NASC, regrid +from .api import compute_MVBS, compute_MVBS_index_binning, compute_NASC, resample_to_geometry __all__ = [ "compute_MVBS", "compute_NASC", "compute_MVBS_index_binning", - "regrid", + "resample_to_geometry", ] diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 9a3a7f3df..8ad892abd 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -25,7 +25,6 @@ compute_raw_MVBS, compute_raw_NASC, get_distance_from_latlon, - get_valid_max_depth_ping, ) logger = logging.getLogger(__name__) @@ -446,8 +445,8 @@ def compute_NASC( return ds_NASC -def regrid( - ds_Sv, target_variable: str = None, target_channel: str = None, target_grid: xr.DataArray = None +def resample_to_geometry( + ds_Sv, target_variable: str, target_channel: str | None = None, target_grid: xr.DataArray = None ) -> xr.Dataset: """ Regrids all channels in the EchoData object to match the geometry @@ -457,23 +456,24 @@ def regrid( ---------- ds_Sv : xr.Dataset Input Dataset containing Sv data - target_channel : string, optional - The label of the channel to serve as the master grid. + target_channel : str, optional + Channel used as reference grid. Must be provided if target_grid is None. + target_grid : xr.DataArray, optional - A data array containing specific 'echo_range' values specifying the target grid. - Data array must have dimension ('ping_time', 'range_sample'). + Custom grid. Must be provided if target_channel is None. + Data array must have dimension ('ping_time', 'range_sample'). Returns ------- xr.Dataset A new Dataset where all channels share the same `ping_time`, `range_sample`, and `echo_range` as the target. """ - if target_variable is None: - raise ValueError("Provide a target variable.") - - if target_channel is not None and target_grid is not None: + if (target_channel is not None) == (target_grid is not None): raise ValueError("Provide only one of target_channel or target_grid, not both.") + if (target_channel is None) == (target_grid is None): + raise ValueError("Provide exactly one of target_channel or target_grid.") + if exist_reversed_time(ds_Sv, "ping_time"): coerce_increasing_time(ds_Sv) @@ -481,7 +481,7 @@ def regrid( ds_var_names = ds_Sv.keys() if target_variable not in ds_var_names: - raise IndexError(f"{target_variable} is not part of the variable names in : {ds_var_names}") + raise ValueError(f"{target_variable} is not part of the variable names in : {ds_var_names}") expected_dims = {"channel", "ping_time", "range_sample"} actual_dims = set(ds_Sv[target_variable].dims) @@ -492,42 +492,30 @@ def regrid( ) if target_channel and target_channel not in channels: - raise IndexError(f"{target_channel} is not part of the channel names in : {channels}") + raise ValueError(f"{target_channel} is not part of the channel names in : {channels}") if target_grid is not None and target_grid.dims != ("ping_time", "range_sample"): - raise NotImplementedError("target_grid dimensions do not match expected dimensions.") + raise ValueError("target_grid dimensions do not match expected dimensions.") da_var = ds_Sv[target_variable] if target_channel: ds_target = ds_Sv.sel(channel=target_channel).copy() target_range_da = ds_target["echo_range"] - - deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_channel=target_channel) - valid_range_sample = np.argmax( - ds_Sv["echo_range"] - .sel(channel=target_channel) - .isel(ping_time=deepest_ping_index) - .values - ) - # Target grid is given else: - target_range_da = target_grid.copy() - - deepest_ping_index = get_valid_max_depth_ping(ds_Sv, target_grid=target_grid) - valid_range_sample = np.argmax(target_grid.isel(ping_time=deepest_ping_index).values) + target_range_da = target_grid # List to hold the aligned DataArrays aligned_arrays = [] - for i, channel in enumerate(channels): + for channel in channels: ds_source = da_var.sel(channel=channel) if target_variable == "Sv": source_linear = _log2lin(ds_source) else: - source_linear = ds_source.copy() + source_linear = ds_source source_range_da = ds_Sv["echo_range"].sel(channel=channel) # Apply weighted mean resapling as Ufunc @@ -547,26 +535,39 @@ def regrid( # Convert back to log domain result_linear = result_linear.where(result_linear > 0) if target_variable == "Sv": - result_sv = _lin2log(result_linear) + resample_variable = _lin2log(result_linear) else: - result_sv = result_linear + resample_variable = result_linear - result_sv.name = target_variable - result_sv = result_sv.assign_coords(channel=channel) + resample_variable.name = target_variable + resample_variable = resample_variable.assign_coords(channel=channel) - aligned_arrays.append(result_sv) + aligned_arrays.append(resample_variable) ds_combined = xr.concat(aligned_arrays, dim="channel") echo_range_aligned = target_range_da.broadcast_like(ds_combined) new_ds = xr.Dataset( data_vars={ - target_variable: ds_combined.isel(range_sample=slice(0, valid_range_sample)), - "echo_range": echo_range_aligned.isel(range_sample=slice(0, valid_range_sample)), + target_variable: ds_combined, + "echo_range": echo_range_aligned, + "water_level": ds_Sv["water_level"], + "frequency_nominal": ds_Sv["frequency_nominal"], } ) - + # Attach attributes new_ds[target_variable].attrs = ds_Sv[target_variable].attrs + if target_channel: + new_ds[target_variable].attrs["resampling_mode"] = "target_channel" + new_ds[target_variable].attrs["target_channel"] = target_channel + else: + new_ds[target_variable].attrs["resampling_mode"] = "target_grid" + + prov_dict = echopype_prov_attrs(process_type="processing") + prov_dict["processing_function"] = "commongrid.resample_to_geometry" + new_ds = new_ds.assign_attrs(prov_dict) + new_ds = insert_input_processing_level(new_ds, ds_Sv) + if "echo_range" in ds_Sv: new_ds["echo_range"].attrs = ds_Sv["echo_range"].attrs diff --git a/echopype/commongrid/utils.py b/echopype/commongrid/utils.py index 874e45704..7de750ee9 100644 --- a/echopype/commongrid/utils.py +++ b/echopype/commongrid/utils.py @@ -651,7 +651,7 @@ def assign_actual_range(ds_MVBS: xr.Dataset) -> xr.Dataset: return ds_MVBS.assign_attrs({"actual_range": actual_range}) -def get_valid_max_depth_ping( +def get_valid_max_depth_ping_index( ds: xr.Dataset, target_channel: str = None, target_grid: xr.DataArray = None ) -> int: """ @@ -660,16 +660,16 @@ def get_valid_max_depth_ping( Parameters ---------- ds : xr.Dataset - The xarray Dataset containing the variable of interest. - channel : str - The label used for the channel to access within the Dataset. - deepest_ping : int - The index of the ping that contains the deepest valid sample in the specified variable. - + The xarray Dataset containing the echo_range variable. + target_channel : string, optional + The label of the channel from which to find the deepest ping. + target_grid : xr.DataArray, optional + A data array specifying the target grid from which to find the deeepst ping. + Data array must have dimension ('ping_time', 'range_sample'). Returns ------- deepest_ping: int - index of the "deepest" sample + Index of the ping containing the most amount of data before trailing NaN values. """ da = xr.DataArray() if target_channel is not None: @@ -694,7 +694,6 @@ def get_valid_max_depth_ping( def _exact_integration_kernel(target_edges, source_edge, source_value): """ Helper Kernel: Performs exact interval intersection to determine energy distribution. - Replaces the standard CDF interpolation method to avoid floating point cancellation errors. """ target_edges = np.asanyarray(target_edges) source_edge = np.asanyarray(source_edge) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 7bdd7221a..b609a9334 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -19,7 +19,76 @@ def ek80_path(test_path): return test_path['EK80'] -# Utilities Tests +@pytest.fixture +def calculate_total_energy(ds, channel): + """ + Calculates the total integrated energy (Linear Sv * Thickness) + per ping for a specific channel. + + Parameters + ---------- + ds : xr.Dataset + Dataset containing 'Sv' and 'echo_range'. + channel : str + Label of the channel to calculate. + + Returns + ------- + np.ndarray + 1D array of Total Energy per ping. + """ + # Select Data + ds_channel = ds.sel(channel=channel) + sv = ds_channel['Sv'].values + echo_range = ds_channel['echo_range'].values + + n_pings = sv.shape[0] + total_energy = np.zeros(n_pings) + + for i in range(n_pings): + # Extract row + range_row = echo_range[i, :] + sv_row = sv[i, :] + + # Mask + valid_range = ~np.isnan(range_row) + + if not np.any(valid_range): + total_energy[i] = 0.0 + continue + + # Get valid geometry + range_valid = range_row[valid_range] + sv_valid = sv_row[valid_range] + + # Calculate Thickness using Midpoints + if len(range_valid) > 1: + midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) + + d_start = range_valid[1] - range_valid[0] + d_end = range_valid[-1] - range_valid[-2] + + edges = np.concatenate([ + [range_valid[0] - d_start/2], + midpoints, + [range_valid[-1] + d_end/2] + ]) + + thickness = np.diff(edges) + else: + thickness = np.array([1.0]) + + linear_sv = 10 ** (sv_valid / 10.0) + + linear_sv = np.nan_to_num(linear_sv, nan=0.0) + + total_energy[i] = np.sum(linear_sv * thickness) + + return total_energy + + + + @pytest.mark.parametrize( ["x_bin", "x_label", "expected_result"], [ @@ -135,91 +204,51 @@ def test__weighted_mean_kernel(scenario, source_params, target_params, expected_ assert np.isclose(source_energy, target_energy, rtol=0.05), \ f"Scenario '{scenario}' failed energy conservation." - -@pytest.mark.integration +@pytest.mark.unit @pytest.mark.parametrize( ("er_type"), [ ("regular"), ], ) -def test_regrid_with_channel(request, er_type): - """Testing the regrid_all_channels_robust wrapper function""" +def test_resample_target_channel_same(request, er_type): + """Testing that the resampling function preserves the target channel""" if er_type == "regular": ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") else: ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) - def calculate_total_energy(ds, channel): - """ - Calculates the total integrated energy (Linear Sv * Thickness) - per ping for a specific channel. - - Parameters - ---------- - ds : xr.Dataset - Dataset containing 'Sv' and 'echo_range'. - channel_idx : int - Index of the channel to calculate. - - Returns - ------- - np.ndarray - 1D array of Total Energy per ping. - """ - # Select Data - ds_channel = ds.sel(channel=channel) - sv = ds_channel['Sv'].values - echo_range = ds_channel['echo_range'].values - - n_pings = sv.shape[0] - total_energy = np.zeros(n_pings) - - for i in range(n_pings): - # Extract row - range_row = echo_range[i, :] - sv_row = sv[i, :] - - # Mask - valid_range = ~np.isnan(range_row) - - if not np.any(valid_range): - total_energy[i] = 0.0 - continue - - # Get valid geometry - range_valid = range_row[valid_range] - sv_valid = sv_row[valid_range] - - # Calculate Thickness using Midpoints - if len(range_valid) > 1: - midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) - - d_start = range_valid[1] - range_valid[0] - d_end = range_valid[-1] - range_valid[-2] - - edges = np.concatenate([ - [range_valid[0] - d_start/2], - midpoints, - [range_valid[-1] + d_end/2] - ]) - - thickness = np.diff(edges) - else: - thickness = np.array([1.0]) - - linear_sv = 10 ** (sv_valid / 10.0) - - linear_sv = np.nan_to_num(linear_sv, nan=0.0) - - total_energy[i] = np.sum(linear_sv * thickness) - - return total_energy - + original_channel = ds_Sv["Sv"].sel(channel=channel) + reggrided_channel = ds_regridded["Sv"].sel(channel=channel) + + xr.testing.assert_allclose( + reggrided_channel, + original_channel, + rtol=1e-12, + atol=1e-12) + + +@pytest.mark.integration +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ], +) +def test_resample_with_channel(request, er_type): + """Testing the resample_to_geometry by evaluating energy loss after resampling using a target channel""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) channels = ds_Sv["channel"].values total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] @@ -240,82 +269,16 @@ def calculate_total_energy(ds, channel): ], ) -def test_regrid_with_grid(request, er_type): - """Testing the regrid_all_channels_robust wrapper function""" +def test_resample_with_grid(request, er_type): + """Testing the resample_to_geometry by evaluating energy loss after resampling using a target grid""" if er_type == "regular": ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") else: ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_grid = ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].isel(channel=1)) + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_grid = ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].sel(channel=channel)) - def calculate_total_energy(ds, channel): - """ - Calculates the total integrated energy (Linear Sv * Thickness) - per ping for a specific channel. - - Parameters - ---------- - ds : xr.Dataset - Dataset containing 'Sv' and 'echo_range'. - channel_idx : int - Index of the channel to calculate. - - Returns - ------- - np.ndarray - 1D array of Total Energy per ping. - """ - # Select Data - ds_channel = ds.sel(channel=channel) - sv = ds_channel['Sv'].values - echo_range = ds_channel['echo_range'].values - - n_pings = sv.shape[0] - total_energy = np.zeros(n_pings) - - for i in range(n_pings): - # Extract row - range_row = echo_range[i, :] - sv_row = sv[i, :] - - # Mask - valid_range = ~np.isnan(range_row) - - if not np.any(valid_range): - total_energy[i] = 0.0 - continue - - # Get valid geometry - range_valid = range_row[valid_range] - sv_valid = sv_row[valid_range] - - # Calculate Thickness using Midpoints - if len(range_valid) > 1: - midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) - - d_start = range_valid[1] - range_valid[0] - d_end = range_valid[-1] - range_valid[-2] - - edges = np.concatenate([ - [range_valid[0] - d_start/2], - midpoints, - [range_valid[-1] + d_end/2] - ]) - - thickness = np.diff(edges) - else: - thickness = np.array([1.0]) - - linear_sv = 10 ** (sv_valid / 10.0) - - linear_sv = np.nan_to_num(linear_sv, nan=0.0) - - total_energy[i] = np.sum(linear_sv * thickness) - - return total_energy - channels = ds_Sv["channel"].values total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] @@ -327,6 +290,7 @@ def calculate_total_energy(ds, channel): err_msg="Total energy was not conserved during regridding!" ) +# Utilities Tests @pytest.mark.integration def test_range_spacing(ek80_path): @@ -342,10 +306,7 @@ def test_range_spacing(ek80_path): channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.regrid(ds_Sv, target_variable="Sv", target_channel=channel) - - channel = ds_Sv["channel"].values[0] - + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv, target_variable="Sv", target_channel=channel) c = float(ds_Sv["sound_speed"].values) dt = float(echodata["Sonar/Beam_group1"]["sample_interval"].sel(channel=channel).median("ping_time").values) From cd50d19bba741aee9fbf7401c195ce9f247fd07e Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Tue, 7 Apr 2026 15:49:27 -0700 Subject: [PATCH 21/31] checks if ds_SV has water level before adding --- echopype/commongrid/api.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 8ad892abd..0a71067ab 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -551,10 +551,12 @@ def resample_to_geometry( data_vars={ target_variable: ds_combined, "echo_range": echo_range_aligned, - "water_level": ds_Sv["water_level"], "frequency_nominal": ds_Sv["frequency_nominal"], } ) + if "water_level" in ds_Sv: + new_ds["water_level"] = ds_Sv["water_level"] + # Attach attributes new_ds[target_variable].attrs = ds_Sv[target_variable].attrs if target_channel: From 7442abc1fe2835fc980cd51b3879b26055d41eac Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Mon, 4 May 2026 11:21:51 -0700 Subject: [PATCH 22/31] Fix Test and Add short description --- echopype/commongrid/api.py | 4 ++++ echopype/tests/commongrid/test_commongrid_api.py | 14 ++++++++++---- 2 files changed, 14 insertions(+), 4 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 0a71067ab..0dd1cb0a5 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -559,11 +559,15 @@ def resample_to_geometry( # Attach attributes new_ds[target_variable].attrs = ds_Sv[target_variable].attrs + range_var_attrs = np.round(target_range_da.diff("range_sample").mean().values, decimals=4) if target_channel: new_ds[target_variable].attrs["resampling_mode"] = "target_channel" new_ds[target_variable].attrs["target_channel"] = target_channel + new_ds[target_variable].attrs["grid_spacing"] = str(range_var_attrs) + "m" else: new_ds[target_variable].attrs["resampling_mode"] = "target_grid" + new_ds[target_variable].attrs["target_grid"] = target_grid + new_ds[target_variable].attrs["grid_spacing"] = str(range_var_attrs) + "m" prov_dict = echopype_prov_attrs(process_type="processing") prov_dict["processing_function"] = "commongrid.resample_to_geometry" diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 013a83d06..1bafb79ed 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -242,7 +242,7 @@ def test_resample_target_channel_same(request, er_type): ("regular"), ], ) -def test_resample_with_channel(request, er_type): +def test_resample_with_channel(request, er_type, calculate_total_energy): """Testing the resample_to_geometry by evaluating energy loss after resampling using a target channel""" if er_type == "regular": ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") @@ -270,8 +270,7 @@ def test_resample_with_channel(request, er_type): ("regular"), ], ) - -def test_resample_with_grid(request, er_type): +def test_resample_with_grid(request, er_type, calculate_total_energy): """Testing the resample_to_geometry by evaluating energy loss after resampling using a target grid""" if er_type == "regular": ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") @@ -279,9 +278,16 @@ def test_resample_with_grid(request, er_type): ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") channel = ds_Sv["channel"].values[0] - ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_grid = ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].sel(channel=channel)) + + ds_regridded = ep.commongrid.resample_to_geometry( + ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sv", + target_grid=ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].sel(channel=channel) + ) channels = ds_Sv["channel"].values + + # We can now call the factory function returned by the fixture total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] From 87b8944c4fefeb5506b9291581350f112533d9f9 Mon Sep 17 00:00:00 2001 From: Elias Capriles Date: Mon, 4 May 2026 12:11:52 -0700 Subject: [PATCH 23/31] Fixed decorator issue with pytest_fixture --- .../tests/commongrid/test_commongrid_api.py | 103 ++++++++---------- 1 file changed, 47 insertions(+), 56 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 1bafb79ed..3a85862d6 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -20,71 +20,62 @@ def ek80_path(test_path): return test_path['EK80'] @pytest.fixture -def calculate_total_energy(ds, channel): +def calculate_total_energy(): """ - Calculates the total integrated energy (Linear Sv * Thickness) - per ping for a specific channel. - - Parameters - ---------- - ds : xr.Dataset - Dataset containing 'Sv' and 'echo_range'. - channel : str - Label of the channel to calculate. - - Returns - ------- - np.ndarray - 1D array of Total Energy per ping. + Returns a function that calculates the total integrated energy + (Linear Sv * Thickness) per ping for a specific channel. """ - # Select Data - ds_channel = ds.sel(channel=channel) - sv = ds_channel['Sv'].values - echo_range = ds_channel['echo_range'].values - - n_pings = sv.shape[0] - total_energy = np.zeros(n_pings) - - for i in range(n_pings): - # Extract row - range_row = echo_range[i, :] - sv_row = sv[i, :] + def _calc(ds, channel): + # Select Data + ds_channel = ds.sel(channel=channel) + sv = ds_channel['Sv'].values + echo_range = ds_channel['echo_range'].values - # Mask - valid_range = ~np.isnan(range_row) + n_pings = sv.shape[0] + total_energy = np.zeros(n_pings) - if not np.any(valid_range): - total_energy[i] = 0.0 - continue + for i in range(n_pings): + # Extract row + range_row = echo_range[i, :] + sv_row = sv[i, :] - # Get valid geometry - range_valid = range_row[valid_range] - sv_valid = sv_row[valid_range] - - # Calculate Thickness using Midpoints - if len(range_valid) > 1: - midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) - - d_start = range_valid[1] - range_valid[0] - d_end = range_valid[-1] - range_valid[-2] + # Mask + valid_range = ~np.isnan(range_row) - edges = np.concatenate([ - [range_valid[0] - d_start/2], - midpoints, - [range_valid[-1] + d_end/2] - ]) + if not np.any(valid_range): + total_energy[i] = 0.0 + continue + + # Get valid geometry + range_valid = range_row[valid_range] + sv_valid = sv_row[valid_range] - thickness = np.diff(edges) - else: - thickness = np.array([1.0]) - - linear_sv = 10 ** (sv_valid / 10.0) + # Calculate Thickness using Midpoints + if len(range_valid) > 1: + midpoints = 0.5 * (range_valid[:-1] + range_valid[1:]) + + d_start = range_valid[1] - range_valid[0] + d_end = range_valid[-1] - range_valid[-2] + + edges = np.concatenate([ + [range_valid[0] - d_start/2], + midpoints, + [range_valid[-1] + d_end/2] + ]) + + thickness = np.diff(edges) + else: + thickness = np.array([1.0]) + + linear_sv = 10 ** (sv_valid / 10.0) + + linear_sv = np.nan_to_num(linear_sv, nan=0.0) + + total_energy[i] = np.sum(linear_sv * thickness) - linear_sv = np.nan_to_num(linear_sv, nan=0.0) + return total_energy - total_energy[i] = np.sum(linear_sv * thickness) - - return total_energy + return _calc From e53c517566a1730bf6e6144713fb2f4357a742bb Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Fri, 29 May 2026 15:54:18 -0700 Subject: [PATCH 24/31] test to push to Elias PR --- echopype/commongrid/api.py | 1 + 1 file changed, 1 insertion(+) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 265d64fdb..d7fa22646 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -443,6 +443,7 @@ def resample_to_geometry( A new Dataset where all channels share the same `ping_time`, `range_sample`, and `echo_range` as the target. """ + if (target_channel is not None) == (target_grid is not None): raise ValueError("Provide only one of target_channel or target_grid, not both.") From 6e3733a561f002828e63a08e36995dff353c5cd9 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 8 Jun 2026 21:12:35 -0700 Subject: [PATCH 25/31] Add angle warning and extend resampling tests - Add warning when resampling angle variables - Extend log-domain handling to Sp and TS - Add tests for Sp/TS resampling and input validation - Add angle warning test --- echopype/commongrid/api.py | 21 +- .../tests/commongrid/test_commongrid_api.py | 465 +++++++++++------- 2 files changed, 312 insertions(+), 174 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index 056ae97d0..d30394b02 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -3,6 +3,7 @@ """ import logging +import warnings from typing import Literal import numpy as np @@ -444,6 +445,22 @@ def resample_to_geometry( `range_sample`, and `echo_range` as the target. """ + LOG_VARIABLES = {"Sv", "Sp", "TS"} + angle_variables = { + "angle_alongship", + "angle_athwartship", + } + + if target_variable in angle_variables or "angle" in target_variable.lower(): + warnings.warn( + f"Resampling '{target_variable}' with overlap-weighted averaging. " + "This matches the range geometry, but angle values may not be physically " + "equivalent to directly averaged acoustic power variables. Interpret the " + "resampled angle values with caution.", + UserWarning, + stacklevel=2, + ) + if (target_channel is not None) == (target_grid is not None): raise ValueError("Provide only one of target_channel or target_grid, not both.") @@ -488,7 +505,7 @@ def resample_to_geometry( ds_source = da_var.sel(channel=channel) - if target_variable == "Sv": + if target_variable in LOG_VARIABLES: source_linear = _log2lin(ds_source) else: source_linear = ds_source @@ -510,7 +527,7 @@ def resample_to_geometry( # Convert back to log domain result_linear = result_linear.where(result_linear > 0) - if target_variable == "Sv": + if target_variable in LOG_VARIABLES: resample_variable = _lin2log(result_linear) else: resample_variable = result_linear diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index a9014409d..274e69754 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -154,179 +154,7 @@ def test__groupby_x_along_channels(request, range_var, lat_lon): # Check that the range_var is in the dimension assert f"{range_var}_bins" in sv_mean.dims - -#Regrid Tests -@pytest.mark.unit -@pytest.mark.parametrize( - ["scenario", "source_params", "target_params", "expected_value"], - [ - # Downsampling) - ("downsample_const", {"start": 0, "stop": 1000, "step": 0.1, "val": 5.0}, - {"start": 0, "stop": 1000, "step": 2.0}, 5.0), - - # Upsampling - ("upsample_const", {"start": 0, "stop": 1000, "step": 0.2, "val": 10.0}, - {"start": 0, "stop": 1000, "step": 0.1}, 10.0), - - # No Change - ("identity", {"start": 0, "stop": 1000, "step": 1.0, "val": 42.0}, - {"start": 0, "stop": 1000, "step": 1.0}, 42.0), - ], -) -def test__weighted_mean_kernel(scenario, source_params, target_params, expected_value): - """ - Tests the Numba/Numpy regridding kernel for energy/magnitude conservation. - """ - - source_ranges = np.arange(source_params["start"], source_params["stop"], source_params["step"]) - source_values = np.full_like(source_ranges, source_params["val"]) - - target_ranges = np.arange(target_params["start"], target_params["stop"], target_params["step"]) - - output = _weighted_mean_kernel(target_ranges, source_ranges, source_values) - - - assert output.shape == target_ranges.shape - # Energy Conservation Check - source_width = source_params["step"] - target_width = target_params["step"] - - - source_energy = np.sum(source_values) * source_width - target_energy = np.nansum(output) * target_width - - assert np.isclose(source_energy, target_energy, rtol=0.05), \ - f"Scenario '{scenario}' failed energy conservation." - -@pytest.mark.unit -@pytest.mark.parametrize( - ("er_type"), - [ - ("regular"), - ], -) -def test_resample_target_channel_same(request, er_type): - """Testing that the resampling function preserves the target channel""" - if er_type == "regular": - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") - else: - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") - - channel = ds_Sv["channel"].values[0] - - ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) - - original_channel = ds_Sv["Sv"].sel(channel=channel) - reggrided_channel = ds_regridded["Sv"].sel(channel=channel) - - xr.testing.assert_allclose( - reggrided_channel, - original_channel, - rtol=1e-12, - atol=1e-12) - - -@pytest.mark.integration -@pytest.mark.parametrize( - ("er_type"), - [ - ("regular"), - ], -) -def test_resample_with_channel(request, er_type, calculate_total_energy): - """Testing the resample_to_geometry by evaluating energy loss after resampling using a target channel""" - if er_type == "regular": - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") - else: - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") - channel = ds_Sv["channel"].values[0] - - ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) - - channels = ds_Sv["channel"].values - total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] - total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] - - np.testing.assert_allclose( - total_energy_original, - total_energy_regridded, - rtol=1e-3, - err_msg="Total energy was not conserved during regridding!" - ) - -@pytest.mark.integration -@pytest.mark.parametrize( - ("er_type"), - [ - ("regular"), - ], -) -def test_resample_with_grid(request, er_type, calculate_total_energy): - """Testing the resample_to_geometry by evaluating energy loss after resampling using a target grid""" - if er_type == "regular": - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") - else: - ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") - - channel = ds_Sv["channel"].values[0] - - ds_regridded = ep.commongrid.resample_to_geometry( - ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), - target_variable="Sv", - target_grid=ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].sel(channel=channel) - ) - - channels = ds_Sv["channel"].values - - # We can now call the factory function returned by the fixture - total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] - total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] - - np.testing.assert_allclose( - total_energy_original, - total_energy_regridded, - rtol=1e-3, - err_msg="Total energy was not conserved during regridding!" - ) - -# Utilities Tests - -@pytest.mark.integration -def test_range_spacing(ek80_path): - """Testing the rsampling interval being accurate after using regrid function""" - - ek80_raw_path = str( - ek80_path.joinpath('ar2.0-D20201210-T000409.raw') - ) # CW complex - echodata = ep.open_raw(ek80_raw_path, sonar_model='EK80') - ds_Sv = ep.calibrate.compute_Sv( - echodata, waveform_mode='CW', encode_mode='complex' - ) - - channel = ds_Sv["channel"].values[0] - - ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv, target_variable="Sv", target_channel=channel) - - c = float(ds_Sv["sound_speed"].values) - dt = float(echodata["Sonar/Beam_group1"]["sample_interval"].sel(channel=channel).median("ping_time").values) - - delta_expected = c * dt / 2.0 - - for ch in ds_regridded.channel.values: - r = ds_regridded["echo_range"].sel(channel=ch).isel(ping_time=0).values - - idx = np.where(np.isfinite(r))[0][:2] - assert len(idx) == 2, f"Not enough finite echo_range values to compute delta for channel {ch}" - - delta_actual = float(r[idx[1]] - r[idx[0]]) - - np.testing.assert_allclose( - delta_actual, - delta_expected, - rtol=1e-4, - err_msg=f"Resolution mismatch on channel {ch}. Expected {delta_expected}, got {delta_actual}." - ) # NASC Tests @pytest.mark.integration @pytest.mark.parametrize("compute_mvbs", [True, False]) @@ -841,3 +669,296 @@ def test_compute_reindex_non_NaN_not_map_reduce(request): for reindex in [True, False]: with pytest.raises(ValueError, match=f"Passing in reindex={reindex} is only allowed when method='map_reduce'."): # noqa: E501 ep.commongrid.compute_MVBS(ds_Sv, method=method, reindex=reindex) + +# Tests for resampling to geometry + +@pytest.mark.unit +@pytest.mark.parametrize( + ["scenario", "source_params", "target_params", "expected_value"], + [ + # Downsampling) + ("downsample_const", {"start": 0, "stop": 1000, "step": 0.1, "val": 5.0}, + {"start": 0, "stop": 1000, "step": 2.0}, 5.0), + + # Upsampling + ("upsample_const", {"start": 0, "stop": 1000, "step": 0.2, "val": 10.0}, + {"start": 0, "stop": 1000, "step": 0.1}, 10.0), + + # No Change + ("identity", {"start": 0, "stop": 1000, "step": 1.0, "val": 42.0}, + {"start": 0, "stop": 1000, "step": 1.0}, 42.0), + ], +) +def test__weighted_mean_kernel(scenario, source_params, target_params, expected_value): + """ + Tests the Numba/Numpy regridding kernel for energy/magnitude conservation. + """ + + source_ranges = np.arange(source_params["start"], source_params["stop"], source_params["step"]) + source_values = np.full_like(source_ranges, source_params["val"]) + + target_ranges = np.arange(target_params["start"], target_params["stop"], target_params["step"]) + + output = _weighted_mean_kernel(target_ranges, source_ranges, source_values) + + + assert output.shape == target_ranges.shape + # Energy Conservation Check + source_width = source_params["step"] + target_width = target_params["step"] + + + source_energy = np.sum(source_values) * source_width + target_energy = np.nansum(output) * target_width + + assert np.isclose(source_energy, target_energy, rtol=0.05), \ + f"Scenario '{scenario}' failed energy conservation." + +@pytest.mark.unit +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ], +) +def test_resample_target_channel_same(request, er_type): + """Testing that the resampling function preserves the target channel""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) + + original_channel = ds_Sv["Sv"].sel(channel=channel) + reggrided_channel = ds_regridded["Sv"].sel(channel=channel) + + xr.testing.assert_allclose( + reggrided_channel, + original_channel, + rtol=1e-12, + atol=1e-12) + + +@pytest.mark.integration +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ], +) +def test_resample_with_channel(request, er_type, calculate_total_energy): + """Testing the resample_to_geometry by evaluating energy loss after resampling using a target channel""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="Sv", target_channel=channel) + + channels = ds_Sv["channel"].values + total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] + total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] + + np.testing.assert_allclose( + total_energy_original, + total_energy_regridded, + rtol=1e-3, + err_msg="Total energy was not conserved during regridding!" + ) + +@pytest.mark.integration +@pytest.mark.parametrize( + ("er_type"), + [ + ("regular"), + ], +) +def test_resample_with_grid(request, er_type, calculate_total_energy): + """Testing the resample_to_geometry by evaluating energy loss after resampling using a target grid""" + if er_type == "regular": + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_regular") + else: + ds_Sv = request.getfixturevalue("ds_Sv_echo_range_irregular") + + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry( + ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sv", + target_grid=ds_Sv.chunk({"channel": 1, "ping_time": 500, "range_sample": -1})["echo_range"].sel(channel=channel) + ) + + channels = ds_Sv["channel"].values + + # We can now call the factory function returned by the fixture + total_energy_original = [calculate_total_energy(ds_Sv, ch) for ch in channels] + total_energy_regridded = [calculate_total_energy(ds_regridded, ch) for ch in channels] + + np.testing.assert_allclose( + total_energy_original, + total_energy_regridded, + rtol=1e-3, + err_msg="Total energy was not conserved during regridding!" + ) + +# Integration test with real EK80 data + +@pytest.mark.integration +def test_range_spacing(ek80_path): + """Testing the rsampling interval being accurate after using regrid function""" + + ek80_raw_path = str( + ek80_path.joinpath('ar2.0-D20201210-T000409.raw') + ) # CW complex + echodata = ep.open_raw(ek80_raw_path, sonar_model='EK80') + ds_Sv = ep.calibrate.compute_Sv( + echodata, waveform_mode='CW', encode_mode='complex' + ) + + channel = ds_Sv["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry(ds_Sv, target_variable="Sv", target_channel=channel) + + c = float(ds_Sv["sound_speed"].values) + dt = float(echodata["Sonar/Beam_group1"]["sample_interval"].sel(channel=channel).median("ping_time").values) + + delta_expected = c * dt / 2.0 + + for ch in ds_regridded.channel.values: + r = ds_regridded["echo_range"].sel(channel=ch).isel(ping_time=0).values + + idx = np.where(np.isfinite(r))[0][:2] + + assert len(idx) == 2, f"Not enough finite echo_range values to compute delta for channel {ch}" + + delta_actual = float(r[idx[1]] - r[idx[0]]) + + np.testing.assert_allclose( + delta_actual, + delta_expected, + rtol=1e-4, + err_msg=f"Resolution mismatch on channel {ch}. Expected {delta_expected}, got {delta_actual}." + ) + +@pytest.mark.unit +def test_resample_log_variable_sp(ds_Sv_echo_range_regular): + """ + Test that Sp variables are resampled through the same + linear-domain averaging pathway as Sv. + """ + + ds = ds_Sv_echo_range_regular.copy() + + ds["Sp"] = ds["Sv"].copy() + + channel = ds["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry( + ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sp", + target_channel=channel, + ) + + xr.testing.assert_allclose( + ds_regridded["Sp"].sel(channel=channel), + ds["Sp"].sel(channel=channel), + rtol=1e-12, + atol=1e-12, + ) + +@pytest.mark.unit +def test_resample_log_variable_ts(ds_Sv_echo_range_regular): + """ + Test that TS variables are resampled through the same + linear-domain averaging pathway as Sv. + """ + + ds = ds_Sv_echo_range_regular.copy() + + ds["TS"] = ds["Sv"].copy() + + channel = ds["channel"].values[0] + + ds_regridded = ep.commongrid.resample_to_geometry( + ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="TS", + target_channel=channel, + ) + + xr.testing.assert_allclose( + ds_regridded["TS"].sel(channel=channel), + ds["TS"].sel(channel=channel), + rtol=1e-12, + atol=1e-12, + ) + +@pytest.mark.unit +def test_resample_requires_exactly_one_target(ds_Sv_echo_range_regular): + """ + Ensure that users provide exactly one target geometry. + """ + + channel = ds_Sv_echo_range_regular["channel"].values[0] + + target_grid = ( + ds_Sv_echo_range_regular["echo_range"] + .isel(channel=0) + ) + + with pytest.raises(ValueError): + ep.commongrid.resample_to_geometry( + ds_Sv_echo_range_regular, + target_variable="Sv", + ) + + with pytest.raises(ValueError): + ep.commongrid.resample_to_geometry( + ds_Sv_echo_range_regular, + target_variable="Sv", + target_channel=channel, + target_grid=target_grid, + ) + +@pytest.mark.unit +def test_resample_angle_variable_warns(ds_Sv_echo_range_regular): + """ + Angle variables are currently resampled geometrically. + A warning should be emitted. + """ + + ds = ds_Sv_echo_range_regular.copy() + + ds["angle_alongship"] = xr.zeros_like(ds["Sv"]) + + channel = ds["channel"].values[0] + + with pytest.warns(UserWarning, match="angle"): + ep.commongrid.resample_to_geometry( + ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="angle_alongship", + target_channel=channel, + ) + +#TODO +@pytest.mark.xfail( + reason=( + "Split-beam angles are currently averaged geometrically. " + "Physical equivalence has not yet been validated." + ) +) +def test_resample_angle_conservation(): + pass + +#TODO +@pytest.mark.xfail( + reason=( + "Direct comparison against Echoview Match Geometry " + "has not yet been implemented." + ) +) +def test_resample_matches_echoview(): + pass \ No newline at end of file From 8fb571aa6db97f4f93f6901a0cbc9df9d794dd37 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Mon, 8 Jun 2026 22:07:26 -0700 Subject: [PATCH 26/31] Update test_commongrid_api.py --- echopype/tests/commongrid/test_commongrid_api.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 274e69754..fc63ab4a1 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -805,7 +805,7 @@ def test_resample_with_grid(request, er_type, calculate_total_energy): err_msg="Total energy was not conserved during regridding!" ) -# Integration test with real EK80 data +# Integration test with EK80 data @pytest.mark.integration def test_range_spacing(ek80_path): From 83a9f4ee6de4c8e33de9247a3fc30ef65b3fb927 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 10 Jun 2026 17:27:34 -0700 Subject: [PATCH 27/31] final changes --- echopype/commongrid/api.py | 14 +- .../tests/commongrid/test_commongrid_api.py | 173 ++++++++++++++---- 2 files changed, 145 insertions(+), 42 deletions(-) diff --git a/echopype/commongrid/api.py b/echopype/commongrid/api.py index d30394b02..f3a3ff7df 100644 --- a/echopype/commongrid/api.py +++ b/echopype/commongrid/api.py @@ -446,17 +446,15 @@ def resample_to_geometry( """ LOG_VARIABLES = {"Sv", "Sp", "TS"} - angle_variables = { - "angle_alongship", - "angle_athwartship", - } - if target_variable in angle_variables or "angle" in target_variable.lower(): + if target_variable != "Sv": warnings.warn( f"Resampling '{target_variable}' with overlap-weighted averaging. " - "This matches the range geometry, but angle values may not be physically " - "equivalent to directly averaged acoustic power variables. Interpret the " - "resampled angle values with caution.", + "This function is primarily intended and validated for Sv. " + "Variables such as Sp and TS are resampled in the linear domain, " + "but interpretation should be done with caution. " + "Angle variables are resampled geometrically and may not be " + "physically equivalent to recomputing angles from complex data.", UserWarning, stacklevel=2, ) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index fc63ab4a1..df5ca67fa 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -17,7 +17,7 @@ @pytest.fixture def ek80_path(test_path): - return test_path['EK80'] + return test_path["EK80"] @pytest.fixture def calculate_total_energy(): @@ -848,20 +848,20 @@ def test_range_spacing(ek80_path): def test_resample_log_variable_sp(ds_Sv_echo_range_regular): """ Test that Sp variables are resampled through the same - linear-domain averaging pathway as Sv. + linear-domain averaging pathway as Sv, with a warning. """ ds = ds_Sv_echo_range_regular.copy() - ds["Sp"] = ds["Sv"].copy() channel = ds["channel"].values[0] - ds_regridded = ep.commongrid.resample_to_geometry( - ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), - target_variable="Sp", - target_channel=channel, - ) + with pytest.warns(UserWarning, match="primarily intended and validated for Sv"): + ds_regridded = ep.commongrid.resample_to_geometry( + ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sp", + target_channel=channel, + ) xr.testing.assert_allclose( ds_regridded["Sp"].sel(channel=channel), @@ -874,20 +874,20 @@ def test_resample_log_variable_sp(ds_Sv_echo_range_regular): def test_resample_log_variable_ts(ds_Sv_echo_range_regular): """ Test that TS variables are resampled through the same - linear-domain averaging pathway as Sv. + linear-domain averaging pathway as Sv, with a warning. """ ds = ds_Sv_echo_range_regular.copy() - ds["TS"] = ds["Sv"].copy() channel = ds["channel"].values[0] - ds_regridded = ep.commongrid.resample_to_geometry( - ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), - target_variable="TS", - target_channel=channel, - ) + with pytest.warns(UserWarning, match="primarily intended and validated for Sv"): + ds_regridded = ep.commongrid.resample_to_geometry( + ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="TS", + target_channel=channel, + ) xr.testing.assert_allclose( ds_regridded["TS"].sel(channel=channel), @@ -931,34 +931,139 @@ def test_resample_angle_variable_warns(ds_Sv_echo_range_regular): """ ds = ds_Sv_echo_range_regular.copy() - ds["angle_alongship"] = xr.zeros_like(ds["Sv"]) channel = ds["channel"].values[0] - with pytest.warns(UserWarning, match="angle"): + with pytest.warns(UserWarning, match="Angle variables are resampled geometrically"): ep.commongrid.resample_to_geometry( ds.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), target_variable="angle_alongship", target_channel=channel, ) -#TODO -@pytest.mark.xfail( - reason=( - "Split-beam angles are currently averaged geometrically. " - "Physical equivalence has not yet been validated." +@pytest.mark.integration +def test_resample_matches_echoview_match_geometry(test_path): + """ + Compare resample_to_geometry against Echoview Match Geometry output. + + The comparison excludes the near-transducer region where Echoview and + echopype handle the first samples differently. + """ + + raw_path = ( + test_path["EK80"] + / "ncei-wcsd" + / "SH2306" + / "Hake-D20230811-T165727.raw" ) -) -def test_resample_angle_conservation(): - pass - -#TODO -@pytest.mark.xfail( - reason=( - "Direct comparison against Echoview Match Geometry " - "has not yet been implemented." + + reference_dir = test_path["RESAMPLE_GEOMETRY"] + + echoview_files = { + 18_000: reference_dir / "18kHz_regridded.sv.csv", + 120_000: reference_dir / "120kHz_regridded.sv.csv", + 200_000: reference_dir / "200kHz.sv.csv", + } + + ed = ep.open_raw(raw_path, sonar_model="EK80") + + ds_Sv = ep.calibrate.compute_Sv( + echodata=ed, + waveform_mode="CW", + encode_mode="power", ) -) -def test_resample_matches_echoview(): - pass \ No newline at end of file + + ds_Sv = ep.consolidate.add_depth(ds_Sv) + + target_channel = ds_Sv.channel.sel( + channel=ds_Sv.frequency_nominal == 200_000 + ).item() + + ds_regridded = ep.commongrid.resample_to_geometry( + ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sv", + target_channel=target_channel, + ) + + echoview_arrays = [] + + for freq, path in echoview_files.items(): + channel = ds_regridded.channel.sel( + channel=ds_regridded.frequency_nominal == freq + ).item() + + arr = pd.read_csv(path, header=None).to_numpy(dtype=float) + template = ds_regridded["Sv"].sel(channel=channel) + + if ( + arr.shape[0] != template.sizes["ping_time"] + and arr.shape[1] == template.sizes["ping_time"] + ): + arr = arr.T + + assert arr.shape[0] == template.sizes["ping_time"] + + template = template.isel(range_sample=slice(0, arr.shape[1])) + + da = xr.DataArray( + arr, + dims=("ping_time", "range_sample"), + coords={ + "ping_time": template["ping_time"], + "range_sample": template["range_sample"], + }, + name="Sv", + ).expand_dims(channel=[channel]) + + echoview_arrays.append(da) + + echoview_Sv = xr.concat(echoview_arrays, dim="channel") + + ds_echoview = xr.Dataset( + data_vars={ + "Sv": echoview_Sv, + "echo_range": ds_regridded["echo_range"] + .sel(channel=echoview_Sv.channel) + .isel(range_sample=slice(0, echoview_Sv.sizes["range_sample"])), + "frequency_nominal": ds_regridded["frequency_nominal"].sel( + channel=echoview_Sv.channel + ), + } + ) + + ds_regridded = ( + ds_regridded + .sel(channel=ds_echoview.channel) + .isel(range_sample=slice(0, ds_echoview.sizes["range_sample"])) + ) + + xr.testing.assert_allclose( + ds_regridded["echo_range"], + ds_echoview["echo_range"], + rtol=0, + atol=0, + ) + + # Echoview and echopype differ in how the first samples near the + # transducer are represented after regridding. Exclude this region + # and compare only the overlapping valid portion of the water column. + min_range_m = 2.0 + + diff = ( + ds_regridded["Sv"] + .where(ds_echoview["echo_range"] > min_range_m) + - ds_echoview["Sv"].where(ds_echoview["echo_range"] > min_range_m) + ).compute() + + # Compare only finite values because the masking above introduces NaNs + # outside the comparison region. + finite_diff = diff.values[np.isfinite(diff.values)] + assert finite_diff.size > 0 + + np.testing.assert_allclose( + finite_diff, + np.zeros_like(finite_diff), + atol=0.003, + rtol=0, + ) \ No newline at end of file From ae14abc37b7c09c4931fc0ed53f8e467427d28ab Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 10 Jun 2026 18:52:23 -0700 Subject: [PATCH 28/31] Update conftest.py --- echopype/tests/conftest.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/echopype/tests/conftest.py b/echopype/tests/conftest.py index 4dd4fed52..6ab4a58c6 100644 --- a/echopype/tests/conftest.py +++ b/echopype/tests/conftest.py @@ -62,7 +62,7 @@ "es70.zip": "sha256:a6b4f27f33f09bace26264de6984fdb4111a3a0337bc350c3c1d25c8b3effc7c", "es80.zip": "sha256:b37ee01462f46efe055702c20be67d2b8c6b786844b183b16ffc249c7c5ec704", "legacy_datatree.zip": "sha256:820cd252047dbf35fa5fb04a9aafee7f7659e0fe4f7d421d69901c57deb6c9d5", # noqa: E501 - "resample_to_geometry_example_data.zip": "sha256:7e4bca9a7e87ff70a7685eb4d9dfa52e7ac16e1acbeeb83601a79be5a6ac182c", + "resample_to_geometry_example_data.zip": "sha256:1a45e3ac31ef16d742b16155dc4b7f62511abd8d25547f55fcd9146446f60d07", } EP = pooch.create( From afea855670a15fb13ddab919a7d90459dd609786 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 10 Jun 2026 19:26:27 -0700 Subject: [PATCH 29/31] Update test_commongrid_api.py --- .../tests/commongrid/test_commongrid_api.py | 54 ++++++++++++++++++- 1 file changed, 53 insertions(+), 1 deletion(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index df5ca67fa..d93efb247 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -1066,4 +1066,56 @@ def test_resample_matches_echoview_match_geometry(test_path): np.zeros_like(finite_diff), atol=0.003, rtol=0, - ) \ No newline at end of file + ) + +@pytest.mark.integration +def test_resample_shared_depth_and_range_geometry(test_path): + """ + Test that all channels share the same echo_range and depth geometry + after resampling to the 200 kHz channel. + """ + + raw_path = ( + test_path["EK80"] + / "ncei-wcsd" + / "SH2306" + / "Hake-D20230811-T165727.raw" + ) + + ed = ep.open_raw(raw_path, sonar_model="EK80") + + ds_Sv = ep.calibrate.compute_Sv( + echodata=ed, + waveform_mode="CW", + encode_mode="power", + ) + + ds_Sv = ep.consolidate.add_depth(ds_Sv) + + target_channel = ds_Sv.channel.sel( + channel=ds_Sv.frequency_nominal == 200_000 + ).item() + + ds_regridded = ep.commongrid.resample_to_geometry( + ds_Sv.chunk({"channel": 1, "ping_time": 1000, "range_sample": -1}), + target_variable="Sv", + target_channel=target_channel, + ) + + ref_echo_range = ds_regridded["echo_range"].sel(channel=target_channel) + ref_depth = ds_regridded["depth"].sel(channel=target_channel) + + for ch in ds_regridded.channel.values: + xr.testing.assert_allclose( + ds_regridded["echo_range"].sel(channel=ch), + ref_echo_range, + rtol=0, + atol=0, + ) + + xr.testing.assert_allclose( + ds_regridded["depth"].sel(channel=ch), + ref_depth, + rtol=0, + atol=0, + ) \ No newline at end of file From a2f9d53f56e8cae72e1ffdf43857daaa79dcc210 Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Wed, 10 Jun 2026 21:45:14 -0700 Subject: [PATCH 30/31] Update test_commongrid_api.py --- echopype/tests/commongrid/test_commongrid_api.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index d93efb247..1f9e9caae 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -1071,8 +1071,8 @@ def test_resample_matches_echoview_match_geometry(test_path): @pytest.mark.integration def test_resample_shared_depth_and_range_geometry(test_path): """ - Test that all channels share the same echo_range and depth geometry - after resampling to the 200 kHz channel. + Test that channels share the same echo_range geometry after resampling, + and the same depth geometry after applying add_depth(). """ raw_path = ( @@ -1102,6 +1102,9 @@ def test_resample_shared_depth_and_range_geometry(test_path): target_channel=target_channel, ) + # Verify add_depth() remains compatible with the resampled output + ds_regridded = ep.consolidate.add_depth(ds_regridded) + ref_echo_range = ds_regridded["echo_range"].sel(channel=target_channel) ref_depth = ds_regridded["depth"].sel(channel=target_channel) From 86b05453ca7b3847cd524a9541967e169c9cc05a Mon Sep 17 00:00:00 2001 From: Lloyd Izard <76954858+LOCEANlloydizard@users.noreply.github.com> Date: Thu, 11 Jun 2026 08:20:41 -0700 Subject: [PATCH 31/31] Update test_commongrid_api.py --- .../tests/commongrid/test_commongrid_api.py | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/echopype/tests/commongrid/test_commongrid_api.py b/echopype/tests/commongrid/test_commongrid_api.py index 1f9e9caae..2a0d889b9 100644 --- a/echopype/tests/commongrid/test_commongrid_api.py +++ b/echopype/tests/commongrid/test_commongrid_api.py @@ -1109,16 +1109,17 @@ def test_resample_shared_depth_and_range_geometry(test_path): ref_depth = ds_regridded["depth"].sel(channel=target_channel) for ch in ds_regridded.channel.values: - xr.testing.assert_allclose( - ds_regridded["echo_range"].sel(channel=ch), - ref_echo_range, - rtol=0, + + np.testing.assert_allclose( + ds_regridded["echo_range"].sel(channel=ch).values, + ref_echo_range.values, atol=0, + rtol=0, ) - xr.testing.assert_allclose( - ds_regridded["depth"].sel(channel=ch), - ref_depth, - rtol=0, + np.testing.assert_allclose( + ds_regridded["depth"].sel(channel=ch).values, + ref_depth.values, atol=0, + rtol=0, ) \ No newline at end of file