Read Sentinel V specific data packets

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

@@ -14,12 +14,9 @@ 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_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 - Done: {'firmware_ver': 51.41, 'inst_make': 'TRDI', 'inst_type': 'ADCP', 'rotate_vars': ['vel'], 'has_imu': 0, 'inst_model': 'Workhorse', 'beam_angle': 20, 'carrier_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
-root - INFO -   1721 ensembles will be produced.
-
 root - INFO - 17 ncells, not BB
 root - DEBUG -   pos 0mb/1mb
 
@@ -99,3 +96,6 @@ 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!
+root - INFO - n 4: 768 0300
+root - DEBUG -   pos: 935, pos_: 0, nbyte: 70, k: 0, byte_offset: -1

mhkit/dolfyn/io/base.py

@@ -298,7 +298,7 @@ def _create_dataset(data):
                         "data": data["data_vars"][key],
                     }
 
-                elif "b5" in tg:
+                elif "b5" in key:
                     ds_dict[key] = {
                         "dims": ("range_b5", "time_b5"),
                         "data": data["data_vars"][key],
@@ -324,7 +324,7 @@ def _create_dataset(data):
                     # "vel_b5" sometimes stored as (1, range_b5, time_b5)
                     ds_dict[key] = {
                         "dims": ("range_b5", "time_b5"),
-                        "data": data["data_vars"][key][0],
+                        "data": data["data_vars"][key].squeeze(),
                     }
                 elif "sl" in key:
                     ds_dict[key] = {
@@ -357,12 +357,12 @@ def _create_dataset(data):
     r_list = [r for r in ds.coords if "range" in r]
     for ky in r_list:
         ds[ky].attrs["units"] = "m"
-        ds[ky].attrs["long_name"] = "Profile Range"
+        ds[ky].attrs["long_name"] = "Profile " + ky.capitalize().replace("_", " ")
         ds[ky].attrs["description"] = "Distance to the center of each depth bin"
     time_list = [t for t in ds.coords if "time" in t]
     for ky in time_list:
         ds[ky].attrs["units"] = "seconds since 1970-01-01 00:00:00"
-        ds[ky].attrs["long_name"] = "Time"
+        ds[ky].attrs["long_name"] = ky.capitalize().replace("_", " ")
         ds[ky].attrs["standard_name"] = "time"
 
     # Set dataset metadata

mhkit/dolfyn/io/rdi.py

@@ -18,9 +18,10 @@ def read_rdi(
     filename,
     userdata=None,
     nens=None,
-    debug_level=-1,
+    debug=0,
     vmdas_search=False,
     winriver=False,
+    search_num=20000,
     **kwargs,
 ) -> xr.Dataset:
     """
@@ -32,11 +33,11 @@ def read_rdi(
       Filename of TRDI file to read.
     userdata : True, False, or string of userdata.json filename
       Whether to read the '<base-filename>.userdata.json' file. Default = True
-    nens : None, int or 2-element tuple (start, stop)
-      Number of pings or ensembles to read from the file.
+    nens : None, int
+      Number of pings or ensembles to read from the file, starting from 0.
       Default is None, read entire file
-    debug_level : int
-      Debug level [0 - 2]. Default = -1
+    debug : int
+      Debug level [0 - 3]. Default = 0
     vmdas_search : bool
       Search from the end of each ensemble for the VMDAS navigation
       block.  The byte offsets are sometimes incorrect. Default = False
@@ -50,7 +51,7 @@ def read_rdi(
       An xarray dataset from the binary instrument data
     """
     # Start debugger logging
-    if debug_level >= 0:
+    if debug > 0:
         for handler in logging.root.handlers[:]:
             logging.root.removeHandler(handler)
         filepath = Path(filename)
@@ -65,7 +66,11 @@ def read_rdi(
     # Reads into a dictionary of dictionaries using netcdf naming conventions
     # Should be easier to debug
     rdr = _RDIReader(
-        filename, debug_level=debug_level, vmdas_search=vmdas_search, winriver=winriver
+        filename,
+        debug=debug,
+        vmdas_search=vmdas_search,
+        winriver=winriver,
+        search_num=search_num,
     )
     datNB, datBB = rdr.load_data(nens=nens)
 
@@ -123,7 +128,7 @@ def read_rdi(
         )
 
     # Close handler
-    if debug_level >= 0:
+    if debug > 0:
         for handler in logging.root.handlers[:]:
             logging.root.removeHandler(handler)
             handler.close()
@@ -159,17 +164,16 @@ def _set_rdi_declination(dat, fname, inplace):
 
 
 class _RDIReader:
-    def __init__(
-        self, fname, navg=1, debug_level=-1, vmdas_search=False, winriver=False
-    ):
+    def __init__(self, fname, debug, vmdas_search, winriver, search_num):
         self.fname = base._abspath(fname)
         print("\nReading file {} ...".format(fname))
-        self._debug_level = debug_level
+        self._debug_level = debug
         self._vmdas_search = vmdas_search
         self._winrivprob = winriver
         self._vm_source = 0
         self._pos = 0
         self.progress = 0
+        self.search_num = search_num
         self._cfac32 = np.float32(180 / 2**31)  # signed 32 to float
         self._cfac16 = np.float32(180 / 2**15)  # unsigned16 to float
         self._fixoffset = 0
@@ -194,11 +198,10 @@ def __init__(
             logging.info("self._BB {}".format(self._BB))
             logging.info("self.cfgBB: {}".format(self.cfgBB))
         self.f.seek(self._pos, 0)
-        self.n_avg = navg
 
-        self.ensemble = lib._ensemble(self.n_avg, self.cfg["n_cells"])
+        self.ensemble = lib._ensemble(self.cfg["n_cells"])
         if self._BB:
-            self.ensembleBB = lib._ensemble(self.n_avg, self.cfgBB["n_cells"])
+            self.ensembleBB = lib._ensemble(self.cfgBB["n_cells"])
 
         self.vars_read = lib._variable_setlist(["time"])
         if self._BB:
@@ -243,7 +246,7 @@ def read_hdrseg(self):
         fd = self.f
         hdr = self.hdr
         hdr["nbyte"] = fd.read_i16(1)
-        spare = fd.read_ui8(1)
+        fd.seek(1, 1)
         ndat = fd.read_ui8(1)
         hdr["dat_offsets"] = fd.read_ui16(ndat)
         self._nbyte = 4 + ndat * 2
@@ -278,29 +281,22 @@ def check_for_double_buffer(self):
 
     def load_data(self, nens=None):
         """Main function run after reader class is initiated."""
+
         if nens is None:
-            # Attempt to overshoot WinRiver2 or *Pro filesize
-            if (self.cfg["coord_sys"] == "ship") or (
-                self.cfg["inst_model"]
-                in [
-                    "RiverPro",
-                    "StreamPro",
-                ]
-            ):
-                self._nens = int(self._filesize / self.hdr["nbyte"] / self.n_avg * 1.1)
-            else:
-                # Attempt to overshoot other instrument filesizes
-                self._nens = int(self._npings / self.n_avg)
+            # Overshoot file size to pre-allocate enough ensembles
+            self._nens = int(self._npings * 10)
         elif nens.__class__ is tuple or nens.__class__ is list:
             raise Exception("    `nens` must be a integer")
         else:
             self._nens = nens
-        if self._debug_level > -1:
-            logging.info("  taking data from pings 0 - %d" % self._nens)
-            logging.info("  %d ensembles will be produced.\n" % self._nens)
+
+        # Pre-allocate data
         self.init_data()
 
-        for iens in range(self._nens):
+        iens = 0
+        while True:
+            if iens == self._nens:
+                break
             if not self.read_buffer():
                 self.remove_end(iens)
                 break
@@ -323,6 +319,9 @@ def load_data(self, nens=None):
                 # reset flag after all variables run
                 self.n_cells_diff = 0
 
+                # b5 clock flag
+                b5 = True if "ping_offset_time_b5" in cfg else False
+
                 # Set clock
                 clock = en.rtc[:, :]
                 if clock[0, 0] < 100:
@@ -340,8 +339,15 @@ def load_data(self, nens=None):
                         )
                     )
                     dat["coords"]["time"][iens] = np.nan
+                    if b5:
+                        dat["coords"]["time_b5"][iens] = np.nan
                 else:
                     dat["coords"]["time"][iens] = np.median(dates)
+                    if b5:
+                        dat["coords"]["time_b5"][iens] = (
+                            np.median(dates) + cfg["ping_offset_time_b5"]
+                        )
+                iens += 1
 
         # Finalize dataset (runs through both NB and BB)
         for dat, cfg in zip(datl, cfgl):
@@ -402,10 +408,37 @@ def read_buffer(self):
             self.ensembleBB.k = -1  # so that k+=1 gives 0 on the first loop.
         self.print_progress()
         hdr = self.hdr
-        while self.ensemble.k < self.ensemble.n_avg - 1:
+        while self.ensemble.k < 0:
             if not self.search_buffer():
                 return False
             startpos = fd.tell() - 2
+
+            noBytesInEnsemble = fd.read_i16(1)
+            # go back to start of ensemble
+            fd.seek(-4, 1)
+            # pack the entire ensemble into a bytearray
+            bytesInEnsemble = bytearray(fd.read_ui8(noBytesInEnsemble))
+            # get checksum (2 bytes unsigned integer)
+            checksum = fd.read_ui16(1)
+            # calculate checksum and check
+            # if the checksum is wrong, back up 100 bytes and search for the next
+            # ensemble
+            if (sum(bytesInEnsemble) & 0xFFFF) != checksum:
+                logging.warning(
+                    "Ensemble starting at startpos {} has a checksum error".format(
+                        startpos
+                    )
+                )
+                logging.warning(
+                    "checksum calculated = %s, actual checksum = %s\n"
+                    % ((sum(bytesInEnsemble) & 0xFFFF), checksum)
+                )
+                fd.seek(-100, 1)
+                self.read_buffer()
+            else:
+                # go back to start of ensemble
+                fd.seek(-noBytesInEnsemble, 1)
+
             self.read_hdrseg()
             if self._debug_level > -1:
                 logging.info("Read Header", hdr)
@@ -464,7 +497,7 @@ def search_buffer(self):
         """
         Check to see if the next bytes indicate the beginning of a
         data block.  If not, search for the next data block, up to
-        _search_num times.
+        search_num times.
         """
         fd = self.f
         id = fd.read_ui8(2)
@@ -479,8 +512,11 @@ def search_buffer(self):
             logging.info("cfgid0: [{:x}, {:x}]".format(*cfgid))
         # If not [127, 127] or if the file ends in the next ensemble
         while (cfgid != [127, 127]) or self.check_eof():
+            if search_cnt == self.search_num:
+                logging.debug(f"Stopped searching at byte position {fd.tell()}")
+                return False
+            # Search for the next header or the end of the file
             if cfgid == [127, 121]:
-                # Search for the next header or the end of the file
                 skipbytes = fd.read_i16(1)
                 fd.seek(skipbytes - 2, 1)
                 id = fd.read_ui8(2)
@@ -497,7 +533,7 @@ def search_buffer(self):
                 cfgid[0] = cfgid[1]
                 cfgid[1] = nextbyte
 
-        if pos_7f79 and self._debug_level > -1:
+        if pos_7f79 and (self._debug_level > -1):
             logging.info("Skipped junk data: [{:x}, {:x}]".format(*[127, 121]))
 
         if search_cnt > 0:
@@ -560,6 +596,10 @@ def read_dat(self, id):
             0: (defs.read_fixed, [False]),
             # 0001 2nd profile fixed leader
             1: (defs.read_fixed, [True]),
+            # 000B Wave parameters
+            11: (defs.skip_Nbyte, [51]),
+            # 000C Wave parameters - sea and swell
+            12: (defs.skip_Nbyte, [44]),
             # 0010 Surface layer fixed leader (RiverPro & StreamPro)
             16: (defs.read_fixed_sl, []),
             # 0080 1st profile variable leader
@@ -612,13 +652,13 @@ def read_dat(self, id):
             1793: (defs.skip_Ncol, [4]),  # 0701 number of pings
             1794: (defs.skip_Ncol, [4]),  # 0702 sum of squared vel
             1795: (defs.skip_Ncol, [4]),  # 0703 sum of velocities
-            2560: (defs.skip_Ncol, []),  # 0A00 Beam 5 velocity
-            2816: (defs.skip_Ncol, []),  # 0B00 Beam 5 correlation
-            3072: (defs.skip_Ncol, []),  # 0C00 Beam 5 amplitude
-            3328: (defs.skip_Ncol, []),  # 0D00 Beam 5 pct_good
+            2560: (defs.read_vel_b5, []),  # 0A00 Beam 5 velocity
+            2816: (defs.read_corr_b5, []),  # 0B00 Beam 5 correlation
+            3072: (defs.read_amp_b5, []),  # 0C00 Beam 5 amplitude
+            3328: (defs.read_prcnt_gd_b5, []),  # 0D00 Beam 5 pct_good
             # Fixed attitude data format for Ocean Surveyor ADCPs
             3000: (defs.skip_Nbyte, [32]),
-            3841: (defs.skip_Nbyte, [38]),  # 0F01 Beam 5 leader
+            3841: (defs.read_vel_b5_leader, []),  # 0F01 Beam 5 leader
             8192: (defs.read_vmdas, []),  # 2000
             # 2013 Navigation parameter data
             8211: (defs.skip_Nbyte, [83]),
@@ -648,8 +688,16 @@ def read_dat(self, id):
             22785: (defs.skip_Nbyte, [65]),
             # 5902 Ping attitude
             22786: (defs.skip_Nbyte, [105]),
-            # 7001 ADC data
-            28673: (defs.skip_Nbyte, [14]),
+            # 7000 Sentinvel V system configuration
+            28672: (defs.read_sentinelv_syscfg, [False]),
+            # 7001 Sentinel V leader
+            28673: (defs.read_sentinelv_ping_setup, [False]),
+            # 7002 ADC data
+            28674: (defs.skip_Nbyte, [14]),
+            # 7003 Sentinel V "Features" (only first ensemble)
+            28675: (defs.skip_Nbyte, [88]),  # min size
+            # 7004 Sentinel V Event Log
+            28676: (defs.read_sentinelv_event_log, []),
         }
         # Call the correct function:
         if self._debug_level > 1:
@@ -787,10 +835,7 @@ def save_profiles(self, dat, cfg, nm, en, iens):
             The updated dataset dictionary with the reformatted profile measurements.
         """
         ds = lib._get(dat, nm)
-        if self.n_avg == 1:
-            bn = en[nm][..., 0]
-        else:
-            bn = np.nanmean(en[nm], axis=-1)
+        bn = en[nm][..., 0]
 
         # If n_cells has changed (RiverPro/StreamPro WinRiver transects)
         if len(ds.shape) == 3:
@@ -873,6 +918,12 @@ def cleanup(self, dat, cfg):
         ).astype(np.float32)
         cfg["range_offset"] = round(bin1_dist - cfg["blank_dist"] - cfg["cell_size"], 3)
 
+        if "n_cells_b5" in cfg:
+            bin1_dist_b5 = cfg.pop("bin1_dist_b5_m")
+            dat["coords"]["range_b5"] = (
+                bin1_dist_b5 + np.arange(cfg["n_cells_b5"]) * cfg["cell_size_b5"]
+            ).astype(np.float32)
+
         # Clean up surface layer profiles
         if "surface_layer" in cfg:  # RiverPro/StreamPro
             # Set SL cell size and range
@@ -986,20 +1037,22 @@ def finalize(self, dat, cfg):
         for nm in set(defs.data_defs.keys()) - self.vars_read:
             lib._pop(dat, nm)
 
-        # VMDAS and WinRiver have different set sampling frequency
-        if ("sourceprog" in cfg) and (
-            cfg["sourceprog"].lower() in ["vmdas", "winriver", "winriver2"]
-        ):
+        # Need to figure out how to differentiate burst mode from averaging mode
+        if (
+            ("source_program" in cfg)
+            and (cfg["source_program"].lower() in ["vmdas", "winriver", "winriver2"])
+        ) or ("sentinelv" in cfg["inst_model"].lower()):
             cfg["fs"] = round(1 / np.median(np.diff(dat["coords"]["time"])), 2)
         else:
             cfg["fs"] = 1 / (cfg["sec_between_ping_groups"] * cfg["pings_per_ensemble"])
 
         # Save configuration data as attributes
         dat["attrs"] = cfg
 
+        # Set 3D variable axes properly (beam, range, time)
         for nm in defs.data_defs:
             shp = defs.data_defs[nm][0]
-            if len(shp) and shp[0] == "nc" and lib._in_group(dat, nm):
+            if (len(shp) == 2) and (shp[0] == "nc") and lib._in_group(dat, nm):
                 lib._setd(dat, nm, np.swapaxes(lib._get(dat, nm), 0, 1))
 
         return dat