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

@@ -205,11 +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
+                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

@@ -1046,20 +1046,20 @@ def read_awac_waves_hdr(self,):
         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['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:
@@ -1068,7 +1068,7 @@ def read_awac_waves_hdr(self,):
             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
         """
@@ -1085,16 +1085,16 @@ def read_awac_waves(self,):
         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)
+        (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)

mhkit/dolfyn/io/nortek2.py

@@ -26,7 +26,7 @@ def read_signature(filename, userdata=True, nens=None, rebuild_index=False,
     userdata : bool
       To search for and use a .userdata.json or not
     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
     rebuild_index : bool
       Force rebuild of dolfyn-written datafile index. Useful for code updates.
@@ -218,14 +218,14 @@ def _init_burst_readers(self, ):
     def init_data(self, ens_start, ens_stop):
         outdat = {}
         nens = int(ens_stop - ens_start)
-        
+
         # ID 26 usually only recorded in first ensemble
         n26 = ((self._index['ID'] == 26) &
                (self._index['ens'] >= ens_start) &
                (self._index['ens'] < ens_stop)).sum()
         if not n26 and 26 in self._burst_readers:
             self._burst_readers.pop(26)
-        
+
         for ky in self._burst_readers:
             if ky == 26:
                 n = n26
@@ -239,7 +239,7 @@ def init_data(self, ens_start, ens_stop):
             outdat[ky]['units'] = self._burst_readers[ky].data_units()
             outdat[ky]['long_name'] = self._burst_readers[ky].data_longnames()
             outdat[ky]['standard_name'] = self._burst_readers[ky].data_stdnames()
-        
+
         return outdat
 
     def _read_hdr(self, do_cs=False):
@@ -279,11 +279,11 @@ def readfile(self, ens_start=0, ens_stop=None):
             except IOError:
                 return outdat
             id = hdr['id']
-            if id in [21, 22, 23, 24, 28]:  # "burst data record" (vel + ast), 
+            if id in [21, 22, 23, 24, 28]:  # "burst data record" (vel + ast),
                 # "avg data record" (vel_avg + ast_avg), "bottom track data record" (bt),
                 # "interleaved burst data record" (vel_b5), "echosounder record" (echo)
                 self._read_burst(id, outdat[id], c)
-            elif id in [26]:  
+            elif id in [26]:
                 # "burst altimeter raw record" (alt_raw) - recorded on nens==0
                 rdr = self._burst_readers[26]
                 if not hasattr(rdr, '_nsamp_index'):
@@ -323,7 +323,7 @@ def readfile(self, ens_start=0, ens_stop=None):
                 outdat[id]['ensemble'][c26] = c
                 c26 += 1
 
-            elif id in [27, 29, 30, 31, 35, 36]: # unknown how to handle
+            elif id in [27, 29, 30, 31, 35, 36]:  # unknown how to handle
                 # "bottom track record", DVL, "altimeter record", "avg altimeter raw record", 
                 # "raw echosounder data record", "raw echosounder transmit data record"
                 if self.debug:
@@ -411,7 +411,7 @@ def _reorg(dat):
     cfg['inst_make'] = 'Nortek'
     cfg['inst_type'] = 'ADCP'
 
-    for id, tag in [(21, ''), (22, '_avg'), (23, '_bt'), 
+    for id, tag in [(21, ''), (22, '_avg'), (23, '_bt'),
                     (24, '_b5'), (26, 'raw'), (28, '_echo')]:
         if id in [24, 26]:
             collapse_exclude = [0]
@@ -478,9 +478,9 @@ def _reorg(dat):
     if 26 in dat:
         for ky in list(outdat['data_vars']):
             if ky.endswith('raw') and not ky.endswith('_altraw'):
-                 outdat['data_vars'].pop(ky)
+                outdat['data_vars'].pop(ky)
         outdat['coords']['time_altraw'] = outdat['coords'].pop('timeraw')
-        outdat['data_vars']['samp_altraw'] =  outdat['data_vars']['samp_altraw'].astype('float32') / 2**8  # convert "signed fractional" to float
+        outdat['data_vars']['samp_altraw'] = outdat['data_vars']['samp_altraw'].astype('float32') / 2**8  # convert "signed fractional" to float
 
         # Read altimeter status
         outdat['data_vars'].pop('status_altraw')
@@ -564,7 +564,7 @@ def _reduce(data):
     --- from different data structures within the same ensemble --- by
     averaging.
     """
-    
+
     dv = data['data_vars']
     dc = data['coords']
     da = data['attrs']
@@ -580,14 +580,14 @@ def _reduce(data):
 
     if 'vel' in dv:
         dc['range'] = ((np.arange(dv['vel'].shape[1])+1) *
-                    da['cell_size'] +
-                    da['blank_dist'])
+                       da['cell_size'] +
+                       da['blank_dist'])
         da['fs'] = da['filehead_config']['BURST']['SR']
         tmat = da['filehead_config']['XFBURST']
     if 'vel_avg' in dv:
         dc['range_avg'] = ((np.arange(dv['vel_avg'].shape[1])+1) *
-                    da['cell_size_avg'] +
-                    da['blank_dist_avg'])
+                           da['cell_size_avg'] +
+                           da['blank_dist_avg'])
         dv['orientmat'] = dv.pop('orientmat_avg')
         tmat = da['filehead_config']['XFAVG']
         da['fs'] = da['filehead_config']['PLAN']['MIAVG']
@@ -614,7 +614,7 @@ def _reduce(data):
     theta = da['filehead_config']['BEAMCFGLIST'][0]
     if 'THETA=' in theta:
         da['beam_angle'] = int(theta[13:15])
-    
+
     tm = np.zeros((tmat['ROWS'], tmat['COLS']), dtype=np.float32)
     for irow in range(tmat['ROWS']):
         for icol in range(tmat['COLS']):

mhkit/dolfyn/io/nortek2_defs.py

@@ -300,7 +300,7 @@ def _calc_echo_struct(config, nc):
                                  'altraw', 'p_gd', 'std']]):
         raise Exception("Echosounder ping contains invalid data?")
     if flags['echo']:
-        dd += [('echo', 'H', [nc], _LinFunc(0.01, dtype=dt32), 'dB', 
+        dd += [('echo', 'H', [nc], _LinFunc(0.01, dtype=dt32), 'dB',
                 'Echo Sounder Acoustic Signal Backscatter', 'acoustic_target_strength_in_sea_water')]
     if flags['ahrs']:
         dd += _ahrs_def
@@ -323,14 +323,14 @@ def _calc_burst_struct(config, nb, nc):
     if flags['le']:
         # There may be a problem here with reading 32bit floats if
         # nb and nc are odd
-        dd += [('le_dist_alt', 'f', [], _LinFunc(dtype=dt32), 'm', 'Altimeter Range Leading Edge Algorithm', 
+        dd += [('le_dist_alt', 'f', [], _LinFunc(dtype=dt32), 'm', 'Altimeter Range Leading Edge Algorithm',
                 'altimeter_range'),
-               ('le_quality_alt', 'H', [], _LinFunc(0.01, dtype=dt32), 'dB', 
+               ('le_quality_alt', 'H', [], _LinFunc(0.01, dtype=dt32), 'dB',
                 'Altimeter Quality Indicator Leading Edge Algorithm'),
                ('status_alt', 'H', [], None, '1', 'Altimeter Status')]
     if flags['ast']:
         dd += [
-            ('ast_dist_alt', 'f', [], _LinFunc(dtype=dt32), 'm', 'Altimeter Range Acoustic Surface Tracking', 
+            ('ast_dist_alt', 'f', [], _LinFunc(dtype=dt32), 'm', 'Altimeter Range Acoustic Surface Tracking',
              'altimeter_range'),
             ('ast_quality_alt', 'H', [], _LinFunc(0.01, dtype=dt32), 'dB',
              'Altimeter Quality Indicator Acoustic Surface Tracking'),

mhkit/dolfyn/io/nortek2_lib.py

@@ -101,8 +101,8 @@ def _create_index(infile, outfile, N_ens, debug):
     fout = open(_abspath(outfile), 'wb')
     fout.write(b'Index Ver:')
     fout.write(struct.pack('<H', _index_version))
-    ids = [21, 22, 23, 24, 26, 28, 
-           27, 29, 30, 31, 35, 36]  
+    ids = [21, 22, 23, 24, 26, 28,
+           27, 29, 30, 31, 35, 36]
     # Saved: burst, avg, bt, vel_b5, alt_raw, echo
     # Not saved: bt record, DVL, alt record, avg alt_raw record, raw echo, raw echo transmit
     ens = dict.fromkeys(ids, 0)
@@ -117,7 +117,7 @@ def _create_index(infile, outfile, N_ens, debug):
             dat = _hdr.unpack(fin.read(_hdr.size))
         except:
             break
-        if dat[2] in ids:  
+        if dat[2] in ids:
             idk = dat[2]
             d_ver, d_off, config = struct.unpack('<BBH', fin.read(4))
             fin.seek(4, 1)
@@ -192,7 +192,7 @@ def _check_index(idx, infile, fix_hw_ens=False):
 
 def _boolarray_firstensemble_ping(index):
     """
-    Return a boolean of the index that indicates only the first ping in 
+    Return a boolean of the index that indicates only the first ping in
     each ensemble.
     """
     dens = np.ones(index['ens'].shape, dtype='bool')
@@ -456,7 +456,7 @@ def _collapse(vec, name=None, exclude=[]):
                           "Values found: {} (counts: {}).\n"
                           "Using the most common value: {}".format(
                               name, list(uniq), list(counts), val))
-    
+
         return val