Nortek Signature Reader Altimeter Update (#280)

* #277 
* #284 - Fix window check
* Test updates plus compression bugfix

examples/data/dolfyn/test_data/BenchFile01.repr.txt

@@ -4,6 +4,7 @@
   . (100 pings @ 2Hz)
   Variables:
   - time ('time',)
+  - time_altraw ('time_altraw',)
   - time_b5 ('time_b5',)
   - vel ('dir', 'range', 'time')
   - vel_b5 ('range_b5', 'time_b5')
@@ -14,6 +15,8 @@
   - roll ('time',)
   - temp ('time',)
   - pressure ('time',)
+  - pressure_alt ('time',)
+  - pressure_altraw ('time_altraw',)
   - amp ('beam', 'range', 'time')
   - amp_b5 ('range_b5', 'time_b5')
   - corr ('beam', 'range', 'time')

mhkit/dolfyn/io/api.py

@@ -76,7 +76,7 @@ def read(fname, userdata=True, nens=None, **kwargs):
     userdata : True, False, or string of userdata.json filename (default ``True``)
       Whether to read the '<base-filename>.userdata.json' file.
     nens : None, int or 2-element tuple (start, stop)
-      Number of pings or ensembles to read from the file. 
+      Number of pings or ensembles to read from the file.
       Default is None, read entire file
     **kwargs : dict
       Passed to instrument-specific parser.
@@ -173,8 +173,8 @@ def save(ds, filename,
     ds.attrs['complex_vars'] = []
     for var in ds.data_vars:
         if np.iscomplexobj(ds[var]):
-            ds[var+'_real'] = ds[var].real
-            ds[var+'_imag'] = ds[var].imag
+            ds[var + '_real'] = ds[var].real
+            ds[var + '_imag'] = ds[var].imag
 
             ds = ds.drop_vars(var)
             ds.attrs['complex_vars'].append(var)
@@ -183,15 +183,25 @@ def save(ds, filename,
         elif ds[var].dtype == np.float64:
             ds[var] = ds[var].astype('float32')
 
-    if compression:
-        enc = dict()
-        for ky in ds.variables:
-            enc[ky] = dict(zlib=True, complevel=1)
-        if 'encoding' in kwargs:
-            # Overwrite ('update') values in enc with whatever is in kwargs['encoding']
-            enc.update(kwargs['encoding'])
-        else:
-            kwargs['encoding'] = enc
+    # Write variable encoding
+    enc = dict()
+    if 'encoding' in kwargs:
+        enc.update(kwargs['encoding'])
+    for ky in ds.variables:
+        # Save prior encoding
+        enc[ky] = ds[ky].encoding
+        # Remove unexpected netCDF4 encoding parameters
+        # https://github.com/pydata/xarray/discussions/5709
+        params = ['szip', 'zstd', 'bzip2', 'blosc', 'contiguous', 'chunksizes']
+        [enc[ky].pop(p) for p in params if p in enc[ky]]
+
+        if compression:
+            # New netcdf4-c cannot compress variable length strings
+            if isinstance(ds[ky].data[0], str):
+                continue
+            enc[ky].update(dict(zlib=True, complevel=1))
+
+    kwargs['encoding'] = enc
 
     # Fix encoding on datetime64 variables.
     ds = _decode_cf(ds)
@@ -220,19 +230,20 @@ def load(filename):
 
     # Convert numpy arrays and strings back to lists
     for nm in ds.attrs:
-        if type(ds.attrs[nm]) == np.ndarray and ds.attrs[nm].size > 1:
+        if isinstance(ds.attrs[nm], np.ndarray) and ds.attrs[nm].size > 1:
             ds.attrs[nm] = list(ds.attrs[nm])
-        elif type(ds.attrs[nm]) == str and nm in ['rotate_vars']:
+        elif isinstance(ds.attrs[nm], str) and nm in ['rotate_vars']:
             ds.attrs[nm] = [ds.attrs[nm]]
 
     # Rejoin complex numbers
-    if hasattr(ds, 'complex_vars') and len(ds.complex_vars):
-        if len(ds.complex_vars[0]) == 1:
-            ds.attrs['complex_vars'] = [ds.complex_vars]
-        for var in ds.complex_vars:
-            ds[var] = ds[var+'_real'] + ds[var+'_imag'] * 1j
-            ds = ds.drop_vars([var+'_real', var+'_imag'])
-    ds.attrs.pop('complex_vars')
+    if hasattr(ds, 'complex_vars'):
+        if len(ds.complex_vars):
+            if len(ds.complex_vars[0]) == 1:
+                ds.attrs['complex_vars'] = [ds.complex_vars]
+            for var in ds.complex_vars:
+                ds[var] = ds[var + '_real'] + ds[var + '_imag'] * 1j
+                ds = ds.drop_vars([var + '_real', var + '_imag'])
+        ds.attrs.pop('complex_vars')
 
     return ds
 
@@ -318,7 +329,7 @@ def load_mat(filename, datenum=True):
     filename : str
       Filename and/or path with the '.mat' extension
     datenum : bool
-      If true, converts time from datenum. If false, converts time from 
+      If true, converts time from datenum. If false, converts time from
       "epoch time".
 
     Returns
@@ -351,12 +362,12 @@ def load_mat(filename, datenum=True):
 
     # Convert numpy arrays and strings back to lists
     for nm in ds.attrs:
-        if type(ds.attrs[nm]) == np.ndarray and ds.attrs[nm].size > 1:
+        if isinstance(ds.attrs[nm], np.ndarray) and ds.attrs[nm].size > 1:
             try:
                 ds.attrs[nm] = [x.strip(' ') for x in list(ds.attrs[nm])]
             except:
                 ds.attrs[nm] = list(ds.attrs[nm])
-        elif type(ds.attrs[nm]) == str and nm in ['time_coords', 'time_data_vars', 'rotate_vars']:
+        elif isinstance(ds.attrs[nm], str) and nm in ['time_coords', 'time_data_vars', 'rotate_vars']:
             ds.attrs[nm] = [ds.attrs[nm]]
 
     if hasattr(ds, 'orientation_down'):

mhkit/dolfyn/io/base.py

@@ -118,7 +118,7 @@ def _create_dataset(data):
     Direction 'dir' coordinates are set in `set_coords`
     """
     ds = xr.Dataset()
-    tag = ['_avg', '_b5', '_echo', '_bt', '_gps', '_ast', '_sl']
+    tag = ['_avg', '_b5', '_echo', '_bt', '_gps', '_altraw', '_sl']
 
     FoR = {}
     try:
@@ -205,6 +205,11 @@ def _create_dataset(data):
                                               'dim_1': 'time_echo'})
                     ds[key] = ds[key].assign_coords({'range_echo': data['coords']['range_echo'],
                                                      'time_echo': data['coords']['time_echo']})
+                elif key == 'samp_altraw':  # raw altimeter samples
+                    ds[key] = ds[key].rename({'dim_0': 'n_altraw',
+                                              'dim_1': 'time_altraw'})
+                    ds[key] = ds[key].assign_coords({'time_altraw': data['coords']['time_altraw']})
+
                 # ADV/ADCP instrument vector data, bottom tracking
                 elif shp[0] == n_beams and not any(val in key for val in tag[:3]):
                     if 'bt' in key and 'time_bt' in data['coords']:

mhkit/dolfyn/io/nortek.py

@@ -165,8 +165,11 @@ class _NortekReader():
                '0x10': 'read_vec_data',
                '0x11': 'read_vec_sysdata',
                '0x12': 'read_vec_hdr',
-               '0x71': 'read_microstrain',
                '0x20': 'read_awac_profile',
+               '0x30': 'read_awac_waves',
+               '0x31': 'read_awac_waves_hdr',
+               '0x36': 'read_awac_waves',  # "SUV"
+               '0x71': 'read_microstrain',
                }
 
     def __init__(self, fname, endian=None, debug=False,
@@ -1026,6 +1029,77 @@ def sci_awac_profile(self,):
         self.data['attrs']['cell_size'] = cs
         self.data['attrs']['blank_dist'] = bd
 
+    def read_awac_waves_hdr(self,):
+        # ID: '0x31'
+        c = self.c
+        if self.debug:
+            print('Reading vector header data (0x31) ping #{} @ {}...'
+                  .format(self.c, self.pos))
+        hdrnow = {}
+        dat = self.data
+        ds = dat['sys']
+        dv = dat['data_vars']
+        if 'time' not in dat['coords']:
+            self._init_data(nortek_defs.waves_hdrdata)
+        byts = self.read(56)
+        # The first two are size, the next 6 are time.
+        tmp = unpack(self.endian + '8x4H3h2HhH4B6H5h', byts)
+        dat['coords']['time'][c] = self.rd_time(byts[2:8])
+        hdrnow['n_records_alt'] = tmp[0]
+        hdrnow['blank_dist_alt'] = tmp[1]  # counts
+        ds['batt_alt'][c] = tmp[2]  # voltage (0.1 V)
+        dv['c_sound_alt'][c] = tmp[3]  # c (0.1 m/s)
+        dv['heading_alt'][c] = tmp[4]  # (0.1 deg)
+        dv['pitch_alt'][c] = tmp[5]  # (0.1 deg)
+        dv['roll_alt'][c] = tmp[6]  # (0.1 deg)
+        dv['pressure1_alt'][c] = tmp[7]  # min pressure previous profile (0.001 dbar)
+        dv['pressure2_alt'][c] = tmp[8]  # max pressure previous profile (0.001 dbar)
+        dv['temp_alt'][c] = tmp[9]  # (0.01 deg C)
+        hdrnow['cell_size_alt'][c] = tmp[10]  # (counts of T3)
+        hdrnow['noise_alt'][c] = tmp[11:15]  # noise amplitude beam 1-4 (counts)
+        hdrnow['proc_magn_alt'][c] = tmp[15:19]  # processing magnitude beam 1-4
+        hdrnow['n_past_window_alt'] = tmp[19]  # number of samples of AST window past boundary
+        hdrnow['n_window_alt'] = tmp[20]  # AST window size (# samples)
+        hdrnow['Spare1'] = tmp[21:]
+        self.checksum(byts)
+        if 'data_header' not in self.config:
+            self.config['data_header'] = hdrnow
+        else:
+            if not isinstance(self.config['data_header'], list):
+                self.config['data_header'] = [self.config['data_header']]
+            self.config['data_header'] += [hdrnow]
+
+    def read_awac_waves(self,):
+        """Read awac wave and suv data
+        """
+        # IDs: 0x30 & 0x36
+        c = self.c
+        dat = self.data
+        if self.debug:
+            print('Reading awac wave data (0x30) ping #{} @ {}...'
+                  .format(self.c, self.pos))
+        if 'dist1_alt' not in dat['data_vars']:
+            self._init_data(nortek_defs.wave_data)
+            self._dtypes += ['wave_data']
+        # The first two are size
+        byts = self.read(20)
+        ds = dat['sys']
+        dv = dat['data_vars']
+        (dv['pressure'][c],  # (0.001 dbar)
+         dv['dist1_alt'][c],  # distance 1 to surface, vertical beam (mm)
+         ds['AnaIn_alt'][c],  # analog input 1
+         dv['vel_alt'][0, c],  # velocity beam 1 (mm/s) East for SUV
+         dv['vel_alt'][1, c],  # North for SUV
+         dv['vel_alt'][2, c],  # Up for SUV
+         dv['dist2_alt'][c],  # distance 2 to surface, vertical beam (mm) or vel 4 for non-AST
+         dv['amp_alt'][0, c],  # amplitude beam 1 (counts)
+         dv['amp_alt'][1, c],  # amplitude beam 2 (counts)
+         dv['amp_alt'][2, c],  # amplitude beam 3 (counts)
+         # AST quality (counts) or amplitude beam 4 for non-AST
+         dv['quality_alt'][c]) = unpack(self.endian + '3H4h4B', byts)
+        self.checksum(byts)
+        self.c += 1
+
     def dat2sci(self,):
         for nm in self._dtypes:
             getattr(self, 'sci_' + nm)()