Turbulence intensity and instrument noise (#293)

* New TI function, fix noise bug in psd

* Update notebook

* Fix test file

* Fix default

Author: jmcvey3

examples/adcp_example.ipynb

@@ -111,8 +111,10 @@ def __init__(
         n_fft_coh : int
           Number of data points to use for coherence and cross-spectra ffts
           Default: `n_fft_coh`=`n_fft`
-        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 `adp.turbulence.doppler_noise_level`.
+          Default: None.
         orientation : str, default='up'
           Instrument's orientation, either 'up' or 'down'
         diff_style : str, default='centered_extended'
@@ -348,9 +350,10 @@ def doppler_noise_level(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",
@@ -869,7 +872,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[0.2, 0.4]):
 
         return m, b
 
-    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
+    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4], noise=None):
         """
         Calculate the TKE dissipation rate from the velocity spectra.
 
@@ -882,6 +885,10 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
         f_range : iterable(2)
           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.doppler_noise_level`.
+          Default: None.
 
         Returns
         -------
@@ -918,6 +925,17 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[0.2, 0.4]):
             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._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]))

examples/data/dolfyn/test_data/Sig1000_IMU_bin.nc

@@ -24,8 +24,10 @@ class ADVBinner(VelBinner):
     n_fft_coh : int
       Number of data points to use for coherence and cross-spectra fft's.
       Optional, default `n_fft_coh` = `n_fft`
-    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.doppler_noise_level`.
+          Default: None.
     """
 
     def __call__(self, ds, freq_units="rad/s", window="hann"):
@@ -263,7 +265,7 @@ def doppler_noise_level(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",
@@ -337,7 +339,7 @@ def check_turbulence_cascade_slope(self, psd, freq_range=[6.28, 12.57]):
 
         return m, b
 
-    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
+    def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57], noise=None):
         """
         Calculate the dissipation rate from the PSD
 
@@ -351,6 +353,10 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
           The range over which to integrate/average the spectrum, in units
           of the psd frequency vector (Hz or rad/s).
           Default = [6.28, 12.57] 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
         -------
@@ -390,6 +396,18 @@ def dissipation_rate_LT83(self, psd, U_mag, freq_range=[6.28, 12.57]):
         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]))
         idx = idx[0]

examples/data/dolfyn/test_data/Sig1000_tidal_bin.nc

@@ -421,7 +421,7 @@ def _psd_base(
 
         for slc in slice1d_along_axis(dat.shape, -1):
             out[slc] = psd_1D(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

examples/data/dolfyn/test_data/vector_data01_bin.nc

@@ -803,6 +803,60 @@ def autocovariance(self, veldat, n_bin=None):
 
         return da
 
+    def turbulence_intensity(self, U_mag, noise=0, thresh=0, detrend=False):
+        """
+        Calculate noise-corrected turbulence intensity.
+
+        Parameters
+        ----------
+        U_mag : xarray.DataArray
+          Raw horizontal velocity magnitude
+        noise : numeric
+          Instrument noise level in same units as velocity. Typically
+          found from `<adv or adp>.turbulence.doppler_noise_level`.
+          Default: None.
+        thresh : numeric
+          Theshold below which TI will not be calculated
+        detrend : bool (default: False)
+          Detrend the velocity data (True), or simply de-mean it
+          (False), prior to computing TI.
+        """
+
+        if "xarray" in type(U_mag).__module__:
+            U = U_mag.values
+        if "xarray" in type(noise).__module__:
+            noise = noise.values
+
+        if detrend:
+            up = self.detrend(U)
+        else:
+            up = self.demean(U)
+
+        # Take RMS and subtract noise
+        u_rms = np.sqrt(np.nanmean(up**2, axis=-1) - noise**2)
+        u_mag = self.mean(U)
+
+        ti = np.ma.masked_where(u_mag < thresh, u_rms / u_mag)
+
+        dims = U_mag.dims
+        coords = {}
+        for nm in U_mag.dims:
+            if "time" in nm:
+                coords[nm] = self.mean(U_mag[nm].values)
+            else:
+                coords[nm] = U_mag[nm].values
+
+        return xr.DataArray(
+            ti.data.astype("float32"),
+            coords=coords,
+            dims=dims,
+            attrs={
+                "units": "% [0,1]",
+                "long_name": "Turbulence Intensity",
+                "comment": f"TI was corrected from a noise level of {noise} m/s",
+            },
+        )
+
     def turbulent_kinetic_energy(self, veldat, noise=None, detrend=True):
         """
         Calculate the turbulent kinetic energy (TKE) (variances
@@ -814,11 +868,12 @@ def turbulent_kinetic_energy(self, veldat, noise=None, detrend=True):
           Velocity data array from ADV or single beam from ADCP.
           The last dimension is assumed to be time.
         noise : float or array-like
-          A vector of the noise levels of the velocity data with
-          the same first dimension as the velocity vector.
+          Instrument noise level in same units as velocity. Typically
+          found from `<adv or adp>.turbulence.doppler_noise_level`.
+          Default: None.
         detrend : bool (default: False)
           Detrend the velocity data (True), or simply de-mean it
-          (False), prior to computing tke. Note: the psd routines
+          (False), prior to computing TKE. Note: the PSD routines
           use detrend, so if you want to have the same amount of
           variance here as there use ``detrend=True``.
 
@@ -892,7 +947,7 @@ def power_spectral_density(
         freq_units="rad/s",
         fs=None,
         window="hann",
-        noise=None,
+        noise=0,
         n_bin=None,
         n_fft=None,
         n_pad=None,
@@ -913,10 +968,9 @@ def power_spectral_density(
         window : string or array
           Specify the window function.
           Options: 1, None, 'hann', 'hamm'
-        noise : float or array-like
-          A vector of the noise levels of the velocity data with
-          the same first dimension as the velocity vector.
-          Default = 0.
+        noise : numeric or array
+          Instrument noise level in same units as velocity.
+          Default: 0 (ADCP) or [0, 0, 0] (ADV).
         n_bin : int (optional)
           The bin-size. Default: from the binner.
         n_fft : int (optional)
@@ -938,8 +992,6 @@ def power_spectral_density(
         n_fft = self._parse_nfft(n_fft)
         if "xarray" in type(veldat).__module__:
             vel = veldat.values
-        if "xarray" in type(noise).__module__:
-            noise = noise.values
         if ("rad" not in freq_units) and ("Hz" not in freq_units):
             raise ValueError("`freq_units` should be one of 'Hz' or 'rad/s'")
 
@@ -964,16 +1016,19 @@ def power_spectral_density(
         ).astype("float32")
 
         # Spectra, if input is full velocity or a single array
-        if len(vel.shape) == 2:
-            if not vel.shape[0] == 3:
-                raise Exception(
+        if len(vel.shape) >= 2:
+            if vel.shape[0] != 3:
+                raise ValueError(
                     "Function can only handle 1D or 3D arrays."
                     " If ADCP data, please select a specific depth bin."
                 )
-            if (noise is not None) and (np.shape(noise)[0] != 3):
-                raise Exception("Noise should have same first dimension as velocity")
+            if np.array(noise).any():
+                if np.size(noise) != 3:
+                    raise ValueError("Noise is expected to be an array of 3 scalars")
             else:
+                # Reset default to list of 3 zeros
                 noise = np.array([0, 0, 0])
+
             out = np.empty(
                 self._outshape_fft(vel[:3].shape, n_fft=n_fft, n_bin=n_bin),
                 dtype=np.float32,
@@ -996,10 +1051,9 @@ def power_spectral_density(
             }
             dims = ["S", "time", "freq"]
         else:
-            if (noise is not None) and (len(np.shape(noise)) > 1):
-                raise Exception("Noise should have same first dimension as velocity")
-            else:
-                noise = np.array(0)
+            if np.array(noise).any() and np.size(noise) > 1:
+                raise ValueError("Noise is expected to be a scalar")
+
             out = self._psd_base(
                 vel,
                 fs=fs,

mhkit/dolfyn/adp/turbulence.py

@@ -92,22 +92,29 @@ def test_adv_turbulence(self):
         dat = tv.dat.copy(deep=True)
         bnr = avm.ADVBinner(n_bin=20.0, fs=dat.fs)
         tdat = bnr(dat)
-        acov = bnr.autocovariance(dat.vel)
+        acov = bnr.autocovariance(dat["vel"])
 
         assert_identical(tdat, avm.turbulence_statistics(dat, n_bin=20.0, fs=dat.fs))
 
-        tdat["stress_detrend"] = bnr.reynolds_stress(dat.vel)
-        tdat["stress_demean"] = bnr.reynolds_stress(dat.vel, detrend=False)
+        tdat["stress_detrend"] = bnr.reynolds_stress(dat["vel"])
+        tdat["stress_demean"] = bnr.reynolds_stress(dat["vel"], detrend=False)
         tdat["csd"] = bnr.cross_spectral_density(
-            dat.vel, freq_units="rad", window="hamm", n_fft_coh=10
+            dat["vel"], freq_units="rad", window="hamm", n_fft_coh=10
         )
-        tdat["LT83"] = bnr.dissipation_rate_LT83(tdat.psd, tdat.velds.U_mag)
-        tdat["SF"] = bnr.dissipation_rate_SF(dat.vel[0], tdat.velds.U_mag)
+        tdat["LT83"] = bnr.dissipation_rate_LT83(tdat["psd"], tdat.velds.U_mag)
+        tdat["noise"] = bnr.doppler_noise_level(tdat["psd"], pct_fN=0.8)
+        tdat["LT83_noise"] = bnr.dissipation_rate_LT83(
+            tdat["psd"], tdat.velds.U_mag, noise=tdat["noise"]
+        )
+        tdat["SF"] = bnr.dissipation_rate_SF(dat["vel"][0], tdat.velds.U_mag)
         tdat["TE01"] = bnr.dissipation_rate_TE01(dat, tdat)
         tdat["L"] = bnr.integral_length_scales(acov, tdat.velds.U_mag)
         slope_check = bnr.check_turbulence_cascade_slope(
             tdat["psd"][-1].mean("time"), freq_range=[10, 100]
         )
+        tdat["psd_noise"] = bnr.power_spectral_density(
+            dat["vel"], freq_units="rad", noise=[0.06, 0.04, 0.01]
+        )
 
         if make_data:
             save(tdat, "vector_data01_bin.nc")
@@ -117,13 +124,22 @@ def test_adv_turbulence(self):
         assert_allclose(tdat, load("vector_data01_bin.nc"), atol=1e-6)
 
     def test_adcp_turbulence(self):
-        dat = tr.dat_sig_i.copy(deep=True)
+        dat = tr.dat_sig_tide.copy(deep=True)
+        dat.velds.rotate2("earth")
+        dat.attrs["principal_heading"] = apm.calc_principal_heading(
+            dat.vel.mean("range")
+        )
         bnr = apm.ADPBinner(n_bin=20.0, fs=dat.fs, diff_style="centered")
         tdat = bnr.bin_average(dat)
-        tdat["dudz"] = bnr.dudz(tdat.vel)
-        tdat["dvdz"] = bnr.dvdz(tdat.vel)
-        tdat["dwdz"] = bnr.dwdz(tdat.vel)
-        tdat["tau2"] = bnr.shear_squared(tdat.vel)
+
+        tdat["dudz"] = bnr.dudz(tdat["vel"])
+        tdat["dvdz"] = bnr.dvdz(tdat["vel"])
+        tdat["dwdz"] = bnr.dwdz(tdat["vel"])
+        tdat["tau2"] = bnr.shear_squared(tdat["vel"])
+        tdat["I"] = tdat.velds.I
+        tdat["ti"] = bnr.turbulence_intensity(dat.velds.U_mag, detrend=False)
+        dat.velds.rotate2("beam")
+
         tdat["psd"] = bnr.power_spectral_density(
             dat["vel"].isel(dir=2, range=len(dat.range) // 2), freq_units="Hz"
         )
@@ -137,13 +153,22 @@ def test_adcp_turbulence(self):
         tdat["tke"] = bnr.total_turbulent_kinetic_energy(
             dat, noise=tdat["noise"], orientation="up", beam_angle=25
         )
+        tdat["ti_noise"] = bnr.turbulence_intensity(
+            dat.velds.U_mag, detrend=False, noise=tdat["noise"]
+        )
         # This is "negative" for this code check
         tdat["wpwp"] = bnr.turbulent_kinetic_energy(dat["vel_b5"], noise=tdat["noise"])
         tdat["dissipation_rate_LT83"] = bnr.dissipation_rate_LT83(
             tdat["psd"],
             tdat.velds.U_mag.isel(range=len(dat.range) // 2),
             freq_range=[0.2, 0.4],
         )
+        tdat["dissipation_rate_LT83_noise"] = bnr.dissipation_rate_LT83(
+            tdat["psd"],
+            tdat.velds.U_mag.isel(range=len(dat.range) // 2),
+            freq_range=[0.2, 0.4],
+            noise=tdat["noise"],
+        )
         (
             tdat["dissipation_rate_SF"],
             tdat["noise_SF"],
@@ -155,10 +180,22 @@ def test_adcp_turbulence(self):
         slope_check = bnr.check_turbulence_cascade_slope(
             tdat["psd"].mean("time"), freq_range=[0.4, 4]
         )
+        tdat["psd_noise"] = bnr.power_spectral_density(
+            dat["vel"].isel(dir=2, range=len(dat.range) // 2),
+            freq_units="Hz",
+            noise=0.01,
+        )
 
         if make_data:
-            save(tdat, "Sig1000_IMU_bin.nc")
+            save(tdat, "Sig1000_tidal_bin.nc")
             return
 
+        with pytest.raises(Exception):
+            bnr.calc_psd(dat["vel"], freq_units="Hz", noise=0.01)
+
+        with pytest.raises(Exception):
+            bnr.calc_psd(dat["vel"][0], freq_units="Hz", noise=0.01)
+
         assert np.round(slope_check[0].values, 4), -1.0682
-        assert_allclose(tdat, load("Sig1000_IMU_bin.nc"), atol=1e-6)
+
+        assert_allclose(tdat, load("Sig1000_tidal_bin.nc"), atol=1e-6)