feat: Add regrid functionality based on match geometry

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, resample_to_geometry
 
 __all__ = [
     "compute_MVBS",
     "compute_NASC",
     "compute_MVBS_index_binning",
+    "resample_to_geometry",
 ]

echopype/commongrid/api.py

@@ -3,21 +3,26 @@
 """
 
 import logging
+import warnings
 from typing import Literal
 
 import numpy as np
 import pandas as pd
 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,
     _get_reduced_positions,
+    _lin2log,
+    _log2lin,
     _parse_x_bin,
     _set_MVBS_attrs,
     _set_var_attrs,
     _setup_and_validate,
+    _weighted_mean_kernel,
     compute_raw_MVBS,
     compute_raw_NASC,
     get_distance_from_latlon,
@@ -416,5 +421,151 @@ def compute_NASC(
     return ds_NASC
 
 
-def regrid():
-    return 1
+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
+    of the specified target channel index.
+
+    Parameters
+    ----------
+    ds_Sv : xr.Dataset
+        Input Dataset containing Sv data
+    target_channel : str, optional
+    Channel used as reference grid. Must be provided if target_grid is None.
+
+    target_grid : xr.DataArray, optional
+    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.
+    """
+
+    LOG_VARIABLES = {"Sv", "Sp", "TS"}
+
+    if target_variable != "Sv":
+        warnings.warn(
+            f"Resampling '{target_variable}' with overlap-weighted averaging. "
+            "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,
+        )
+
+    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)
+
+    channels = ds_Sv.channel.values
+    ds_var_names = ds_Sv.keys()
+
+    if target_variable not 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)
+    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}"
+        )
+
+    if target_channel and target_channel not 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 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"]
+    # Target grid is given
+    else:
+        target_range_da = target_grid
+
+    # List to hold the aligned DataArrays
+    aligned_arrays = []
+
+    for channel in channels:
+
+        ds_source = da_var.sel(channel=channel)
+
+        if target_variable in LOG_VARIABLES:
+            source_linear = _log2lin(ds_source)
+        else:
+            source_linear = ds_source
+        source_range_da = ds_Sv["echo_range"].sel(channel=channel)
+
+        # 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)
+        if target_variable in LOG_VARIABLES:
+            resample_variable = _lin2log(result_linear)
+        else:
+            resample_variable = result_linear
+
+        resample_variable.name = target_variable
+        resample_variable = resample_variable.assign_coords(channel=channel)
+
+        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,
+            "echo_range": echo_range_aligned,
+            "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
+    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"
+    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
+
+    return new_ds

echopype/commongrid/utils.py

@@ -696,3 +696,170 @@ def ping_time_bin_parsing_and_conversion(ping_time_bin: str):
     )
     ping_time_bin_resunit_label = timedelta_units[ping_time_bin_resunit]["unitstr"]
     return ping_time_bin_resvalue, ping_time_bin_resunit_label
+
+
+def get_valid_max_depth_ping_index(
+    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.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        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 ping containing the most amount of data before trailing NaN 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
+
+    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())
+        if max_range_sample == row_indices.max():
+            deepest_ping = i
+    return deepest_ping
+
+
+def _exact_integration_kernel(target_edges, source_edge, source_value):
+    """
+    Helper Kernel: Performs exact interval intersection to determine energy distribution.
+    """
+    target_edges = np.asanyarray(target_edges)
+    source_edge = np.asanyarray(source_edge)
+    source_value = np.asanyarray(source_value)
+
+    target_min, target_max = target_edges[0], target_edges[-1]
+
+    idx_start = np.searchsorted(source_edge, target_min, side="right") - 1
+    idx_end = np.searchsorted(source_edge, target_max, side="left") + 1
+
+    idx_start = max(0, idx_start)
+    idx_end = min(len(source_edge), idx_end)
+
+    source_edge_sub = source_edge[idx_start:idx_end]
+    source_val_sub = source_value[idx_start : max(idx_start, idx_end - 1)]
+
+    if len(source_edge_sub) < 2:
+        return np.full(len(target_edges) - 1, np.nan)
+
+    all_edges = np.concatenate([target_edges, source_edge_sub])
+    union_edges = np.unique(all_edges)
+
+    union_centers = 0.5 * (union_edges[:-1] + union_edges[1:])
+    union_widths = np.diff(union_edges)
+
+    source_indices = np.searchsorted(source_edge_sub, union_centers, side="right") - 1
+
+    valid_source = (source_indices >= 0) & (source_indices < len(source_val_sub))
+
+    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)
+
+    valid_width = target_widths > 0
+    output_values[valid_width] = total_energy_per_target[valid_width] / target_widths[valid_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
+    ----------
+    target_ranges: np.ndarray
+        1D array of target bin centers to which the source data will be resampled.
+    source_ranges: np.ndarray
+        1D array of source bin centers representing the original data geometry.
+    source_values: np.ndarray
+        1D array of values corresponding to each source bin.
+
+    Returns
+    -------
+    output: np.ndarray
+        1D array of resampled values at the target bin centers.
+    """
+
+    valid_source_range_mask = ~np.isnan(source_ranges)
+    valid_target_range_mask = ~np.isnan(target_ranges)
+
+    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_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_edge = 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_edge = 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])
+
+    source_value = np.nan_to_num(source_value, nan=0.0)
+
+    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