Allow external arrays in add depth; Add universal ping time alignment…

… function (#1369)

* move depth tests outside of consolidate_integration

* add space

* add align_to_ping_time function, and replace the other interp functions

* modify test to no longer check for ffill condition

* add test for align_ping_time_values

* allow external depth_offset and tilt arrays to add_depth

* add test for multidim depth offset and tilt value errors

* fix test_add_location

* remove redundant aligning to ping time logic

* add case where external_da time is None; be more lenient on dropping NaNs in harmonize env param times

* add tests for irregular expected_da inputs

* simplify glider test

* check that the correct extrapolation was done

* fix glider azfp test

* separate out interp from nmea selection

* fix loc name

* fix loc names

* fix estimate background noise upsampling

* Update echopype/utils/align.py

Co-authored-by: Wu-Jung Lee <leewujung@gmail.com>

* Update echopype/utils/align.py

Co-authored-by: Wu-Jung Lee <leewujung@gmail.com>

* Update echopype/consolidate/api.py

Co-authored-by: Wu-Jung Lee <leewujung@gmail.com>

* Update echopype/consolidate/api.py

Co-authored-by: Wu-Jung Lee <leewujung@gmail.com>

* set 3 senarios to 4 scenarios

---------

Co-authored-by: Wu-Jung Lee <leewujung@gmail.com>

echopype/calibrate/env_params.py

@@ -6,6 +6,7 @@
 
 from ..echodata import EchoData
 from ..utils import uwa
+from ..utils.align import align_to_ping_time
 from .cal_params import param2da
 
 ENV_PARAMS = (
@@ -63,13 +64,9 @@ def harmonize_env_param_time(
                 f"ping_time needs to be provided for comparison or interpolating {p.name}"
             )
 
-        # Direct assignment if all timestamps are identical (EK60 data)
-        if np.array_equal(p["time1"].data, ping_time.data):
-            return p.rename({"time1": "ping_time"})
-
-        # Interpolate `p` to `ping_time`
-        return (
-            p.dropna(dim="time1").interp(time1=ping_time).ffill(dim="ping_time").drop_vars("time1")
+        # Align array to ping time
+        return align_to_ping_time(
+            p.dropna(dim="time1", how="all"), "time1", ping_time, method="linear"
         )
     return p
 

echopype/clean/api.py

@@ -410,6 +410,7 @@ def estimate_background_noise(
 
     # Align noise `ping_time` to the first index of each coarsened `ping_time` bin
     noise = noise.assign_coords(ping_time=ping_num * np.arange(len(noise["ping_time"])))
+    power_cal = power_cal.assign_coords(ping_time=np.arange(len(power_cal["ping_time"])))
 
     # Limit max noise level
     noise = (
@@ -420,7 +421,9 @@ def estimate_background_noise(
 
     # Upsample noise to original ping time dimension
     Sv_noise = (
-        noise.reindex({"ping_time": power_cal["ping_time"]}, method="ffill")
+        noise.reindex({"ping_time": power_cal["ping_time"]}, method="ffill").assign_coords(
+            ping_time=ds_Sv["ping_time"]
+        )
         + spreading_loss
         + absorption_loss
     )
@@ -434,7 +437,7 @@ def remove_background_noise(
     ping_num: int,
     range_sample_num: int,
     background_noise_max: str = None,
-    SNR_threshold: float = 3.0,
+    SNR_threshold: float = "3.0dB",
 ) -> xr.Dataset:
     """
     Remove noise by using estimates of background noise

echopype/consolidate/api.py

@@ -1,5 +1,6 @@
 import datetime
 import pathlib
+from numbers import Number
 from pathlib import Path
 from typing import Optional, Union
 
@@ -9,6 +10,7 @@
 from ..calibrate.ek80_complex import get_filter_coeff
 from ..echodata import EchoData
 from ..echodata.simrad import retrieve_correct_beam_group
+from ..utils.align import align_to_ping_time
 from ..utils.io import get_file_format, open_source
 from ..utils.log import _init_logger
 from ..utils.prov import add_processing_level
@@ -17,7 +19,7 @@
     ek_use_platform_angles,
     ek_use_platform_vertical_offsets,
 )
-from .loc_utils import check_loc_time_dim_duplicates, check_loc_vars_validity, sel_interp
+from .loc_utils import check_loc_time_dim_duplicates, check_loc_vars_validity, sel_nmea
 from .split_beam_angle import get_angle_complex_samples, get_angle_power_samples
 
 logger = _init_logger(__name__)
@@ -65,8 +67,8 @@ def swap_dims_channel_frequency(ds: Union[xr.Dataset, str, pathlib.Path]) -> xr.
 def add_depth(
     ds: Union[xr.Dataset, str, pathlib.Path],
     echodata: Optional[Union[EchoData, str, pathlib.Path]] = None,
-    depth_offset: Optional[float] = None,
-    tilt: Optional[float] = None,
+    depth_offset: Optional[Union[Number, xr.DataArray]] = None,
+    tilt: Optional[Union[Number, xr.DataArray]] = None,
     downward: bool = True,
     use_platform_vertical_offsets: bool = False,
     use_platform_angles: bool = False,
@@ -82,11 +84,15 @@ def add_depth(
         Source Sv dataset to which a depth variable will be added.
     echodata : Optional[EchoData, str, pathlib.Path], default `None`
         `EchoData` object from which the `Sv` dataset originated.
-    depth_offset : Optional[float], default None
+    depth_offset : Optional[Union[Number, xr.DataArray]], default None
         Offset along the vertical (depth) dimension to account for actual transducer
         position in water, since `echo_range` is counted from transducer surface.
-    tilt : Optional[float], default None.
+        If provided as an xr.DataArray, it must have a single dimension corresponding to time.
+        The data array should represent the depth offset in meters at each time step.
+    tilt : Optional[Union[Number, xr.DataArray]], default None.
         Transducer tilt angle [degree]. 0 corresponds to a transducer pointing vertically.
+        If provided as an xr.DataArray, it must have a single dimension corresponding to time.
+        The data array should represent the tilt angle in degrees at each time step.
     downward : bool, default `True`
         The transducers point downward.
     use_platform_vertical_offsets: bool, default `False`
@@ -155,9 +161,21 @@ def add_depth(
 
     # Initialize transducer depth to 0.0 (no effect on depth)
     transducer_depth = 0.0
-    if depth_offset:
+    if isinstance(depth_offset, Number):
         # Set transducer depth to user specified depth offset
         transducer_depth = depth_offset
+    if isinstance(depth_offset, xr.DataArray):
+        # Will not accept data variables that don't have a single dimension
+        if len(depth_offset.dims) != 1:
+            raise ValueError(
+                "If depth_offset is passed in as an xr.DataArray, "
+                "it must contain a single dimension."
+            )
+        else:
+            # Align `depth_offset` to `ping_time`
+            transducer_depth = align_to_ping_time(
+                depth_offset, depth_offset.dims[0], ds["ping_time"]
+            )
     elif echodata and sonar_model in ["EK60", "EK80"]:
         if use_platform_vertical_offsets:
             # Compute transducer depth in EK systems using platform vertical offset data
@@ -167,9 +185,20 @@ def add_depth(
 
     # Initialize echo range scaling to 1.0 (no effect on depth)
     echo_range_scaling = 1.0
-    if tilt:
+    if isinstance(tilt, Number):
         # Set echo range scaling (angular scaling) using user specified tilt
         echo_range_scaling = np.cos(np.deg2rad(tilt))
+    if isinstance(tilt, xr.DataArray):
+        # Will not accept data variables that don't have a single dimension
+        if len(tilt.dims) != 1:
+            raise ValueError(
+                "If tilt is passed in as an xr.DataArray, it must contain a single dimension."
+            )
+        else:
+            # Align `tilt` to `ping_time`
+            echo_range_scaling = np.cos(
+                np.deg2rad(align_to_ping_time(tilt, tilt.dims[0], ds["ping_time"]))
+            )
     elif echodata and sonar_model in ["EK60", "EK80"]:
         if use_platform_angles:
             # Compute echo range scaling in EK systems using platform angle data
@@ -275,23 +304,18 @@ def add_location(
     # Copy dataset
     interp_ds = ds.copy()
 
-    # Interpolate location variables and place into `interp_ds`
-    interp_ds["latitude"] = sel_interp(
-        ds=ds,
-        echodata=echodata,
-        loc_name=lat_name,
-        time_dim_name=time_dim_name,
-        nmea_sentence=nmea_sentence,
-        datagram_type=datagram_type,
-    )
-    interp_ds["longitude"] = sel_interp(
-        ds=ds,
-        echodata=echodata,
-        loc_name=lon_name,
-        time_dim_name=time_dim_name,
-        nmea_sentence=nmea_sentence,
-        datagram_type=datagram_type,
-    )
+    # Select NMEA subset (if applicable) and interpolate location variables and place
+    # into `interp_ds`.
+    for loc_name, interp_loc_name in [(lat_name, "latitude"), (lon_name, "longitude")]:
+        loc_var = sel_nmea(
+            echodata=echodata,
+            loc_name=loc_name,
+            nmea_sentence=nmea_sentence,
+            datagram_type=datagram_type,
+        )
+        interp_ds[interp_loc_name] = align_to_ping_time(
+            loc_var, time_dim_name, ds["ping_time"], "linear"
+        )
 
     # Most attributes are attached automatically via interpolation
     # here we add the history

echopype/consolidate/ek_depth_utils.py

@@ -2,6 +2,7 @@
 import xarray as xr
 from scipy.spatial.transform import Rotation as R
 
+from ..utils.align import align_to_ping_time
 from ..utils.log import _init_logger
 
 logger = _init_logger(__name__)
@@ -26,31 +27,6 @@ def _check_and_log_nans(
             )
 
 
-def _var_time2_to_ping_time(var_with_time2, ping_time_da):
-    """
-    If `time2` does not differ from `var`, we rename `time2` to 'ping_time',
-    else interpolate `transducer_depth`'s `time2` dimension to `ping_time`.
-
-    Useful for handling EK60 and EK80 platform data:
-
-    EK80 will have platform variables with time2 dimension that does not match
-    Beam group ping time values, while EK60 will have time2 dimension that
-    matches Beam group ping time values.
-    """
-    if not ping_time_da.equals(var_with_time2["time2"].rename({"time2": "ping_time"})):
-        var_with_ping_time = var_with_time2.interp(
-            {"time2": ping_time_da},
-            method="nearest",
-            # More details for `fill_value` and `extrapolate` can be found here:
-            # https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.interp1d.html # noqa
-            kwargs={"fill_value": "extrapolate"},
-        ).drop_vars("time2")
-    else:
-        var_with_ping_time = var_with_time2.rename({"time2": "ping_time"})
-
-    return var_with_ping_time
-
-
 def ek_use_platform_vertical_offsets(
     platform_ds: xr.Dataset,
     ping_time_da: xr.DataArray,
@@ -72,7 +48,7 @@ def ek_use_platform_vertical_offsets(
     # Compute z translation for transducer position vs water level
     transducer_depth = transducer_offset_z - (water_level + vertical_offset)
 
-    return _var_time2_to_ping_time(transducer_depth, ping_time_da)
+    return align_to_ping_time(transducer_depth, "time2", ping_time_da)
 
 
 def ek_use_platform_angles(platform_ds: xr.Dataset, ping_time_da: xr.DataArray) -> xr.DataArray:
@@ -95,7 +71,7 @@ def ek_use_platform_angles(platform_ds: xr.Dataset, ping_time_da: xr.DataArray)
         echo_range_scaling, dims="time2", coords={"time2": platform_ds["time2"]}
     )
 
-    return _var_time2_to_ping_time(echo_range_scaling, ping_time_da)
+    return align_to_ping_time(echo_range_scaling, "time2", ping_time_da)
 
 
 def ek_use_beam_angles(

echopype/consolidate/loc_utils.py

@@ -118,16 +118,14 @@ def check_loc_time_dim_duplicates(echodata: EchoData, time_dim_name: str) -> Non
         )
 
 
-def sel_interp(
-    ds: xr.Dataset,
+def sel_nmea(
     echodata: EchoData,
     loc_name: str,
-    time_dim_name: str,
     nmea_sentence: Union[str, None] = None,
     datagram_type: Union[str, None] = None,
 ) -> xr.DataArray:
     """
-    Selection and interpolation for a location variable.
+    Select location subset for a location variable based on NMEA sentence.
 
     The selection logic is as follows, with 4 possible scenarios:
     Note here datagram_type = None is equivalent to datagram_type = NMEA
@@ -137,29 +135,15 @@ def sel_interp(
     3) If datagram_type is not None and nmea_sentence is None, then do nothing.
     4) If datagram_type is not None and nmea_sentence is not None, then raise ValueError since NMEA
        sentence selection can only be used on location variables from the NMEA datagram.
-
-    After selection logic, the location variable is then interpolated time-wise to match
-    that of the input dataset's time dimension.
     """
     # NMEA sentence selection if datagram_type is None (NMEA corresponds to None)
     if nmea_sentence and datagram_type is None:
-        position_var = echodata["Platform"][loc_name][
+        return echodata["Platform"][loc_name][
             echodata["Platform"]["sentence_type"] == nmea_sentence
         ]
     elif nmea_sentence and datagram_type is not None:
         raise ValueError(
             "If datagram_type is not `None`, then `nmea_sentence` cannot be specified."
         )
     else:
-        position_var = echodata["Platform"][loc_name]
-
-    if len(position_var) == 1:
-        # Propagate single, fixed-location coordinate
-        return xr.DataArray(
-            data=position_var.values[0] * np.ones(len(ds["ping_time"]), dtype=np.float64),
-            dims=["ping_time"],
-            attrs=position_var.attrs,
-        )
-    else:
-        # Values may be nan if there are ping_time values outside the time_dim_name range
-        return position_var.interp(**{time_dim_name: ds["ping_time"]})
+        return echodata["Platform"][loc_name]

echopype/tests/calibrate/test_calibrate.py

@@ -310,53 +310,12 @@ def test_compute_Sv_combined_ed_ping_time_extend_past_time1():
         "pH",
     ]
 
-    # Grab time variables
-    time1 = ed_combined["Environment"]["time1"]
-    ping_time = ed_combined["Sonar/Beam_group1"]["ping_time"]
-
     # Iterate through vars
     for env_var_name in environment_related_variable_names:
         env_var = ds_Sv[env_var_name]
         # Check that no NaNs exist
         assert not np.any(np.isnan(env_var.data))
 
-        # Check that all values past the max of time1 are ffilled with value
-        # that is time-wise closest to max of time1
-        if "channel" not in env_var.dims:
-            assert np.allclose(
-                np.unique(
-                    env_var.sel(
-                        ping_time=slice(
-                            time1.max(),
-                            ping_time.max()
-                        )
-                    ).data
-                ),
-                env_var.sel(ping_time=slice(time1.max())).data[-1]
-            )
-        else:
-            # Iterate through environment variable channels to do the same
-            # check as above per channel
-            for channel_index in range(len(env_var["channel"])):
-                assert np.allclose(
-                    np.unique(
-                        env_var.isel(
-                            channel=channel_index
-                        )
-                        .sel(
-                            ping_time=slice(
-                                time1.max(),
-                                ping_time.max()
-                            )
-                        ).data
-                    ),
-                    env_var.isel(
-                        channel=channel_index
-                    ).sel(
-                        ping_time=slice(time1.max())
-                    ).data[-1]
-                )
-
 
 @pytest.mark.parametrize(
     "raw_path, sonar_model, xml_path, waveform_mode, encode_mode",

echopype/tests/calibrate/test_env_params.py

@@ -30,6 +30,7 @@ def ek80_cal_path(test_path):
     return test_path["EK80_CAL"]
 
 
+@pytest.mark.unit
 def test_harmonize_env_param_time():
     # Scalar
     p = 10.05
@@ -67,16 +68,6 @@ def test_harmonize_env_param_time():
     assert (p_new.data == p.data).all()
 
     # ping_time target requires actual interpolation
-    ping_time_target = xr.DataArray(
-        data=[1, 2],
-        coords={
-            "ping_time": np.array(
-                ["2017-06-20T01:00:15", "2017-06-20T01:00:30"], dtype="datetime64[ns]"
-            )
-        },
-        dims=["ping_time"],
-    )
-
     ping_time_target = xr.DataArray(
         data=[1],
         coords={
@@ -133,6 +124,7 @@ def test_harmonize_env_param_time():
     # .all computes dask array under the hood
     assert (p_new.data == [0.5, 2880.5]).all()
 
+
 @pytest.mark.unit
 def test_harmonize_env_param_time_only_one_non_NaN_along_time1():
     # Create data array with time1 dimension with only 1 non-NaN value
@@ -145,6 +137,7 @@ def test_harmonize_env_param_time_only_one_non_NaN_along_time1():
     assert output_da == 1, "Output data array should just be 1"
     assert 'time1' not in output_da.dims, "```harmonize_env_param_time``` should have dropped 'time1' dimension"
 
+
 @pytest.mark.parametrize(
     ("user_dict", "channel", "out_dict"),
     [