Fix ever-changing RDI RiverPro depth bin ranges

examples/data/dolfyn/test_data/RDI_withBT.dolfyn.log

@@ -3,8 +3,9 @@ root - INFO - cfgid0: [7f, 7f]
 root - INFO - {'nbyte': 579, 'dat_offsets': array([ 20,  79, 144, 282, 352, 422, 492])}
 root - INFO - pos 20
 root - INFO - id 0 offset 20
-root - INFO - Number of cells set to 17
-root - INFO - Cell size set to 1.0
+root - DEBUG - Number of cells set to 17
+root - DEBUG - Cell size set to 1.0
+root - DEBUG - Bin 1 distance set to 2.09
 root - INFO - Read Config
 root - INFO - Read Fixed
 root - INFO - id 128 offset 79
@@ -13,7 +14,7 @@ root - INFO - id 512 offset 282
 root - INFO - id 768 offset 352
 root - INFO - id 1024 offset 422
 root - INFO - id 1536 offset 492
-root - INFO - Done: {'prog_ver': 51.41, 'inst_model': 'Workhorse', 'beam_angle': 20, 'freq': 600, 'beam_pattern': 'convex', 'orientation': 'down', 'n_beams': 4, 'n_cells': 17, 'pings_per_ensemble': 1, 'cell_size': 1.0, 'blank_dist': 0.88, 'profiling_mode': 1, 'min_corr_threshold': 64, 'n_code_reps': 5, 'min_prcnt_gd': 0, 'max_error_vel': 2.0, 'sec_between_ping_groups': 0.5, 'coord_sys': 'earth', 'use_pitchroll': 'yes', 'use_3beam': 'yes', 'bin_mapping': 'yes', 'heading_misalign_deg': 0.0, 'magnetic_var_deg': 0.0, 'sensors_src': '01111101', 'sensors_avail': '00111101', 'bin1_dist_m': 2.09, 'transmit_pulse_m': 1.18, 'water_ref_cells': [1, 5], 'false_target_threshold': 50, 'transmit_lag_m': 0.24, 'bandwidth': 0, 'power_level': 255, 'serialnum': 18655}
+root - INFO - Done: {'prog_ver': 51.41, 'inst_make': 'TRDI', 'inst_type': 'ADCP', 'rotate_vars': ['vel'], 'has_imu': 0, 'inst_model': 'Workhorse', 'beam_angle': 20, 'freq': 600, 'beam_pattern': 'convex', 'orientation': 'down', 'n_beams': 4, 'n_cells': 17, 'pings_per_ensemble': 1, 'cell_size': 1.0, 'blank_dist': 0.88, 'profiling_mode': 1, 'min_corr_threshold': 64, 'n_code_reps': 5, 'min_prcnt_gd': 0, 'max_error_vel': 2.0, 'sec_between_ping_groups': 0.5, 'coord_sys': 'earth', 'use_pitchroll': 'yes', 'use_3beam': 'yes', 'bin_mapping': 'yes', 'heading_misalign_deg': 0.0, 'magnetic_var_deg': 0.0, 'sensors_src': '01111101', 'sensors_avail': '00111101', 'bin1_dist_m': 2.09, 'transmit_pulse_m': 1.18, 'water_ref_cells': [1, 5], 'false_target_threshold': 50, 'transmit_lag_m': 0.24, 'bandwidth': 0, 'power_level': 255, 'serialnum': 18655}
 root - INFO - self._bb False
 root - INFO - self.cfgbb: {}
 root - INFO -   taking data from pings 0 - 1721
@@ -98,4 +99,3 @@ root - DEBUG -   pos: 865, pos_: 0, nbyte: 138, k: 0, byte_offset: -1
 root - DEBUG - Trying to Read 512
 root - INFO -   Reading code 0x200...
 root - INFO - Read Corr
-root - INFO -     success!

mhkit/acoustics/__init__.py

@@ -1,17 +1,17 @@
 """
 The passive acoustics module provides a set of functions
-for analyzing and visualizing passive acoustic monitoring 
+for analyzing and visualizing passive acoustic monitoring
 data deployed in water bodies. This package reads in raw
-*.wav* files and conducts basic acoustics analysis and 
+*.wav* files and conducts basic acoustics analysis and
 visualization.
 
 To start using the module, import it directly from MHKiT:
 ``from mhkit import acoustics``. The analysis functions
-are available directly from the main import, while the 
-I/O and graphics submodules are available from 
+are available directly from the main import, while the
+I/O and graphics submodules are available from
 ``acoustics.io`` and  ``acoustics.graphics``, respectively.
-The base functions are intended to be used on top of the I/O submodule, and 
-include functionality to calibrate data, create spectral densities, sound 
+The base functions are intended to be used on top of the I/O submodule, and
+include functionality to calibrate data, create spectral densities, sound
 pressure levels, and time or band aggregate spectral data.
 """
 

mhkit/acoustics/analysis.py

@@ -337,7 +337,7 @@ def sound_pressure_spectral_density_level(spsd: xr.DataArray) -> xr.DataArray:
 
 
 def _validate_method(
-    method: Union[str, Dict[str, Union[float, int]]]
+    method: Union[str, Dict[str, Union[float, int]]],
 ) -> Tuple[str, Optional[Union[float, int]]]:
     """
     Validates the 'method' parameter and returns the method name and its argument (if any)

mhkit/acoustics/graphics.py

@@ -1,6 +1,6 @@
 """
-This submodule provides essential plotting functions for visualizing passive acoustics 
-data. The functions allow for customizable plotting of sound pressure spectral density 
+This submodule provides essential plotting functions for visualizing passive acoustics
+data. The functions allow for customizable plotting of sound pressure spectral density
 levels across time and frequency dimensions.
 
 Each plotting function leverages the flexibility of Matplotlib, allowing for passthrough
@@ -11,12 +11,12 @@
 -------------
 1. **plot_spectrogram**:
 
-   - Generates a spectrogram plot from sound pressure spectral density level data, 
+   - Generates a spectrogram plot from sound pressure spectral density level data,
      with a logarithmic frequency scale by default for improved readability of acoustic data.
 
 2. **plot_spectra**:
 
-   - Produces a spectral density plot with a log-transformed x-axis, allowing for clear 
+   - Produces a spectral density plot with a log-transformed x-axis, allowing for clear
      visualization of spectral density across frequency bands.
 """
 

mhkit/acoustics/io.py

@@ -1,7 +1,7 @@
 """
 This submodule provides input/output functions for passive acoustics data,
 focusing on hydrophone recordings stored in WAV files. The main functionality
-includes reading and processing hydrophone data from various manufacturers 
+includes reading and processing hydrophone data from various manufacturers
 and exporting audio files for easy playback and analysis.
 
 Supported Hydrophone Models
@@ -14,28 +14,28 @@
 
 1. **Data Reading**:
 
-   - `read_hydrophone`: Main function to read a WAV file from a hydrophone and 
-     convert it to either a voltage or pressure time series, depending on the 
+   - `read_hydrophone`: Main function to read a WAV file from a hydrophone and
+     convert it to either a voltage or pressure time series, depending on the
      availability of sensitivity data.
 
-   - `read_soundtrap`: Wrapper for reading Ocean Instruments SoundTrap hydrophone 
+   - `read_soundtrap`: Wrapper for reading Ocean Instruments SoundTrap hydrophone
      files, automatically using appropriate metadata.
 
-   - `read_iclisten`: Wrapper for reading Ocean Sonics icListen hydrophone files, 
-     including metadata processing to apply hydrophone sensitivity for direct 
+   - `read_iclisten`: Wrapper for reading Ocean Sonics icListen hydrophone files,
+     including metadata processing to apply hydrophone sensitivity for direct
      sound pressure calculation.
 
 2. **Audio Export**:
 
-   - `export_audio`: Converts processed sound pressure data back into a WAV file 
+   - `export_audio`: Converts processed sound pressure data back into a WAV file
      format, with optional gain adjustment to improve playback quality.
 
 3. **Data Extraction**:
 
-   - `_read_wav_metadata`: Extracts metadata from a WAV file, including bit depth 
+   - `_read_wav_metadata`: Extracts metadata from a WAV file, including bit depth
      and other header information.
 
-   - `_calculate_voltage_and_time`: Converts raw WAV data into voltage values and 
+   - `_calculate_voltage_and_time`: Converts raw WAV data into voltage values and
      generates a time index based on the sampling frequency.
 """
 

mhkit/dolfyn/adv/clean.py

@@ -1,5 +1,4 @@
-"""Module containing functions to clean data
-"""
+"""Module containing functions to clean data"""
 
 import warnings
 import numpy as np

mhkit/dolfyn/io/rdi.py

@@ -177,7 +177,9 @@ def __init__(
         self.n_cells_diff = 0
         self.n_cells_sl = 0
         self.cs_diff = 0
+        self.cs_sl_diff = 0
         self.cs = []
+        self.cs_sl = []
         self.cfg = {}
         self.cfgbb = {}
         self.hdr = {}
@@ -344,47 +346,32 @@ def load_data(self, nens=None):
         # Finalize dataset (runs through both nb and bb)
         for dat, cfg in zip(datl, cfgl):
             dat, cfg = self.cleanup(dat, cfg)
-            dat = self.finalize(dat)
+            dat = self.finalize(dat, cfg)
             if "vel_bt" in dat["data_vars"]:
-                dat["attrs"]["rotate_vars"].append("vel_bt")
+                cfg["rotate_vars"].append("vel_bt")
 
         datbb = self.outdBB if self._bb else None
-        return self.outd, datbb
+        return dat, datbb
 
     def init_data(self):
         """Initiate data structure"""
         outd = {
             "data_vars": {},
             "coords": {},
-            "attrs": {},
             "units": {},
             "long_name": {},
             "standard_name": {},
             "sys": {},
         }
-        outd["attrs"]["inst_make"] = "TRDI"
-        outd["attrs"]["inst_type"] = "ADCP"
-        outd["attrs"]["rotate_vars"] = [
-            "vel",
-        ]
-        # Currently RDI doesn't use IMUs
-        outd["attrs"]["has_imu"] = 0
         if self._bb:
             outdbb = {
                 "data_vars": {},
                 "coords": {},
-                "attrs": {},
                 "units": {},
                 "long_name": {},
                 "standard_name": {},
                 "sys": {},
             }
-            outdbb["attrs"]["inst_make"] = "TRDI"
-            outdbb["attrs"]["inst_type"] = "ADCP"
-            outdbb["attrs"]["rotate_vars"] = [
-                "vel",
-            ]
-            outdbb["attrs"]["has_imu"] = 0
 
         # Preallocate variables and data sizes
         for nm in defs.data_defs:
@@ -828,6 +815,9 @@ def save_profiles(self, dat, nm, en, iens):
         if self.cs_diff:
             self.cs.append([iens, self.cfg["cell_size"]])
             self.cs_diff = 0
+        if self.cs_sl_diff:
+            self.cs_sl.append([iens, self.cfg["cell_size_sl"]])
+            self.cs_sl_diff = 0
 
         # Then copy the ensemble to the dataset.
         ds[..., iens] = bn
@@ -859,53 +849,69 @@ def cleanup(self, dat, cfg):
         """
         # Clean up changing cell size, if necessary
         cs = np.array(self.cs, dtype=np.float32)
-        cell_sizes = cs[:, 1]
+        cs_sl = np.array(self.cs_sl, dtype=np.float32)
 
         # If cell sizes change, depth-bin average the smaller cell sizes
         if len(self.cs) > 1:
-            bins_to_merge = cell_sizes.max() / cell_sizes
-            idx_start = cs[:, 0].astype(int)
-            idx_end = np.append(cs[1:, 0], self._nens).astype(int)
-
             dv = dat["data_vars"]
-            for var in dv:
-                if (len(dv[var].shape) == 3) and ("_sl" not in var):
-                    # Create a new NaN var to save data in
-                    new_var = (np.zeros(dv[var].shape) * np.nan).astype(dv[var].dtype)
-                    # For each cell size change, reshape and bin-average
-                    for id1, id2, b in zip(idx_start, idx_end, bins_to_merge):
-                        array = np.transpose(dv[var][..., id1:id2])
-                        bin_arr = np.transpose(np.mean(self.reshape(array, b), axis=-1))
-                        new_var[: len(bin_arr), :, id1:id2] = bin_arr
-                    # Reset data. This often leaves nan data at farther ranges
-                    dv[var] = new_var
+            self.merge_bins(cs, dv, sl=False)
+        if len(self.cs_sl) > 1:
+            dv = dat["data_vars"]
+            self.merge_bins(cs_sl, dv, sl=True)
 
         # Set cell size and range
         cfg["n_cells"] = self.ensemble["n_cells"]
-        cfg["cell_size"] = round(cell_sizes.max(), 3)
+        cfg["cell_size"] = round(cs[:, 1].max(), 3)
+        bin1_dist = cfg.pop("bin1_dist_m")
         dat["coords"]["range"] = (
-            cfg["bin1_dist_m"] + np.arange(cfg["n_cells"]) * cfg["cell_size"]
+            bin1_dist + np.arange(cfg["n_cells"]) * cfg["cell_size"]
         ).astype(np.float32)
-
-        # Save configuration data as attributes
-        for nm in cfg:
-            dat["attrs"][nm] = cfg[nm]
+        cfg["range_offset"] = round(bin1_dist - cfg["blank_dist"] - cfg["cell_size"], 3)
 
         # Clean up surface layer profiles
         if "surface_layer" in cfg:  # RiverPro/StreamPro
+            # Set SL cell size and range
+            cfg["cell_size_sl"] = round(cs_sl[:, 1].max(), 3)
+            cfg["n_cells_sl"] = self.n_cells_sl
+            bin1_dist_sl = cfg.pop("bin1_dist_m_sl")
+            # Blank distance not recorded
+            cfg["blank_dist_sl"] = round(bin1_dist_sl - cfg["cell_size_sl"], 3)
+            # Range offset not added in "bin1_dist_m_sl" for some reason
+            bin1_dist_sl += cfg["range_offset"]
             dat["coords"]["range_sl"] = (
-                cfg["bin1_dist_m_sl"]
-                + np.arange(0, self.n_cells_sl) * cfg["cell_size_sl"]
+                bin1_dist_sl + np.arange(0, self.n_cells_sl) * cfg["cell_size_sl"]
             )
             # Trim off extra nan data
             dv = dat["data_vars"]
             for var in dv:
                 if "sl" in var:
                     dv[var] = dv[var][: self.n_cells_sl]
-            dat["attrs"]["rotate_vars"].append("vel_sl")
+            cfg["rotate_vars"].append("vel_sl")
 
         return dat, cfg
 
+    def merge_bins(self, cs, dv, sl=False):
+        cell_sizes = cs[:, 1]
+        bins_to_merge = cell_sizes.max() / cell_sizes
+        idx_start = cs[:, 0].astype(int)
+        idx_end = np.append(cs[1:, 0], self._nens).astype(int)
+
+        for var in dv:
+            if not sl:
+                flag = "_sl" not in var
+            elif sl:
+                flag = "_sl" in var
+            if (len(dv[var].shape) == 3) and flag:
+                # Create a new NaN var to save data in
+                new_var = (np.zeros(dv[var].shape) * np.nan).astype(dv[var].dtype)
+                # For each cell size change, reshape and bin-average
+                for id1, id2, b in zip(idx_start, idx_end, bins_to_merge):
+                    array = np.transpose(dv[var][..., id1:id2])
+                    bin_arr = np.transpose(np.mean(self.reshape(array, b), axis=-1))
+                    new_var[: len(bin_arr), :, id1:id2] = bin_arr
+                # Reset data. This often leaves nan data at farther ranges
+                dv[var] = new_var
+
     def reshape(self, arr, n_bin=None):
         """
         Reshapes the input array `arr` to a shape of (..., n, n_bin).
@@ -949,7 +955,7 @@ def reshape(self, arr, n_bin=None):
 
         return out
 
-    def finalize(self, dat):
+    def finalize(self, dat, cfg):
         """
         This method cleans up the dataset by removing any attributes that were
         defined but not loaded, updates configuration attributes, and sets the
@@ -968,19 +974,21 @@ def finalize(self, dat):
         dict
             The finalized dataset dictionary with cleaned attributes and added metadata.
         """
+
+        # Drop empty data variables
         for nm in set(defs.data_defs.keys()) - self.vars_read:
             lib._pop(dat, nm)
-        for nm in self.cfg:
-            dat["attrs"][nm] = self.cfg[nm]
 
         # VMDAS and WinRiver have different set sampling frequency
-        da = dat["attrs"]
-        if ("sourceprog" in da) and (
-            da["sourceprog"].lower() in ["vmdas", "winriver", "winriver2"]
+        if ("sourceprog" in cfg) and (
+            cfg["sourceprog"].lower() in ["vmdas", "winriver", "winriver2"]
         ):
-            da["fs"] = round(1 / np.median(np.diff(dat["coords"]["time"])), 2)
+            cfg["fs"] = round(1 / np.median(np.diff(dat["coords"]["time"])), 2)
         else:
-            da["fs"] = 1 / (da["sec_between_ping_groups"] * da["pings_per_ensemble"])
+            cfg["fs"] = 1 / (cfg["sec_between_ping_groups"] * cfg["pings_per_ensemble"])
+
+        # Save configuration data as attributes
+        dat["attrs"] = cfg
 
         for nm in defs.data_defs:
             shp = defs.data_defs[nm][0]