Merge 164e254 into 184a5ab

mhkit/dolfyn/adp/turbulence.py

@@ -646,7 +646,7 @@ def stress_tensor_5beam(
         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
@@ -663,7 +663,7 @@ def stress_tensor_5beam(
 
         # 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
@@ -763,49 +763,6 @@ def stress_tensor_5beam(
 
             return tke_vec, stress_vec
 
-    def total_turbulent_kinetic_energy(
-        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
-        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.stress_tensor_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
@@ -1139,7 +1096,7 @@ def friction_velocity(self, ds_avg, upwp_, z_inds=slice(1, 5), H=None):
             raise TypeError("`z_inds` must be an instance of `slice(int,int)`.")
 
         if not H:
-            H = ds_avg.depth.values
+            H = ds_avg["depth"].values
         z = ds_avg["range"].values
         upwp_ = upwp_.values
 
@@ -1151,6 +1108,6 @@ def friction_velocity(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"},
         )

mhkit/dolfyn/adv/turbulence.py

@@ -54,7 +54,7 @@ def __call__(self, ds, freq_units="rad/s", window="hann"):
     def reynolds_stress(self, veldat, detrend=True):
         """
         Calculate the specific Reynolds stresses
-        (cross-covariances of u,v,w in m^2/s^2)
+        (covariances of u,v,w in m^2/s^2)
 
         Parameters
         ----------

mhkit/tests/dolfyn/test_analysis.py

@@ -130,14 +130,16 @@ def test_adcp_turbulence(self):
             dat.vel.mean("range")
         )
         bnr = apm.ADPBinner(n_bin=20.0, fs=dat.fs, diff_style="centered")
+        U_mag = dat.velds.U_mag
+        dat["U_mag"] = U_mag
         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["I"] = tdat.velds.I
-        tdat["ti"] = bnr.turbulence_intensity(dat.velds.U_mag, detrend=False)
+        tdat["ti"] = bnr.turbulence_intensity(U_mag, detrend=False)
         dat.velds.rotate2("beam")
 
         tdat["psd"] = bnr.power_spectral_density(
@@ -150,41 +152,44 @@ def test_adcp_turbulence(self):
         tdat["tke_vec5"], tdat["stress_vec5"] = bnr.stress_tensor_5beam(
             dat, noise=tdat["noise"], orientation="up", beam_angle=25, tke_only=False
         )
-        tdat["tke"] = bnr.total_turbulent_kinetic_energy(
-            dat, noise=tdat["noise"], orientation="up", beam_angle=25
-        )
+        # Back in "inst" coordinate frame now
+        dat.velds.rotate2("beam")
+
+        # This needs U_mag in principal direction
         tdat["ti_noise"] = bnr.turbulence_intensity(
-            dat.velds.U_mag, detrend=False, noise=tdat["noise"]
+            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),
+            tdat["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),
+            tdat["U_mag"].isel(range=len(dat["range"]) // 2),
             freq_range=[0.2, 0.4],
             noise=tdat["noise"],
         )
         (
             tdat["dissipation_rate_SF"],
             tdat["noise_SF"],
             tdat["D_SF"],
-        ) = bnr.dissipation_rate_SF(dat.vel.isel(dir=2), r_range=[1, 5])
-        tdat["friction_vel"] = bnr.friction_velocity(
-            tdat, upwp_=tdat["stress_vec5"].sel(tau="upwp_"), z_inds=slice(1, 5), H=50
-        )
+        ) = bnr.dissipation_rate_SF(dat["vel"].isel(dir=2), r_range=[1, 5])
+
         slope_check = bnr.check_turbulence_cascade_slope(
             tdat["psd"].mean("time"), freq_range=[0.4, 4]
         )
+        # Check noise subtraction in psd function
         tdat["psd_noise"] = bnr.power_spectral_density(
-            dat["vel"].isel(dir=2, range=len(dat.range) // 2),
+            dat["vel"].isel(dir=2, range=len(dat["range"]) // 2),
             freq_units="Hz",
             noise=0.01,
         )
+        tdat["friction_vel"] = bnr.friction_velocity(
+            tdat, upwp_=tdat["stress_vec5"].sel(tau="upwp_"), z_inds=slice(1, 5), H=50
+        )
 
         if make_data:
             save(tdat, "Sig1000_tidal_bin.nc")