Created turbulence.py

Created a place for turbulence related calculations to live

Removed TKE from kinematics.py in calc

Added calculations for getting perturbation time series and turbulence kinetic energy

Added tests for tke, perturbation time series

Created docs for turbulence.py

docs/source/api/calc.rst

@@ -11,3 +11,4 @@ Calculations
    basic
    kinematics
    thermo
+   turbulence

docs/source/api/turbulence.rst

@@ -0,0 +1,11 @@
+***********************************
+Turbulence Time Series Calculations
+***********************************
+
+
+:mod:`metpy.calc.turbulence`
+============================
+
+.. automodule:: metpy.calc.turbulence
+   :members:
+   :undoc-members:

metpy/calc/__init__.py

@@ -1,6 +1,7 @@
 from .basic import *  # noqa
 from .kinematics import *  # noqa
 from .thermo import *  # noqa
+from . import turbulence  # noqa
 __all__ = []
 __all__.extend(basic.__all__)  # pylint: disable=undefined-variable
 __all__.extend(kinematics.__all__)  # pylint: disable=undefined-variable

metpy/calc/kinematics.py

@@ -198,29 +198,3 @@ def geostrophic_wind(heights, f, dx, dy, geopotential=False):
     grad = np.gradient(heights, *deltas)
     dx, dy = grad[0], grad[1]  # This throws away unused gradient components
     return -norm_factor * dy, norm_factor * dx
-
-
-@exporter.export
-def tke(u, v, w):
-    r'''Compute the turbulence kinetic energy (tke) from the time series of the
-    velocity components.
-
-    Parameters
-    ----------
-    u : array_like
-        The wind component along the x-axis
-    v : array_like
-        The wind component along the y-axis
-    w : array_like
-        The wind componennt along the z-axis
-
-    Returns
-    -------
-    array_like
-        The corresponding tke value(s)
-    '''
-
-    up = u - u.mean()
-    vp = v - v.mean()
-    wp = w - w.mean()
-    return np.sqrt(up * up + vp * vp + wp * wp)

metpy/calc/tests/test_turbulence.py

@@ -0,0 +1,134 @@
+import numpy as np
+from numpy.testing import (TestCase, assert_array_equal)
+from metpy.calc.turbulence import *  # noqa
+
+
+class TestTurbulenceKineticEnergy(TestCase):
+    def get_uvw_and_known_tke(self):
+        u = np.array([-2, -1, 0, 1, 2])
+        v = -u
+        w = 2 * u
+        #  0.5 * sqrt(2 + 2 + 8)
+        e_true = np.sqrt(12) / 2.
+        return u, v, w, e_true
+
+    def test_no_tke_1d(self):
+        observations = 5
+        # given all the values are the same, there should not be any tke
+        u = np.ones(observations)
+        v = np.ones(observations)
+        w = np.ones(observations)
+        e_zero = 0
+        assert_array_equal(e_zero, tke(u, v, w))
+
+    def test_no_tke_2d_axis_last(self):
+        observations = 5
+        instruments = 2
+        # given all the values are the same, there should not be any tke
+        u = np.ones((instruments, observations))
+        v = np.ones((instruments, observations))
+        w = np.ones((instruments, observations))
+        e_zero = np.zeros(instruments)
+        assert_array_equal(e_zero, tke(u, v, w, axis=-1))
+
+    def test_no_tke_2d_axis_first(self):
+        observations = 5
+        instruments = 2
+        # given all the values are the same, there should not be any tke
+        u = np.ones((observations, instruments))
+        v = np.ones((observations, instruments))
+        w = np.ones((observations, instruments))
+        e_zero = np.zeros(instruments)
+        assert_array_equal(e_zero, tke(u, v, w, axis=0))
+
+    def test_known_tke(self):
+        u, v, w, e_true = self.get_uvw_and_known_tke()
+        assert_array_equal(e_true, tke(u, v, w))
+
+    def test_known_tke_2d_axis_last(self):
+        '''test array with shape (3, 5) [pretend time axis is -1]'''
+        u, v, w, e_true = self.get_uvw_and_known_tke()
+        u = np.array([u, u, u])
+        v = np.array([v, v, v])
+        w = np.array([w, w, w])
+        e_true = e_true * np.ones(3)
+        assert_array_equal(e_true, tke(u, v, w, axis=-1))
+
+    def test_known_tke_2d_axis_first(self):
+        '''test array with shape (5, 3) [pretend time axis is 0]'''
+        u, v, w, e_true = self.get_uvw_and_known_tke()
+        u = np.array([u, u, u]).transpose()
+        v = np.array([v, v, v]).transpose()
+        w = np.array([w, w, w]).transpose()
+        e_true = e_true * np.ones(3).transpose()
+        assert_array_equal(e_true, tke(u, v, w, axis=0))
+        assert_array_equal(e_true, tke(u, v, w, axis=0, perturbation=True))
+
+
+class TestGetPerturbation(TestCase):
+    def get_pert_from_zero_mean(self):
+        ts = np.array([-2, -1, 0, 1, 2])
+        #  ts.meam() = 0
+        pert_true = ts.copy()
+        return ts, pert_true
+
+    def get_pert_from_non_zero_mean(self):
+        ts = np.array([-2, 0, 2, 4, 6])
+        # ts.mean() = 2
+        pert_true = np.array([-4, -2, 0, 2, 4])
+        return ts, pert_true
+
+    def test_no_perturbation_1d(self):
+        observations = 5
+        # given all the values are the same, there should not be perturbations
+        ts = np.ones(observations)
+        pert_zero = 0
+        assert_array_equal(pert_zero, get_perturbation(ts))
+
+    def test_no_perturbation_2d_axis_last(self):
+        observations = 5
+        instruments = 2
+        # given all the values are the same, there should not be perturbations
+        ts = np.ones((instruments, observations))
+        pert_zero = np.zeros((instruments, observations))
+        assert_array_equal(pert_zero, get_perturbation(ts, axis=-1))
+
+    def test_no_tke_2d_axis_first(self):
+        observations = 5
+        instruments = 2
+        # given all the values are the same, there should not be perturbations
+        ts = np.ones((observations, instruments))
+        pert_zero = np.zeros((observations, instruments))
+        assert_array_equal(pert_zero, get_perturbation(ts, axis=0))
+
+    def test_known_perturbation_zero_mean_1d(self):
+        ts, pert_known = self.get_pert_from_zero_mean()
+        assert_array_equal(pert_known, get_perturbation(ts))
+
+    def test_known_perturbation_zero_mean_2d_axis_last(self):
+        ts, pert_known = self.get_pert_from_zero_mean()
+        ts = np.array([ts, ts, ts])
+        pert_known = np.array([pert_known, pert_known, pert_known])
+        assert_array_equal(pert_known, get_perturbation(ts, axis=-1))
+
+    def test_known_perturbation_zero_mean_2d_axis_first(self):
+        ts, pert_known = self.get_pert_from_zero_mean()
+        ts = np.array([ts, ts, ts]).transpose()
+        pert_known = np.array([pert_known, pert_known, pert_known]).transpose()
+        assert_array_equal(pert_known, get_perturbation(ts, axis=0))
+
+    def test_known_perturbation_non_zero_mean_1d(self):
+        ts, pert_known = self.get_pert_from_non_zero_mean()
+        assert_array_equal(pert_known, get_perturbation(ts))
+
+    def test_known_perturbation_non_zero_mean_2d_axis_last(self):
+        ts, pert_known = self.get_pert_from_non_zero_mean()
+        ts = np.array([ts, ts, ts])
+        pert_known = np.array([pert_known, pert_known, pert_known])
+        assert_array_equal(pert_known, get_perturbation(ts, axis=-1))
+
+    def test_known_perturbation_non_zero_mean_2d_axis_first(self):
+        ts, pert_known = self.get_pert_from_non_zero_mean()
+        ts = np.array([ts, ts, ts]).transpose()
+        pert_known = np.array([pert_known, pert_known, pert_known]).transpose()
+        assert_array_equal(pert_known, get_perturbation(ts, axis=0))

metpy/calc/turbulence.py

@@ -0,0 +1,106 @@
+import numpy as np
+from ..package_tools import Exporter
+
+exporter = Exporter(globals())
+
+
+@exporter.export
+def get_perturbation(ts, axis=-1):
+    r"""Compute the perturbation from the mean of a time series.
+
+    Parameters
+    ----------
+    ts : array_like
+         The time series from which you wish to find the perturbation
+         time series (perturbation from the mean).
+
+    Returns
+    -------
+    array_like
+        The perturbation time series.
+
+    Other Parameters
+    ----------------
+    axis : int
+           The index of the time axis. Default is -1
+
+    Notes
+    -----
+    The perturbation time series produced by this funcation is defined as
+    the perturbations about the mean:
+
+    .. math:: x(t)^{\prime} = x(t) - \overline{x(t)}
+
+    """
+    slices = [slice(None)] * ts.ndim
+    slices[axis] = None
+    return ts - ts.mean(axis=axis)[slices]
+
+
+@exporter.export
+def tke(u, v, w, perturbation=False, axis=-1):
+    r"""Compute the turbulence kinetic energy (e) from the time series of the
+    velocity components.
+
+    Parameters
+    ----------
+    u : array_like
+        The wind component along the x-axis
+    v : array_like
+        The wind component along the y-axis
+    w : array_like
+        The wind component along the z-axis
+
+    perturbation : {False, True}, optional
+                   True if the u, v, and w components of wind speed supplied to
+                   the function are perturbation velocities. If False,
+                   perturbation velocities will be calculated by removing
+                   the mean value from each component.
+
+    Returns
+    -------
+    array_like
+        The corresponding turbulence kinetic energy value
+
+    Other Parameters
+    ----------------
+    axis : int
+           The index of the time axis. Default is -1
+
+    See Also
+    --------
+    get_perturbation : Used to compute perturbations if `perturbation`
+                       is False.
+
+    Notes
+    -----
+    Turbulence Kinetic Energy is computed as:
+
+    .. math:: e = 0.5 \sqrt{\overline{u^{\prime2}} +
+                            \overline{v^{\prime2}} +
+                            \overline{w^{\prime2}}},
+
+    where the velocity components
+
+    .. math:: u^{\prime}, v^{\prime}, u^{\prime}
+
+    are perturbation velocities. Fore more information on the subject, please
+    see [1]_.
+
+    References
+    ----------
+    .. [1] Garratt, J.R., 1994: The Atmospheric Boundary Layer. Cambridge
+           University Press, 316 pp.
+
+    """
+
+    if not perturbation:
+        u = get_perturbation(u, axis=axis)
+        v = get_perturbation(v, axis=axis)
+        w = get_perturbation(w, axis=axis)
+
+    u_cont = np.mean(u * u, axis=axis)
+    v_cont = np.mean(v * v, axis=axis)
+    w_cont = np.mean(w * w, axis=axis)
+
+    return 0.5 * np.sqrt(u_cont + v_cont + w_cont)