Added support for split-beam angle computation with mixed channel beam types in complex mode

Author: beatfactor

echopype/consolidate/split_beam_angle.py

@@ -8,9 +8,28 @@
 import numpy as np
 import xarray as xr
 
+from ..utils.log import _init_logger
 from ..calibrate.ek80_complex import compress_pulse, get_norm_fac, get_transmit_signal
 
 
+logger = _init_logger(__name__)
+
+# Beam type identifiers
+BEAM_TYPE_SPLIT_4_SECTOR = 1      # 4-sector split-beam (common Simrad type)
+BEAM_TYPE_SPLIT_3_SECTOR = 17     # 3-sector
+BEAM_TYPE_SPLIT_3_PLUS_CENTER = 49  # 3-sector + center element
+BEAM_TYPE_SPLIT_VARIANT_65 = 65   # Another 3+1 variant (vendor-specific)
+BEAM_TYPE_SPLIT_VARIANT_81 = 81   # Another 3+1 variant (vendor-specific)
+
+SUPPORTED_BEAM_TYPES = [
+    BEAM_TYPE_SPLIT_4_SECTOR,
+    BEAM_TYPE_SPLIT_3_SECTOR,
+    BEAM_TYPE_SPLIT_3_PLUS_CENTER,
+    BEAM_TYPE_SPLIT_VARIANT_65,
+    BEAM_TYPE_SPLIT_VARIANT_81,
+]
+
+
 def _compute_angle_from_complex(
     bs: xr.DataArray, beam_type: int, sens: List[xr.DataArray], offset: List[xr.DataArray]
 ) -> Tuple[xr.DataArray, xr.DataArray]:
@@ -48,7 +67,7 @@ def _compute_angle_from_complex(
     """
 
     # 4-sector transducer
-    if beam_type == 1:
+    if beam_type == BEAM_TYPE_SPLIT_4_SECTOR:
         bs_fore = (bs.isel(beam=2) + bs.isel(beam=3)) / 2  # forward
         bs_aft = (bs.isel(beam=0) + bs.isel(beam=1)) / 2  # aft
         bs_star = (bs.isel(beam=0) + bs.isel(beam=3)) / 2  # starboard
@@ -60,9 +79,14 @@ def _compute_angle_from_complex(
         phi = np.arctan2(np.imag(bs_phi), np.real(bs_phi)) / np.pi * 180
 
     # 3-sector transducer with or without center element
-    elif beam_type in [17, 49, 65, 81]:
+    elif beam_type in [
+        BEAM_TYPE_SPLIT_3_SECTOR,
+        BEAM_TYPE_SPLIT_3_PLUS_CENTER,
+        BEAM_TYPE_SPLIT_VARIANT_65,
+        BEAM_TYPE_SPLIT_VARIANT_81,
+    ]:
         # 3-sector
-        if beam_type == 17:
+        if beam_type == BEAM_TYPE_SPLIT_3_SECTOR:
             bs_star = bs.isel(beam=0)
             bs_port = bs.isel(beam=1)
             bs_fore = bs.isel(beam=2)
@@ -209,8 +233,14 @@ def get_angle_complex_samples(
         )
     else:
         # beam_type different for some channels, process each channel separately
-        theta, phi = [], []
+        theta_list, phi_list, valid_channels = [], [], []
         for ch_id in bs["channel"].data:
+            beam_type_ch = ds_beam["beam_type"].sel(channel=ch_id).item()
+
+            if beam_type_ch not in SUPPORTED_BEAM_TYPES:
+                logger.warning(f"Skipping channel {ch_id}: unsupported beam_type {beam_type_ch}")
+                continue
+
             theta_ch, phi_ch = _compute_angle_from_complex(
                 bs=bs.sel(channel=ch_id),
                 # beam_type is not time-varying
@@ -224,25 +254,18 @@ def get_angle_complex_samples(
                     angle_params["angle_offset_athwartship"].sel(channel=ch_id),
                 ],
             )
-            theta.append(theta_ch)
-            phi.append(phi_ch)
+            theta_list.append(theta_ch)
+            phi_list.append(phi_ch)
+            valid_channels.append(ch_id)
+
+        if not theta_list:
+            raise ValueError("No valid channels found for angle computation.")
 
         # Combine angles from all channels
-        theta = xr.DataArray(
-            data=theta,
-            coords={
-                "channel": bs["channel"],
-                "ping_time": bs["ping_time"],
-                "range_sample": bs["range_sample"],
-            },
-        )
-        phi = xr.DataArray(
-            data=phi,
-            coords={
-                "channel": bs["channel"],
-                "ping_time": bs["ping_time"],
-                "range_sample": bs["range_sample"],
-            },
-        )
+        theta = xr.concat(theta_list, dim="channel")
+        theta = theta.assign_coords(channel=("channel", valid_channels))
+
+        phi = xr.concat(phi_list, dim="channel")
+        phi = phi.assign_coords(channel=("channel", valid_channels))
 
     return theta, phi

echopype/tests/consolidate/test_consolidate_integration.py

@@ -335,5 +335,34 @@ def test_add_splitbeam_angle_BB_pc(test_path):
     assert np.allclose(pyel_vals, ep_vals, atol=1e-6)
 
 
+def test_add_splitbeam_angle_partial_valid_channels(test_path):
+    raw_file = test_path["EK80_CAL"] / "2018115-D20181213-T094600.raw"
+    ed = ep.open_raw(raw_file, sonar_model="EK80")
+
+    # Manually override beam_type for one channel to simulate single-beam
+    beam_types = ed["Sonar/Beam_group1"]["beam_type"].values
+    total_channels = beam_types.shape[0]
+
+    # Force the first channel to an unsupported beam_type
+    beam_types[1] = 0
+    ed["Sonar/Beam_group1"]["beam_type"].data[:] = beam_types
+
+    # Compute Sv
+    ds_Sv = ep.calibrate.compute_Sv(ed, waveform_mode="CW", encode_mode="complex")
+
+    # Add split-beam angles
+    ds_Sv = ep.consolidate.add_splitbeam_angle(source_Sv=ds_Sv, echodata=ed, waveform_mode="CW",
+                                               encode_mode="complex", pulse_compression=False, to_disk=False)
+
+    valid_angle_channels = [
+        ch for ch in ds_Sv.channel.values
+        if not np.all(np.isnan(ds_Sv["angle_alongship"].sel(channel=ch)))
+    ]
+
+    # Verify angles exist for valid channel only
+    assert "angle_alongship" in ds_Sv
+    assert "angle_athwartship" in ds_Sv
+    assert len(valid_angle_channels) == total_channels - 1
+
 # TODO: need a test for power/angle data, with mock EchoData object
 # containing some channels with single-beam data and some channels with split-beam data