Fixes for Nortek Altimeter Measurements (lkilcher#125)

* Updates

* Nortek AWAC AST support (#14)

* Update awac ast functions

* sample conversion

* Save bug

* Actually fix save bug

* Fix wave variables and update burst_config

* Fix wave variables and update burst_config part 2

* Except raw altimeter data from averaging

* Reset test

* Update changelog

changelog.md

@@ -12,8 +12,11 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 		- Fix for instruments that did not record water velocity
 		- Fix netCDF4 compression encoding
 		- Retain prior netCDF4 variable encoding
+		- Fix bug in reading raw Nortek Signature altimeter data
+
 	- API/Useability
 	    - Updates to support python 3.10 and 3.11
+		- Added ability to read Nortek AWAC waves data
 
 ## Version 1.3.0
     - Bugfixes

dolfyn/io/api.py

@@ -229,13 +229,15 @@ def load(filename):
             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
 

dolfyn/io/base.py

@@ -116,7 +116,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:
@@ -203,6 +203,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']:

dolfyn/io/nortek.py

@@ -163,8 +163,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,
@@ -1021,7 +1024,78 @@ def sci_awac_profile(self,):
         self.data['coords']['range'] = r
         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)()

dolfyn/io/nortek2.py

@@ -216,9 +216,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
@@ -232,6 +237,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):
@@ -271,13 +277,16 @@ def readfile(self, ens_start=0, ens_stop=None):
             except IOError:
                 return outdat
             id = hdr['id']
-            if id in [21, 22, 23, 24, 28]:  # vel, bt, vel_b5, echo
+            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]:  # alt_raw (altimeter burst)
+            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'):
                     first_pass = True
-                    tmp_idx = rdr._nsamp_index = rdr._names.index('altraw_nsamp')  # noqa
+                    tmp_idx = rdr._nsamp_index = rdr._names.index('nsamp_alt')
                     shift = rdr._nsamp_shift = calcsize(
                         defs._format(rdr._format[:tmp_idx],
                                      rdr._N[:tmp_idx]))
@@ -299,9 +308,9 @@ def readfile(self, ens_start=0, ens_stop=None):
                     rdr._cs_struct = defs.Struct(
                         '<' + '{}H'.format(int(rdr.nbyte // 2)))
                     # Initialize the array
-                    outdat[26]['altraw_samp'] = defs._nans(
+                    outdat[26]['samp_alt'] = defs._nans(
                         [rdr._N[tmp_idx],
-                         len(outdat[26]['altraw_samp'])],
+                         len(outdat[26]['samp_alt'])],
                         dtype=np.uint16)
                 else:
                     if sz != rdr._N[tmp_idx]:
@@ -312,8 +321,9 @@ 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]:  # bt record, DVL,
-                # alt record, avg alt_raw record, raw echo, raw echo transmit
+            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:
                     logging.debug(
                         "Skipped ID: 0x{:02X} ({:02d})\n".format(id, id))
@@ -399,7 +409,7 @@ def _reorg(dat):
     cfg['inst_type'] = 'ADCP'
 
     for id, tag in [(21, ''), (22, '_avg'), (23, '_bt'), 
-                    (24, '_b5'), (26, '_ast'), (28, '_echo')]:
+                    (24, '_b5'), (26, 'raw'), (28, '_echo')]:
         if id in [24, 26]:
             collapse_exclude = [0]
         else:
@@ -451,10 +461,10 @@ def _reorg(dat):
                 outdat['standard_name'][ky + tag] = 'number_of_observations'
 
         for ky in ['vel', 'amp', 'corr', 'prcnt_gd', 'echo', 'dist',
-                   'orientmat', 'angrt', 'quaternions', 'ast_pressure',
-                   'alt_dist', 'alt_quality', 'alt_status',
-                   'ast_dist', 'ast_quality', 'ast_offset_time',
-                   'altraw_nsamp', 'altraw_dsamp', 'altraw_samp',
+                   'orientmat', 'angrt', 'quaternions', 'pressure_alt',
+                   'le_dist_alt', 'le_quality_alt', 'status_alt',
+                   'ast_dist_alt', 'ast_quality_alt', 'ast_offset_time_alt',
+                   'nsamp_alt', 'dsamp_alt', 'samp_alt',
                    'status0', 'fom', 'temp_press', 'press_std',
                    'pitch_std', 'roll_std', 'heading_std', 'xmit_energy',
                    ]:
@@ -463,20 +473,19 @@ def _reorg(dat):
 
     # Move 'altimeter raw' data to its own down-sampled structure
     if 26 in dat:
-        ard = outdat['altraw']
         for ky in list(outdat['data_vars']):
-            if ky.endswith('_ast'):
-                grp = ky.split('.')[0]
-                if '.' in ky and grp not in ard:
-                    ard[grp] = {}
-                ard[ky.rstrip('_ast')] = outdat['data_vars'].pop(ky)
+            if ky.endswith('raw') and not ky.endswith('_altraw'):
+                 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
 
         # Read altimeter status
-        alt_status = lib._alt_status2data(outdat['data_vars']['alt_status'])
-        for ky in alt_status:
+        outdat['data_vars'].pop('status_altraw')
+        status_alt = lib._alt_status2data(outdat['data_vars']['status_alt'])
+        for ky in status_alt:
             outdat['attrs'][ky] = lib._collapse(
-                alt_status[ky].astype('uint8'), name=ky)
-        outdat['data_vars'].pop('alt_status')
+                status_alt[ky].astype('uint8'), name=ky)
+        outdat['data_vars'].pop('status_alt')
 
         # Power level index
         power = {0: 'high', 1: 'med-high', 2: 'med-low', 3: 'low'}

dolfyn/io/nortek2_defs.py

@@ -296,8 +296,8 @@ def _calc_echo_struct(config, nc):
     flags = lib._headconfig_int2dict(config)
     dd = copy(_burst_hdr)
     dd[19] = ('blank_dist', 'H', [], _LinFunc(0.001))  # m
-    if any([flags[nm] for nm in ['vel', 'amp', 'corr', 'alt', 'ast',
-                                 'alt_raw', 'p_gd', 'std']]):
+    if any([flags[nm] for nm in ['vel', 'amp', 'corr', 'le', 'ast',
+                                 'altraw', 'p_gd', 'std']]):
         raise Exception("Echosounder ping contains invalid data?")
     if flags['echo']:
         dd += [('echo', 'H', [nc], _LinFunc(0.01, dtype=dt32), 'dB', 
@@ -320,30 +320,33 @@ def _calc_burst_struct(config, nb, nc):
     if flags['corr']:
         dd.append(('corr', 'B', [nb, nc], None, '%', 'Acoustic Signal Correlation',
                   'beam_consistency_indicator_from_multibeam_acoustic_doppler_velocity_profiler_in_sea_water'))
-    if flags['alt']:
+    if flags['le']:
         # There may be a problem here with reading 32bit floats if
         # nb and nc are odd?
-        dd += [('alt_dist', 'f', [], _LinFunc(dtype=dt32), 'm', 'Altimeter Range', 'altimeter_range'),
-               ('alt_quality', 'H', [], _LinFunc(0.01, dtype=dt32), '1', 'Altimeter Quality Indicator'),
-               ('alt_status', 'H', [], None, '1', 'Altimeter Status')]
+        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', 
+                'Altimeter Quality Indicator Leading Edge Algorithm'),
+               ('status_alt', 'H', [], None, '1', 'Altimeter Status')]
     if flags['ast']:
         dd += [
-            ('ast_dist', 'f', [], _LinFunc(dtype=dt32), 'm', 'Acoustic Surface Tracking Range'),
-            ('ast_quality', 'H', [], _LinFunc(0.01, dtype=dt32), '1',
-             'Acoustic Surface Tracking Quality Indicator'),
-            ('ast_offset_time', 'h', [], _LinFunc(0.0001, dtype=dt32),
+            ('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'),
+            ('ast_offset_time_alt', 'h', [], _LinFunc(0.0001, dtype=dt32),
              's', 'Acoustic Surface Tracking Time Offset to Velocity Ping'),
-            ('ast_pressure', 'f', [], None, 'dbar', 'Pressure measured during AST ping',
+            ('pressure_alt', 'f', [], None, 'dbar', 'Pressure measured during AST ping',
              'sea_water_pressure'),
             # This use of 'x' here is a hack
-            ('ast_spare', 'B7x', [], None),
+            ('spare', 'B7x', [], None),
         ]
-    if flags['alt_raw']:
+    if flags['altraw']:
         dd += [
-            ('altraw_nsamp', 'I', [], None, '1', 'Number of Altimeter Samples'),
-            ('altraw_dsamp', 'H', [], _LinFunc(0.0001, dtype=dt32), 'm',
+            ('nsamp_alt', 'I', [], None, '1', 'Number of Altimeter Samples'),
+            ('dsamp_alt', 'H', [], _LinFunc(0.0001, dtype=dt32), 'm',
              'Altimeter Distance between Samples'),
-            ('altraw_samp', 'h', [], None),
+            ('samp_alt', 'h', [], None, '1', 'Altimeter Samples'),
         ]
     if flags['ahrs']:
         dd += _ahrs_def

dolfyn/io/nortek2_lib.py

@@ -336,8 +336,8 @@ def _headconfig_int2dict(val, mode='burst'):
             vel=_getbit(val, 5),
             amp=_getbit(val, 6),
             corr=_getbit(val, 7),
-            alt=_getbit(val, 8),
-            alt_raw=_getbit(val, 9),
+            le=_getbit(val, 8),
+            altraw=_getbit(val, 9),
             ast=_getbit(val, 10),
             echo=_getbit(val, 11),
             ahrs=_getbit(val, 12),