shear.py + tests

wetb/wind/shear.py

+'''
+Created on 16/06/2014
+
+@author: MMPE
+'''
+
+from scipy.optimize.optimize import fmin
+import numpy as np
+
+
+def _z_u(z_u_lst):
+    z = np.array([z for z, _ in z_u_lst])
+    u = np.array([np.mean(np.array([u])[:]) for _, u in z_u_lst])
+    return z, u
+
+def power_shear(alpha, z_ref, u_ref, z):
+    """Power shear
+
+    Parameters
+    ----------
+    alpha : int or float
+        The alpha shear parameter
+    z_ref : int or float
+        The reference height
+    z : int, float or array_like
+        The heights of interest
+    u_ref : int or float
+        The wind speed of the reference height
+
+    Returns
+    -------
+    power_shear : array-like
+        Wind speeds at height z
+
+    Excample
+    --------
+    >>> power_shear(.5, 70, [20,50,70], 9)
+    [ 4.81070235  7.60638829  9.        ]
+    """
+    return u_ref * (np.array(z) / z_ref) ** alpha
+
+
+def fit_power_shear(z_u_lst):
+    """Estimate power shear parameter, alpha, from a mean wind of two heights
+
+    Parameters
+    ----------
+    z_u_lst : [(z_ref, u_z_ref), (z1, u_z1),...]
+        - z_ref: Reference height\n
+        - u_z_ref: Wind speeds or mean wind speed at z_ref
+        - z1: another height
+        - u_z1: Wind speeds or mean wind speeds at z1
+
+    Returns
+    -------
+    alpha : float
+        power shear parameter
+
+    Excample
+    --------
+    >>> fit_power_shear([(85, 9.0), (21, 4.47118)])
+    0.50036320835
+    """
+    z, u = _z_u(z_u_lst)
+    z_hub, u_hub = z[0], u[0]
+    alpha, _ = np.polyfit(np.log(z / z_hub), np.log(u / u_hub), 1)
+    return alpha
+
+def fit_power_shear_ref(z_u_lst, z_ref):
+    """Estimate power shear parameter, alpha, for specified reference height
+
+    Parameters
+    ----------
+    z_u_lst : [(z1, u_z1), (z2, u_z2),...]
+        - z1: Some height
+        - z1_ref: Wind speeds or mean wind speed at z1
+        - z2: another height
+        - u_z2: Wind speeds or mean wind speeds at z1
+    z_ref : float or int
+        Reference height
+
+    Returns
+    -------
+    alpha : float
+        power shear parameter
+    u_ref : float
+        Wind speed at reference height
+
+    Excample
+    --------
+    >>> fit_power_shear_ref([(85, 8.88131), (21, 4.41832)],  87.13333)
+    [ 0.49938238  8.99192568]
+    """
+    def shear_error(x, z_u_lst, z_ref):
+        alpha, u_ref = x
+        return np.sum([(np.mean(u) - u_ref * (z / z_ref) ** alpha) ** 2 for z, u in z_u_lst])
+    return fmin(shear_error, (.5, 10), (z_u_lst, z_ref), disp=False)
+
+
+
+def log_shear(u_star, z0, z, Obukhov_length=None):
+    """logarithmic shear
+
+    Parameters
+    ----------
+    u_star : int or float
+        Friction velocity
+    z0 : int or float
+        Surface roughness [m]
+    z : int, float or array_like
+        The heights of interest
+
+    Returns
+    -------
+    log_shear : array-like
+        Wind speeds at height z
+
+    Excample
+    --------
+    >>> log_shear(1, 10, [20, 50, 70])
+    [ 1.73286795  4.02359478  4.86477537]
+    """
+    K = 0.4  # von Karmans constant
+    if Obukhov_length is None:
+        return u_star / K * (np.log(np.array(z) / z0))
+    else:
+        return u_star / K * (np.log(z / z0) - stability_term(z, z0, Obukhov_length))
+
+def stability_term(z, z0, L):
+    """Calculate the stability term for the log shear
+
+    Not validated!!!
+    """
+    zL = z / L
+    phi_us = lambda zL : (1 + 16 * np.abs(zL)) ** (-1 / 4)  #unstable
+    phi_s = lambda zL : 1 + 5 * zL  # stable
+    phi = np.zeros_like(zL) + np.nan
+
+    for m, f in [(((-2 <= zL) & (zL <= 0)), phi_us), (((0 < zL) & (zL <= 1)), phi_s)]:
+        phi[m] = f(zL[m])
+    f = L / z * (1 - phi)
+    psi = np.cumsum(np.r_[0, (f[1:] + f[:-1]) / 2 * np.diff(z / L)])
+    return psi
+
+
+def fit_log_shear(z_u_lst):
+    """Estimate log shear parameter, u_star and z0
+
+    Parameters
+    ----------
+    z_u_lst : [(z1, u_z1), (z2, u_z2),...]
+        - z1: Some height
+        - z1_ref: Wind speeds or mean wind speed at z1
+        - z2: another height
+        - u_z2: Wind speeds or mean wind speeds at z1
+
+    Returns
+    -------
+    u_star : float
+        friction velocity
+    z0 : float
+        Surface roughness [m]
+
+    Excample
+    --------
+    >>> fit_log_shear([(85, 8.88131), (21, 4.41832)],  87.13333)
+    [ 0.49938238  8.99192568]
+    """
+    def shear_error(x, z_u_lst):
+        u_star, z0 = x
+        return np.sum([(np.mean(u) - log_shear(u_star, z0, z)) ** 2 for z, u in z_u_lst])
+    return fmin(shear_error, (1, 1), (z_u_lst,), disp=False)
+
+

wetb/wind/tests/test_Shear.py

+'''
+Created on 05/06/2012
+
+@author: Mads
+'''
+import os
+from wetb.functions.geometry import xyz2uvw
+import wetb.gtsdf
+from wetb.wind.shear import power_shear, fit_power_shear, fit_power_shear_ref, \
+    log_shear, fit_log_shear, stability_term
+from pylab import *
+
+
+import unittest
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+
+all = True
+class TestMannTurbulence(unittest.TestCase):
+
+    def setUp(self):
+        self.testfilepath = "test_files/"
+
+        """
+        Sensor list of hdf5 files
+        0: WSP gl. coo.,Abs_vhor 85
+        1: WSP gl. coo.,Vdir_hor 85
+        2: WSP gl. coo.,Abs_vhor 53
+        3: WSP gl. coo.,Vdir_hor 53
+        4: WSP gl. coo.,Abs_vhor 21
+        5: WSP gl. coo.,Vdir_hor 21
+        6: WSP gl. coo.,Vx 85
+        7: WSP gl. coo.,Vy 85
+        8: WSP gl. coo.,Vz 85
+        9: WSP gl. coo.,Vx 53
+        10: WSP gl. coo.,Vy 53
+        11: WSP gl. coo.,Vz 53
+        12: WSP gl. coo.,Vx 21
+        13: WSP gl. coo.,Vy 21
+        14: WSP gl. coo.,Vz 21
+        """
+
+
+    def test_power_shear(self):
+        if all:
+            _, data, _ = wetb.gtsdf.load(self.testfilepath + 'wind3.hdf5')
+            u20 = data[:, 4]
+            u70 = data[:, 6]
+            z_u_lst1 = [(70, u70), (20, u20)]
+            z_u_lst2 = [ (80, 6.93121004105), (100, 7.09558010101), (60, 6.58919000626), (40, 6.12744998932)]
+
+            alpha = fit_power_shear(z_u_lst2)
+            self.assertAlmostEqual(alpha, 0.163432280539)
+            if 0:
+                z = np.arange(10, 100)
+                plt.plot([np.mean(u20), np.mean(u70)], [20, 70], 'bo')
+                alpha = fit_power_shear(z_u_lst1)
+                plt.plot(power_shear(alpha, 70, z, np.mean(u70)), z, 'b')
+
+                z, u = np.array(z_u_lst2).T
+                plt.plot(u, z, 'ro')
+                alpha = fit_power_shear(z_u_lst2)
+                z = np.sort(z)
+                plt.plot(power_shear(alpha, 80, z, 6.9), z, 'r')
+                plt.show()
+
+
+
+
+    def test_fit_power_shear1(self):
+
+        time, data, info = wetb.gtsdf.load(self.testfilepath + 'shear_noturb_85.hdf5')  #shear, alpha = 0.5
+
+        wsp85 = data[:, 0]
+        wsp53 = data[:, 2]
+        wsp21 = data[:, 4]
+        u85 = data[:, 7]
+        u21 = data[:, 13]
+
+        self.assertAlmostEqual(fit_power_shear([(85, wsp85), (53, wsp53), (21, wsp21)]), 0.5, delta=0.001)
+        self.assertAlmostEqual(fit_power_shear([(85, wsp85), (21, wsp21)]), 0.5, delta=0.001)
+        self.assertAlmostEqual(fit_power_shear([(85, u85), (21, u21)]), 0.5, delta=0.001)
+
+    def test_fit_power_shear_ref(self):
+
+        time, data, info = wetb.gtsdf.load(self.testfilepath + 'shear_noturb.hdf5')  #shear, alpha = 0.5
+
+        wsp85 = data[:, 0]
+        wsp53 = data[:, 2]
+        wsp21 = data[:, 4]
+
+
+        alpha, u_ref = fit_power_shear_ref([(85, wsp85), (21, wsp21)], 87.13333)
+        self.assertAlmostEqual(alpha, .5, delta=.001)
+        self.assertAlmostEqual(u_ref, 9, delta=.01)
+
+        alpha, u_ref = fit_power_shear_ref([(85, wsp85), (53, wsp53), (21, wsp21)], 87.13333)
+        self.assertAlmostEqual(alpha, .5, delta=.001)
+        self.assertAlmostEqual(u_ref, 9, delta=.01)
+
+
+    def test_fit_power_shear3(self):
+
+        time, data, info = wetb.gtsdf.load(self.testfilepath + 'shear.hdf5')  #shear, alpha = 0.5
+        #print "\n".join(["%d: %s" % (i, n) for i, n in enumerate(info['attribute_names'])])
+        """
+        0: WSP gl. coo.,Abs_vhor 85
+        1: WSP gl. coo.,Vdir_hor 85
+        2: WSP gl. coo.,Abs_vhor 53
+        3: WSP gl. coo.,Vdir_hor 53
+        4: WSP gl. coo.,Abs_vhor 21
+        5: WSP gl. coo.,Vdir_hor 21
+        6: WSP gl. coo.,Vx 85
+        7: WSP gl. coo.,Vy 85
+        8: WSP gl. coo.,Vz 85
+        9: WSP gl. coo.,Vx 53
+        10: WSP gl. coo.,Vy 53
+        11: WSP gl. coo.,Vz 53
+        12: WSP gl. coo.,Vx 21
+        13: WSP gl. coo.,Vy 21
+        14: WSP gl. coo.,Vz 21
+        """
+        wsp85 = data[:, 0]
+        wsp53 = data[:, 2]
+        wsp21 = data[:, 4]
+        u85 = data[:, 7]
+        u21 = data[:, 13]
+
+        alpha, u_ref = fit_power_shear_ref([(85, wsp85), (21, wsp21)], 87.13333)
+        self.assertAlmostEqual(alpha, .5, delta=.2)
+        self.assertAlmostEqual(u_ref, 9, delta=.01)
+
+
+    def test_log_shear(self):
+        u = log_shear(2, 3, 9)
+        self.assertAlmostEqual(u, 5.49306144)
+
+    def test_fit_log_shear(self):
+        zu = [(85, 8.88131), (21, 4.41832)]
+        u_star, z0 = fit_log_shear(zu)
+        if 0:
+            for z, u in zu:
+                plt.plot(u, z, 'r.')
+        z = np.arange(10, 100)
+        for _zu, b in zip(zu, log_shear(u_star, z0, [85, 21])):
+            self.assertAlmostEqual(_zu[1], b, 4)
+
+
+    def test_show_log_shear(self):
+        if 0:
+            for ustar in [1, 2]:
+                for z0 in [1, 10]:
+                    z = np.arange(z0, 200)
+                    plot(log_shear(ustar, z0, z), z, label="z0=%d, u*=%d" % (z0, ustar))
+            yscale('log')
+            legend()
+            show()
+
+    def test_show_log_shear_stability(self):
+        if 0:
+            z0 = 1
+            ustar = 1
+            z = np.arange(z0, 200)
+            for L in [-2000, -100, 100, 2000]:
+                plot(log_shear(ustar, z0, z, L), z, label="L=%d" % (L))
+            #yscale('log')
+            legend()
+            show()
+
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()