Support for EK80 Broadband Complex Multiplex Data (#1302)

* add modifications to work with ek80 bb complex multiplex data

* let convolve choose fft vs dicrete, use drop_
var

* add test for bb complex multiplex ek80 data

* use drop

* add test for env params when dropna has only 1 value

* revert small change

* add changes for convolving nans efficiently and handling nans in gain compensation

* use floats not ints in convolve, move convolve test to appropriate test file

* NaN to zeros -> zero out NaNs

* remove mark test

* unit test assertion message correction

* use the original lambda function and do the nan zero nan operation in compress pulse, place the atol 1e-5 back

* fix lazy loading

Author: ctuguinay

echopype/calibrate/calibrate_ek.py

@@ -473,6 +473,10 @@ def _get_B_theta_phi_m(self):
         ) ** 2
         B_theta_phi_m = 0.5 * 6.0206 * (fac_along + fac_athwart - 0.18 * fac_along * fac_athwart)
 
+        # Zero out NaNs that appear due to 1) multiplex ping patterns and 2) single beam transducer
+        # systems that don't contain angle offset information.
+        B_theta_phi_m = B_theta_phi_m.fillna(0)
+
         return B_theta_phi_m
 
     def _cal_complex_samples(self, cal_type: str) -> xr.Dataset:

echopype/calibrate/ek80_complex.py

@@ -226,11 +226,10 @@ def get_transmit_signal(
     # but keeping this here for use as standalone function
     if waveform_mode == "BB" and np.all(beam["transmit_type"] == "CW"):
         raise TypeError("File does not contain BB mode complex samples!")
-
     # Generate all transmit replica
     y_all = {}
     y_time_all = {}
-    # TODO: expand to deal with the case with varying tx param across ping_time
+    # TODO: expand to deal with the case with varying non-NaN tx param across ping_time
     tx_param_names = [
         "transmit_duration_nominal",
         "slope",
@@ -240,16 +239,18 @@ def get_transmit_signal(
     for ch in beam["channel"].values:
         tx_params = {}
         for p in tx_param_names:
-            tx_params[p] = np.unique(beam[p].sel(channel=ch))
+            # Extract beam values and filter out NaNs
+            beam_values = np.unique(beam[p].sel(channel=ch))
+            # Filter out NaN values
+            beam_values_without_nan = beam_values[~np.isnan(beam_values)]
+            tx_params[p] = beam_values_without_nan
             if tx_params[p].size != 1:
                 raise TypeError("File contains changing %s!" % p)
         fs_chan = fs.sel(channel=ch).data if isinstance(fs, xr.DataArray) else fs
         tx_params["fs"] = fs_chan
         y_ch, _ = tapered_chirp(**tx_params)
-
         # Filter and decimate chirp template
         y_ch, y_tmp_time = filter_decimate_chirp(coeff_ch=coeff[ch], y_ch=y_ch, fs=fs_chan)
-
         # Fill into output dict
         y_all[ch] = y_ch
         y_time_all[ch] = y_tmp_time
@@ -278,9 +279,21 @@ def compress_pulse(backscatter: xr.DataArray, chirp: Dict) -> xr.DataArray:
         # Select channel `chan` and drop the specific beam dimension if all of the values are nan.
         backscatter_chan = backscatter.sel(channel=chan).dropna(dim="beam", how="all")
 
+        # Create NaN mask
+        # If `backscatter_chan` is lazy loaded, then `nan_mask` too will be lazy loaded.
+        nan_mask = np.isnan(backscatter_chan)
+
+        # Zero out backscatter NaN values
+        # If `nan_mask` is lazy loaded, then resulting `backscatter_chan` will be lazy loaded.
+        backscatter_chan = xr.where(nan_mask, 0.0 + 0j, backscatter_chan)
+
+        # Extract transmit values
         tx = chirp[str(chan.values)]
+
+        # Compute complex conjugate of transmit values and reverse order of elements
         replica = np.flipud(np.conj(tx))
 
+        # Apply convolve on backscatter (along range sample dimension) and replica
         pc = xr.apply_ufunc(
             lambda m: (signal.convolve(m, replica, mode="full")[replica.size - 1 :]),
             backscatter_chan,
@@ -291,6 +304,11 @@ def compress_pulse(backscatter: xr.DataArray, chirp: Dict) -> xr.DataArray:
             output_dtypes=[np.complex64],
         ).compute()
 
+        # Restore NaN values in the resulting array.
+        # Computing of `nan_mask` here is necessary in the case when `nan_mask` is lazy loaded
+        # or else the resulting `pc` will also be lazy loaded.
+        pc = xr.where(nan_mask.compute(), np.nan, pc)
+
         pc_all.append(pc)
 
     pc_all = xr.concat(pc_all, dim="channel")

echopype/calibrate/env_params.py

@@ -50,10 +50,13 @@ def harmonize_env_param_time(
         if "time1" not in p.coords:
             return p
 
-        # If there's only 1 time1 value,
-        # or if after dropping NaN there's only 1 time1 value
-        if p["time1"].size == 1 or p.dropna(dim="time1").size == 1:
-            return p.dropna(dim="time1").squeeze(dim="time1").drop("time1")
+        # If there's only 1 time1 value:
+        if p["time1"].size == 1:
+            return p.squeeze(dim="time1").drop_vars("time1")
+
+        # If after dropping NaN along time1 dimension there's only 1 time1 value:
+        if p.dropna(dim="time1").size == 1:
+            return p.dropna(dim="time1").squeeze(dim="time1").drop_vars("time1")
 
         if ping_time is None:
             raise ValueError(

echopype/tests/calibrate/test_calibrate_ek80.py

@@ -177,7 +177,6 @@ def test_ek80_BB_range(ek80_cal_path, ek80_ext_path):
     pyel_vals = pyel_BB_p_data["range"]
     assert np.allclose(pyel_vals, ep_vals)
 
-
 @pytest.mark.parametrize(
     ("raw_data_path,raw_file_name,pyecholab_data_path,pyecholab_file_path, dask_array"),
     [
@@ -253,7 +252,6 @@ def test_ek80_BB_power_from_complex(
     )
     assert np.allclose(pyel_vals[idx_to_cmp], ep_vals[idx_to_cmp])
 
-
 @pytest.mark.parametrize(
     ("raw_data_path,raw_file_name,pyecholab_data_path,pyecholab_file_path, dask_array"),
     [
@@ -371,10 +369,9 @@ def test_ek80_BB_power_echoview(ek80_path):
     # Below only compare the first ping
     ev_vals = df_real.values[:, :]
     ep_vals = pc_mean.values.real[:, :]
-    assert np.allclose(ev_vals[:, 69:8284], ep_vals[:, 69:], atol=1e-4)
-    assert np.allclose(ev_vals[:, 90:8284], ep_vals[:, 90:], atol=1e-5)
+    assert np.allclose(ev_vals[:, 69:], ep_vals[:, 69:], atol=1e-4)
+    assert np.allclose(ev_vals[:, 90:], ep_vals[:, 90:], atol=1e-5)
 
-    
 def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path):
     ek80_raw_path = str(ek80_path.joinpath("D20230804-T083032.raw"))
     ed = ep.open_raw(ek80_raw_path, sonar_model="EK80")
@@ -390,3 +387,67 @@ def test_ek80_CW_complex_Sv_receiver_sampling_freq(ek80_path):
     # receiver_sampling_frequency substituted with default value in compute_Sv
     assert ds_Sv["receiver_sampling_frequency"] is not None
     assert np.allclose(ds_Sv["receiver_sampling_frequency"].data, 1500000)
+
+
+@pytest.mark.parametrize(
+    ("raw_data_path,target_channel_ping_pattern"),
+    [
+        (
+            "NYOS2105-D20210525-T213648.raw",
+            # (3,4), 2, 1, (3,4), 2, 1, etc...
+            np.array(
+                [
+                    [np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan],
+                    [np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan],
+                    [ 1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1.],
+                    [ 1., np.nan, np.nan,  1., np.nan, np.nan,  1., np.nan, np.nan,  1.]
+                ]
+            ),
+        ),
+        (
+            "DRIX08-D20231003-T120051.raw",
+            # 2, 1, 2, 1, 2, 1, etc...
+            np.array(
+                [
+                    [np.nan,  1., np.nan, 1,  np.nan, 1, np.nan,  1., np.nan, 1],
+                    [1, np.nan,  1., np.nan, 1,  np.nan, 1, np.nan,  1., np.nan],
+                ]
+            ),
+        ),
+    ],
+)
+@pytest.mark.integration
+def test_ek80_BB_complex_multiplex_NaNs_and_non_NaNs(raw_data_path, target_channel_ping_pattern):
+    # Extract bb complex multiplex EK80 data
+    ed = ep.open_raw(f"echopype/test_data/ek80_bb_complex_multiplex/{raw_data_path}", sonar_model="EK80")
+
+    # Compute Sv
+    ds_Sv = ep.calibrate.compute_Sv(ed,waveform_mode='BB',encode_mode='complex')
+    
+    # Extract mask for both real backscatter and calibrated Sv:
+    
+    # EchoData Real (component) Backscatter mask should be 1 if all values along 
+    # a specific range sample and beam dimension are non-NaN; else NaN.
+    ed_backscatter_r_mask = xr.where(
+        ~np.isnan(ed["Sonar/Beam_group1"]["backscatter_r"]
+    ).all(dim=['range_sample', "beam"]), 1, np.nan)
+    # Calibrated Sv mask should be 1 if all values along 
+    # a specific range sample are non-NaN; else NaN.
+    calibrated_Sv_mask = xr.where(
+        ~np.isnan(ds_Sv["Sv"]
+    ).all(dim=['range_sample']), 1, np.nan)
+
+    # Check that they are equal
+    assert np.array_equal(
+        ed_backscatter_r_mask.data,
+        calibrated_Sv_mask.data,
+        equal_nan=True
+    )
+
+    # Check that a slice of calibrated Sv mask is equal to the following array
+    # with the target channel ping pattern:
+    assert np.array_equal(
+        calibrated_Sv_mask.isel(ping_time=slice(0,10)).data,
+        target_channel_ping_pattern,
+        equal_nan=True
+    )

echopype/tests/calibrate/test_ek80_complex.py

@@ -70,4 +70,4 @@ def test_get_vend_filter_EK80(ch_num, filter_len, has_nan):
 
             assert sel_vend[var_df].values == get_vend_filter_EK80(
                 vend, channel_id=ch, filter_name=filter_name, param_type="decimation"
-            )
+            )

echopype/tests/calibrate/test_env_params.py

@@ -3,6 +3,7 @@
 import dask.array
 import numpy as np
 import xarray as xr
+import pandas as pd
 
 import echopype as ep
 from echopype.calibrate.env_params import (
@@ -132,6 +133,17 @@ 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
+    data = np.array([1, np.nan, np.nan])
+    time = pd.date_range(start='2024-01-01', periods=3, freq='D')
+    da = xr.DataArray(data, dims='time1', coords={'time1': time})
+
+    # Check that output da has only 1 value and no dimension
+    output_da = harmonize_env_param_time(da, None)
+    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"),