Fix problems in applying compute_Sv on combined EchoData object (#…

…1328)

* remove redundancy in code, one has no effects, one causing error in harmonize env params

* ffill after interp time1 to ping_time, occurring for combined EchoData objects

* add bottleneck to requirements

* add test and drop time1

---------

Co-authored-by: ctuguinay <cmtuguinay2000@gmail.com>

echopype/calibrate/calibrate_ek.py

@@ -28,8 +28,6 @@ def __init__(self, echodata: EchoData, env_params, cal_params, ecs_file, **kwarg
 
         self.ed_beam_group = None  # will be assigned in child class
 
-        self.ed_beam_group = None  # will be assigned in child class
-
     def compute_echo_range(self, chan_sel: xr.DataArray = None):
         """
         Compute echo range for EK echosounders.

echopype/calibrate/env_params.py

@@ -68,7 +68,9 @@ def harmonize_env_param_time(
             return p.rename({"time1": "ping_time"})
 
         # Interpolate `p` to `ping_time`
-        return p.dropna(dim="time1").interp(time1=ping_time)
+        return (
+            p.dropna(dim="time1").interp(time1=ping_time).ffill(dim="ping_time").drop_vars("time1")
+        )
     return p
 
 

echopype/calibrate/range.py

@@ -170,12 +170,6 @@ def compute_range_EK(
         else echodata[ed_beam_group].sel(channel=chan_sel)
     )
 
-    # Harmonize sound_speed time1 and Beam_groupX ping_time
-    sound_speed = harmonize_env_param_time(
-        p=sound_speed,
-        ping_time=beam.ping_time,
-    )
-
     # Range in meters, not modified for TVG compensation
     range_meter = beam["range_sample"] * beam["sample_interval"] * sound_speed / 2
     # make order of dims conform with the order of backscatter data

echopype/tests/calibrate/test_calibrate.py

@@ -269,6 +269,95 @@ def test_compute_Sv_ek80_CW_complex_BB_complex(ek80_cal_path, ek80_path):
 
 
 @pytest.mark.integration
+def test_compute_Sv_combined_ed_ping_time_extend_past_time1():
+    """
+    Test computing combined Echodata object when ping time dimension in Beam group
+    extends past time1 dimension in Environment group.
+    The output Sv dataset should not have any NaN values within any of the time1
+    variables derived from the Environment group. Additionally, the output Sv dataset
+    should not contain the time1 dimension.
+    """
+    # Parse RAW files and combine Echodata objects
+    raw_list = [
+        "echopype/test_data/ek80/pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T001150-0.raw",
+        "echopype/test_data/ek80/pifsc_saildrone/SD_TPOS2023_v03-Phase0-D20230530-T002350-0.raw",
+    ]
+    ed_list = []
+    for raw_file in raw_list:
+        ed = ep.open_raw(raw_file, sonar_model="EK80")
+        # Modify environment variables so that they are non-uniform across Echodata objects
+        ed["Environment"]["acidity"].values = [np.random.uniform(low=7.9, high=8.1)]
+        ed["Environment"]["salinity"].values = [np.random.uniform(low=34.0, high=35.0)]
+        ed["Environment"]["temperature"].values = [np.random.uniform(low=25.0, high=26.0)]
+        ed_list.append(ed)
+    ed_combined = ep.combine_echodata(ed_list)
+
+    # Compute Sv
+    ds_Sv = ep.calibrate.compute_Sv(
+        ed_combined,
+        waveform_mode="CW",
+        encode_mode="complex"
+    )
+
+    # Check that Sv doesn't have time1 coordinate
+    assert "time1" not in ds_Sv.coords
+
+    # Define environment related variables
+    environment_related_variable_names = [
+        "sound_absorption",
+        "temperature",
+        "salinity",
+        "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",
     [

requirements.txt

@@ -1,3 +1,4 @@
+bottleneck
 dask[array,distributed]
 jinja2
 netCDF4>1.6