Merge 425806a into 146049a

changelog.md

@@ -13,10 +13,16 @@ and this project adheres to [Semantic Versioning](http://semver.org/spec/v2.0.0.
 		- Fix netCDF4 compression encoding
 		- Retain prior netCDF4 variable encoding
 		- Fix bug in reading raw Nortek Signature altimeter data
+		- Fix bug where noise input wasn't being subtracted from auto-spectra
+		- Fix bug that would error out when entering custom FFT window
 
 	- API/Useability
 	    - Updates to support python 3.10 and 3.11
 		- Added ability to read Nortek AWAC waves data
+		- Added ability to subtract Doppler noise in TKE dissipation rate functions
+		- Added function to calculate turbulence intensity and remove noise
+		- Add ability to read Nortek dual profiling instruments
+		- Add ability to read ID 31 (initial altimeter scan for averaged altimeter measurements)
 
 	- Nortek Vectrino (.vno)
 		- Add support for Nortek Vectrino (.vno) files.

docs/apidoc/dolfyn.binners.rst

@@ -21,12 +21,12 @@ Below is a list of functions that can be called from `VelBinner`.
 	~dolfyn.binned.TimeBinner.std
 	~dolfyn.velocity.VelBinner.do_avg
 	~dolfyn.velocity.VelBinner.do_var
+	~dolfyn.velocity.VelBinner.calc_ti
+	~dolfyn.velocity.VelBinner.calc_psd
 	~dolfyn.velocity.VelBinner.calc_coh
 	~dolfyn.velocity.VelBinner.calc_phase_angle
 	~dolfyn.velocity.VelBinner.calc_acov
 	~dolfyn.velocity.VelBinner.calc_xcov
-	~dolfyn.velocity.VelBinner.calc_tke
-	~dolfyn.velocity.VelBinner.calc_psd
 	~dolfyn.binned.TimeBinner.calc_freq
 	~dolfyn.binned.TimeBinner.calc_psd_base
 	~dolfyn.binned.TimeBinner.calc_csd_base
@@ -43,6 +43,7 @@ Functions for analyzing ADV data via the `ADVBinner` class, beyond those describ
 	~dolfyn.adv.turbulence.ADVBinner
 	~dolfyn.adv.turbulence.calc_turbulence
 	~dolfyn.adv.turbulence.ADVBinner.calc_csd
+	~dolfyn.velocity.VelBinner.calc_tke
 	~dolfyn.adv.turbulence.ADVBinner.calc_stress
 	~dolfyn.adv.turbulence.ADVBinner.calc_doppler_noise
 	~dolfyn.adv.turbulence.ADVBinner.check_turbulence_cascade_slope
@@ -64,7 +65,6 @@ Functions for analyzing ADCP data via the `ADPBinner` class, beyond those descri
 	~dolfyn.adp.turbulence.ADPBinner.calc_doppler_noise
 	~dolfyn.adp.turbulence.ADPBinner.calc_stress_4beam
 	~dolfyn.adp.turbulence.ADPBinner.calc_stress_5beam
-	~dolfyn.adp.turbulence.ADPBinner.calc_total_tke
 	~dolfyn.adp.turbulence.ADPBinner.check_turbulence_cascade_slope
 	~dolfyn.adp.turbulence.ADPBinner.calc_dissipation_LT83
 	~dolfyn.adp.turbulence.ADPBinner.calc_dissipation_SF

dolfyn/adp/turbulence.py

@@ -238,9 +238,10 @@ def calc_doppler_noise(self, psd, pct_fN=0.8):
         N2 = psd.sel(freq=f_range) * psd.freq.sel(freq=f_range)
         noise_level = np.sqrt(N2.mean(dim='freq'))
 
+        time_coord = psd.dims[0]  # no reason this shouldn't be time or time_b5
         return xr.DataArray(
             noise_level.values.astype('float32'),
-            dims=['time'],
+            coords={time_coord: psd.coords[time_coord]},
             attrs={'units': 'm s-1',
                    'long_name': 'Doppler Noise Level',
                    'description': 'Doppler noise level calculated '
@@ -439,7 +440,7 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t
         in pitch and roll. u'v'_ cannot be directly calculated by a 5-beam ADCP,
         so it is approximated by the covariance of `u` and `v`. The uncertainty
         introduced by using this approximation is small if deviations from pitch
-        and roll are small (< 10 degrees).
+        and roll are small (<= 5 degrees).
 
         Dewey, R., and S. Stringer. "Reynolds stresses and turbulent kinetic
         energy estimates from various ADCP beam configurations: Theory." J. of
@@ -455,7 +456,7 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t
 
         # Run through warnings
         b_angle, noise = self._stress_func_warnings(
-            ds, beam_angle, noise, tilt_thresh=10)
+            ds, beam_angle, noise, tilt_thresh=5)
 
         # Fetch beam order
         beam_order, phi2, phi3 = self._check_orientation(
@@ -522,47 +523,6 @@ def calc_stress_5beam(self, ds, noise=None, orientation=None, beam_angle=None, t
 
             return tke_vec, stress_vec
 
-    def calc_total_tke(self, ds, noise=None, orientation=None, beam_angle=None):
-        """Calculate magnitude of turbulent kinetic energy from 5-beam ADCP. 
-
-        Parameters
-        ----------
-        ds : xarray.Dataset
-          Raw dataset in beam coordinates
-        ds_avg : xarray.Dataset
-          Binned dataset in final coordinate reference frame
-        noise : int or xarray.DataArray, default=0 (time)
-          Doppler noise level in units of m/s
-        orientation : str, default=ds.attrs['orientation']
-          Direction ADCP is facing ('up' or 'down')
-        beam_angle : int, default=ds.attrs['beam_angle']
-          ADCP beam angle in units of degrees
-
-        Returns
-        -------
-        tke : xarray.DataArray
-          Turbulent kinetic energy magnitude
-
-        Notes
-        -----
-        This function is a wrapper around 'calc_stress_5beam' that then
-        combines the TKE components.
-
-        Warning: the integral length scale of turbulence captured by the
-        ADCP measurements (i.e. the size of turbulent structures) increases 
-        with increasing range from the instrument.
-        """
-
-        tke_vec = self.calc_stress_5beam(
-            ds, noise, orientation, beam_angle, tke_only=True)
-
-        tke = tke_vec.sum('tke') / 2
-        tke.attrs['units'] = 'm2 s-2'
-        tke.attrs['long_name'] = 'TKE Magnitude',
-        tke.attrs['standard_name'] = 'specific_turbulent_kinetic_energy_of_sea_water'
-
-        return tke.astype('float32')
-
     def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):
         """This function calculates the slope of the PSD, the power spectra 
         of velocity, within the given frequency range. The purpose of this
@@ -623,7 +583,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):
 
         return m, b
 
-    def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
+    def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
         """Calculate the TKE dissipation rate from the velocity spectra.
 
         Parameters
@@ -633,8 +593,12 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
         U_mag : xarray.DataArray (time)
           The bin-averaged horizontal velocity (a.k.a. speed) from a single depth bin (range)
         f_range : iterable(2)
-          The range over which to integrate/average the spectrum, in units 
+          The range over which to integrate/average the spectrum, in units
           of the psd frequency vector (Hz or rad/s)
+        noise : float or array-like
+          Instrument noise level in same units as velocity. Typically
+          found from `adp.turbulence.calc_doppler_noise`. 
+          Default: None.
 
         Returns
         -------
@@ -669,6 +633,20 @@ def calc_dissipation_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
             raise Exception('PSD should be 2-dimensional (time, frequency)')
         if len(U_mag.shape) != 1:
             raise Exception('U_mag should be 1-dimensional (time)')
+        if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
+            raise ValueError("`freq_range` must be an iterable of length 2.")
+        if noise is not None:
+            if np.shape(noise)[0] != np.shape(psd)[0]:
+                raise Exception(
+                    'Noise should have same first dimension as PSD')
+        else:
+            noise = np.array(0)
+
+        # Noise subtraction from binner.TimeBinner.calc_psd_base
+        psd = psd.copy()
+        if noise is not None:
+            psd -= noise**2 / (self.fs / 2)
+            psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)
 
         freq = psd.freq
         idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))
@@ -853,7 +831,7 @@ def calc_ustar_fit(self, ds_avg, upwp_, z_inds=slice(1, 5), H=None):
         """
 
         if not H:
-            H = ds_avg.depth.values
+            H = ds_avg["depth"].values
         z = ds_avg['range'].values
         upwp_ = upwp_.values
 
@@ -863,6 +841,6 @@ def calc_ustar_fit(self, ds_avg, upwp_, z_inds=slice(1, 5), H=None):
 
         return xr.DataArray(
             u_star.astype('float32'),
-            coords={'time': ds_avg.time},
+            coords={'time': ds_avg["time"]},
             attrs={'units': 'm s-1',
                    'long_name': 'Friction Velocity'})

dolfyn/adv/turbulence.py

@@ -21,8 +21,10 @@ class ADVBinner(VelBinner):
       The length of the FFT for computing spectra (must be <= n_bin)
     n_fft_coh : int (optional, default: `n_fft_coh`=`n_fft`)
       Number of data points to use for coherence and cross-spectra ffts
-    noise : float, list or numpy.ndarray
-      Instrument's doppler noise in same units as velocity
+    noise : float or array-like
+      Instrument noise level in same units as velocity. Typically
+      found from `adv.turbulence.calc_doppler_noise`. 
+      Default: None.
     """
 
     def __call__(self, ds, freq_units='rad/s', window='hann'):
@@ -49,7 +51,7 @@ def __call__(self, ds, freq_units='rad/s', window='hann'):
 
     def calc_stress(self, veldat, detrend=True):
         """
-        Calculate the stresses (covariances of u,v,w)
+        Calculate the stresses (covariances of u,v,w in m^2/s^2)
 
         Parameters
         ----------
@@ -165,7 +167,7 @@ def calc_csd(self, veldat,
                                    'time': time,
                                    'coh_freq': coh_freq},
                            dims=['C', 'time', 'coh_freq'],
-                           attrs={'units': units, 
+                           attrs={'units': units,
                                   'n_fft_coh': n_fft,
                                   'long_name': 'Cross Spectral Density'})
         csd['coh_freq'].attrs['units'] = freq_units
@@ -230,7 +232,7 @@ def calc_doppler_noise(self, psd, pct_fN=0.8):
 
         return xr.DataArray(
             noise_level.values.astype('float32'),
-            dims=['dir', 'time'],
+            coords={'S': psd['S'], 'time': psd['time']},
             attrs={'units': 'm/s',
                    'long_name': 'Doppler Noise Level',
                    'description': 'Doppler noise level calculated '
@@ -296,7 +298,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[6.28, 12.57]):
 
         return m, b
 
-    def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
+    def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None):
         """Calculate the dissipation rate from the PSD
 
         Parameters
@@ -308,6 +310,10 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
         freq_range : iterable(2) (default: [6.28, 12.57])
           The range over which to integrate/average the spectrum, in units 
           of the psd frequency vector (Hz or rad/s)
+        noise : float or array-like
+          Instrument noise level in same units as velocity. Typically
+          found from `adv.turbulence.calc_doppler_noise`. 
+          Default: None.
 
         Returns
         -------
@@ -339,8 +345,23 @@ def calc_epsilon_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
         """
 
         # Ensure time has been averaged
-        if len(psd.time)!=len(U_mag.time):
+        if len(psd.time) != len(U_mag.time):
             raise Exception("`U_mag` should be from ensembled-averaged dataset")
+        if not hasattr(freq_range, "__iter__") or len(freq_range) != 2:
+            raise ValueError("`freq_range` must be an iterable of length 2.")
+
+        if noise is not None:
+            if np.shape(noise)[0] != 3:
+                raise Exception(
+                    'Noise should have same first dimension as velocity')
+        else:
+            noise = np.array([0, 0, 0])[:, None, None]
+
+        # Noise subtraction from binner.TimeBinner.calc_psd_base
+        psd = psd.copy()
+        if noise is not None:
+            psd -= (noise**2 / (self.fs / 2))
+            psd = psd.where(psd > 0, np.min(np.abs(psd)) / 100)
 
         freq = psd.freq
         idx = np.where((freq_range[0] < freq) & (freq < freq_range[1]))

dolfyn/binned.py

@@ -401,7 +401,7 @@ def calc_psd_base(self, dat, fs=None, window='hann', noise=0,
         for slc in slice1d_along_axis(dat.shape, -1):
             out[slc] = psd(dat[slc], n_fft, fs,
                            window=window, step=step)
-        if noise != 0:
+        if np.any(noise):
             out -= noise**2 / (fs/2)
             # Make sure all values of the PSD are >0 (but still small):
             out[out < 0] = np.min(np.abs(out)) / 100

dolfyn/io/api.py

@@ -192,7 +192,7 @@ def save(ds, filename,
 
         if compression:
             # New netcdf4-c cannot compress variable length strings
-            if isinstance(ds[ky].data[0], str):
+            if ds[ky].size <= 1 or isinstance(ds[ky].data[0], str):
                 continue
             enc[ky].update(dict(zlib=True, complevel=1))
 
@@ -219,6 +219,12 @@ def load(filename):
     """
 
     filename = _check_file_ext(filename, 'nc')
+
+    file_type = _get_filetype(filename)
+    if file_type == '<GIT-LFS pointer>':
+        raise IOError("File '{}' looks like a git-lfs pointer. You may need to "
+                      "install and initialize git-lfs. See https://git-lfs.github.com"
+                      " for details.".format(filename))
 
     ds = xr.load_dataset(filename, engine='netcdf4')