wetb/wind/tests/test_Shear.py

@@ -3,13 +3,6 @@ Created on 05/06/2012
 
 @author: Mads
 '''
-from __future__ import unicode_literals
-from __future__ import print_function
-from __future__ import division
-from __future__ import absolute_import
-from builtins import zip
-from future import standard_library
-standard_library.install_aliases()
 import os
 from wetb.wind.utils import xyz2uvw
 import wetb.gtsdf

wetb/wind/tests/test_weibull.py

@@ -8,8 +8,8 @@ import matplotlib.pyplot as plt
 import numpy as np
 from wetb.wind import weibull
 
-class TestWeibull(unittest.TestCase):
 
+class TestWeibull(unittest.TestCase):
 
     def test_weibull(self):
         wb = weibull.pdf(4, 2)
@@ -21,37 +21,39 @@ class TestWeibull(unittest.TestCase):
 
     def test_random_weibull(self):
         y = weibull.random(4, 2, 1000000)
-        pdf, x = np.histogram(y, 100, normed=True)
+        pdf, x = np.histogram(y, 100, density=True)
         x = (x[1:] + x[:-1]) / 2
         self.assertLess(sum((pdf - weibull.pdf(4, 2)(x)) ** 2), 0.0001)
         if 0:
             plt.plot(x, pdf)
             plt.plot(x, weibull.pdf(4, 2)(x))
             plt.show()
- 
+
     def test_fit_weibull(self):
+        np.random.seed(1)
         y = weibull.random(4, 2, 1000000)
         A, k = weibull.fit(y)
         self.assertAlmostEqual(A, 4, delta=0.01)
         self.assertAlmostEqual(k, 2, delta=0.01)
         if 0:
-            plt.hist(y, 100, normed=True)
+            plt.hist(y, 100, density=True)
             x = np.arange(0, 20, .1)
             plt.plot(x, weibull.pdf(4, 2)(x))
             plt.show()
-             
+
     def test_weibull_cdf(self):
         wb = weibull.pdf(4, 2)
         x = np.arange(0, 20, .01)
-        cdf = weibull.cdf(4,2)
-        self.assertEqual(cdf(999),1)
-        np.testing.assert_array_almost_equal(np.cumsum(wb(x))/100, cdf(x), 3)
-         
+        cdf = weibull.cdf(4, 2)
+        self.assertEqual(cdf(999), 1)
+        np.testing.assert_array_almost_equal(np.cumsum(wb(x)) / 100, cdf(x), 3)
+
         if 0:
-            plt.plot(x, np.cumsum(wb(x))*.01)
+            plt.plot(x, np.cumsum(wb(x)) * .01)
             plt.plot(x, cdf(x))
             plt.show()
-             
+
+
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testweibull']
     unittest.main()

wetb/wind/turbulence/__init__.py

+

wetb/wind/turbulence/constraint_file.py

 import numpy as np
 import os
 from wetb.hawc2.htc_file import HTCFile
+
+
 class ConstraintFile(object):
-    
+
     def __init__(self, center_gl_xyz, box_transport_speed, no_grid_points=(4096, 32, 32), box_size=(6000, 100, 100)):
         """Generate list of constraints for turbulence constraint simulator
 
@@ -23,22 +25,21 @@ class ConstraintFile(object):
         self.box_size = box_size
         self.dxyz = [d / (n - 1) for d, n in zip(self.box_size, self.no_grid_points)]
 
-        self.constraints = {'u':[], 'v':[], 'w':[]}
-        
+        self.constraints = {'u': [], 'v': [], 'w': []}
+
     def load(self, filename):
         """Load constraint from existing constraint file (usable for plotting constraints via time_series()"""
         with open(filename) as fid:
             lines = fid.readlines()
         for l in lines:
             values = l.split(";")
-            mx,my,mz= map(int, values[:3])
-            comp = values[3]#{"100":'u','010':'v','001':'w'}["".join(values[3:6])]
-            self.constraints[comp].append((mx,my,mz,float(values[-1])))
-            
+            mx, my, mz = map(int, values[:3])
+            comp = values[3]  # {"100":'u','010':'v','001':'w'}["".join(values[3:6])]
+            self.constraints[comp].append((mx, my, mz, float(values[-1])))
 
     def add_constraints(self, glpos, tuvw, subtract_mean=True, fail_outside_box=True, nearest=True):
         """Add constraints to constraint file
-        
+
         parameters
         ----------
         glpos : (float or array_like, float or array_like, float or array_like)
@@ -54,7 +55,7 @@ class ConstraintFile(object):
             if False, mann coordinates modulo N is used
         nearest : boolean, optional
             if True, the wsp at the position closest to the mann grid point are applied\n
-            if False, the average of all the wsp, that map to the mann grid point are applied 
+            if False, the average of all the wsp, that map to the mann grid point are applied
         """
         nx, ny, nz = self.no_grid_points
         dx, dy, dz = self.dxyz
@@ -62,50 +63,53 @@ class ConstraintFile(object):
 
         time, u, v, w = (list(tuvw.T) + [None, None])[:4]
         x, y, z = glpos
-        mxs = np.round((time * self.box_transport_speed + (center_y - y)) / dx).astype(np.int)
-        if not np.array(y).shape==mxs.shape:
-            y = np.zeros_like(mxs)+y
-        if not np.array(x).shape==mxs.shape:
-            x = np.zeros_like(mxs)+x 
-        if not np.array(z).shape==mxs.shape:
-            z = np.zeros_like(mxs)+z
-        x_err = y-(-(mxs*dx-time*self.box_transport_speed-center_y))
-        mys = np.round(((ny - 1) / 2. + (-x+center_x) / dy)).astype(np.int)
-        mzs = np.round(((nz - 1) / 2. + (center_z - z) / dz)).astype(np.int)
+        mxs = np.round((time * self.box_transport_speed + (center_y - y)) / dx).astype(int)
+        if not np.array(y).shape == mxs.shape:
+            y = np.zeros_like(mxs) + y
+        if not np.array(x).shape == mxs.shape:
+            x = np.zeros_like(mxs) + x
+        if not np.array(z).shape == mxs.shape:
+            z = np.zeros_like(mxs) + z
+        x_err = y - (-(mxs * dx - time * self.box_transport_speed - center_y))
+        mys = np.round(((ny - 1) / 2. + (-x + center_x) / dy)).astype(int)
+        mzs = np.round(((nz - 1) / 2. + (center_z - z) / dz)).astype(int)
 
-        y_err = x-(-((mys-(ny-1)/2)*dy-center_x))
-        z_err = z-(-(mzs-(nz-1)/2)*dz+center_z)
-        pos_err = np.sqrt(x_err**2+y_err**2+z_err**2)
+        y_err = x - (-((mys - (ny - 1) / 2) * dy - center_x))
+        z_err = z - (-(mzs - (nz - 1) / 2) * dz + center_z)
+        pos_err = np.sqrt(x_err**2 + y_err**2 + z_err**2)
         if fail_outside_box:
-            for ms, n in zip([mxs,mys,mzs],self.no_grid_points):
+            for ms, n in zip([mxs, mys, mzs], self.no_grid_points):
                 if ms.min() + 1 < 1:
                     i = np.argmin(ms)
-                    raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz))
+                    raise ValueError(
+                        "At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" %
+                        (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx, ny, nz))
                 if ms.max() + 1 > n:
                     i = np.argmax(ms)
-                    raise ValueError("At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" % (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx , ny , nz))
+                    raise ValueError(
+                        "At time, t=%s, global position (%s,%s,%s) maps to grid point (%d,%d,%d) whichs is outside turbulence box, range (1..%d,1..%d,1..%d)" %
+                        (time[i], x[i], y[i], z[i], mxs[i] + 1, mys[i] + 1, mzs[i] + 1, nx, ny, nz))
         else:
-            mxs%=nx
-            mys%=ny
-            mzs%=nz
-        mann_pos = ["%06d%02d%02d"%(mx,my,mz) for mx,my,mz in zip(mxs,mys,mzs)]
+            mxs %= nx
+            mys %= ny
+            mzs %= nz
+        mann_pos = ["%06d%02d%02d" % (mx, my, mz) for mx, my, mz in zip(mxs, mys, mzs)]
         i = 0
         N = len(mann_pos)
-        uvw_lst = [uvw for uvw in [u,v,w] if uvw is not None]
+        uvw_lst = [uvw for uvw in [u, v, w] if uvw is not None]
         if subtract_mean:
             uvw_lst = [uvw - uvw.mean() for uvw in uvw_lst]
-        while i<N:
-            for j in range(i+1,N+1):
-                if j==N or mann_pos[i]!=mann_pos[j]:
+        while i < N:
+            for j in range(i + 1, N + 1):
+                if j == N or mann_pos[i] != mann_pos[j]:
                     break
             for comp, wsp in zip(['u', 'v', 'w'], uvw_lst):
                 if nearest:
                     value = wsp[i:j][np.argmin(pos_err[i:j])]
-                else: #mean
+                else:  # mean
                     value = wsp[i:j].mean()
                 self.constraints[comp].append((mxs[i] + 1, mys[i] + 1, mzs[i] + 1, value))
-            i=j
-                    
+            i = j
 
     def __str__(self):
         return "\n".join(["%d;%d;%d;%s;%.10f" % (mx, my, mz, comp, value) for comp in ['u', 'v', 'w']
@@ -117,15 +121,14 @@ class ConstraintFile(object):
             os.makedirs(path, exist_ok=True)
         with open(os.path.join(path, "%s.con" % name), 'w') as fid:
             fid.write(str(self))
-            
+
     def time_series(self, y_position=0):
-        time = (np.arange(self.no_grid_points[0])*float(self.dxyz[0]) + y_position)/ self.box_transport_speed
-        u,v,w = [(np.zeros_like(time)+ np.nan)]*3
-        for uvw, constr in zip([u,v,w], [np.array(self.constraints[uvw]) for uvw in 'uvw']):
-            if constr.shape[0]>0:
-                uvw[constr[:,0].astype(np.int)-1] =constr[:,3] 
-        return time, np.array([u,v,w]).T
-    
+        time = (np.arange(self.no_grid_points[0]) * float(self.dxyz[0]) + y_position) / self.box_transport_speed
+        u, v, w = [(np.zeros_like(time) + np.nan)] * 3
+        for uvw, constr in zip([u, v, w], [np.array(self.constraints[uvw]) for uvw in 'uvw']):
+            if constr.shape[0] > 0:
+                uvw[constr[:, 0].astype(int) - 1] = constr[:, 3]
+        return time, np.array([u, v, w]).T
 
     def simulation_cmd(self, ae23, L, Gamma, seed, name, folder="./constraints/"):
         assert isinstance(seed, int) and seed > 0, "seed must be a positive integer"
@@ -133,17 +136,13 @@ class ConstraintFile(object):
         cmd += "%.6f %.2f %.2f %d " % (ae23, L, Gamma, seed)
         cmd += "./turb/%s_s%s %s/%s.con" % (name, seed, folder, name)
         return cmd
-    
+
     def hawc2_mann_section(self, name, seed):
         htc = HTCFile()
         mann = htc.wind.add_section("mann")
         for uvw in 'uvw':
-            setattr(mann, 'filename_%s'%uvw, "./turb/%s_s%04d%s.bin" % ( name, int(seed),uvw))
-        for uvw, l,n in zip('uvw', self.box_size, self.no_grid_points):
-            setattr(mann, 'box_dim_%s'%uvw, [n,l/n])
+            setattr(mann, 'filename_%s' % uvw, "./turb/%s_s%04d%s.bin" % (name, int(seed), uvw))
+        for uvw, l, n in zip('uvw', self.box_size, self.no_grid_points):
+            setattr(mann, 'box_dim_%s' % uvw, [n, l / n])
         mann.dont_scale = 1
         return mann
-        
-        
-
-

wetb/wind/turbulence/mann_parameters.py

@@ -14,8 +14,6 @@ from wetb.wind.turbulence.spectra import spectra, logbin_spectra, plot_spectra,\
     detrend_wsp
 
 
-
-
 sp1, sp2, sp3, sp4 = np.load(os.path.dirname(__file__).replace("library.zip", '') + "/mann_spectra_data.npy")
 
 yp = np.arange(-3, 3.1, 0.1)
@@ -26,8 +24,7 @@ RBS3 = RectBivariateSpline(xp, yp, sp3)
 RBS4 = RectBivariateSpline(xp, yp, sp4)
 
 
-
-#def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
+# def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
 #    if isinstance(fs, (int, float)):
 #        fs = [fs] * len(u_lst)
 #    if v_lst is None:
@@ -41,10 +38,6 @@ RBS4 = RectBivariateSpline(xp, yp, sp4)
 #        return [np.mean(np.array([kuvw[i] for kuvw in uvw_lst]), 0) for i in range(5)]
 
 
-
-
-
-
 def get_mann_model_spectra(ae, L, G, k1):
     """Mann model spectra
 
@@ -73,15 +66,11 @@ def get_mann_model_spectra(ae, L, G, k1):
     xq = np.log10(L * k1)
     yq = (np.zeros_like(xq) + G)
     f = L ** (5 / 3) * ae
-    uu = f * RBS1.ev(yq, xq)
-    vv = f * RBS2.ev(yq, xq)
-    ww = f * RBS3.ev(yq, xq)
-    uw = f * RBS4.ev(yq, xq)
+    m = (yq >= 0) & (yq <= 5) & (xq >= -3) & (xq <= 3)
+    uu, vv, ww, uw = [f * np.where(m, RBS.ev(yq, xq), np.nan) for RBS in [RBS1, RBS2, RBS3, RBS4]]
     return uu, vv, ww, uw
 
 
-
-
 def _local_error(x, k1, uu, vv, ww=None, uw=None):
 
     ae, L, G = x
@@ -95,7 +84,9 @@ def _local_error(x, k1, uu, vv, ww=None, uw=None):
             val += np.sum((k1 * ww - k1 * tmpww) ** 2) + np.sum((k1 * uw - k1 * tmpuw) ** 2)
     return val
 
-def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2, min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
+
+def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
+                           min_bin_count=2, start_vals_for_optimisation=(0.01, 50, 3.3), plt=False):
     """Fit a mann model to the spectra
 
     Bins the spectra, into logarithmic sized bins and find the mann model parameters,
@@ -139,7 +130,8 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
     >>> ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v))
     """
     from scipy.optimize import fmin
-    x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False)
+    x = fmin(_local_error, start_vals_for_optimisation, logbin_spectra(
+        k1, uu, vv, ww, uw, log10_bin_size, min_bin_count), disp=False)
 
     if plt:
         if not hasattr(plt, 'plot'):
@@ -147,7 +139,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
 #         plot_spectra(k1, uu, vv, ww, uw, plt=plt)
 #         plot_mann_spectra(*x, plt=plt)
         ae, L, G = x
-        plot_fit(ae, L, G, k1, uu, vv, ww, uw,  log10_bin_size=log10_bin_size, plt=plt)
+        plot_fit(ae, L, G, k1, uu, vv, ww, uw, log10_bin_size=log10_bin_size, plt=plt)
         plt.title('ae:%.3f, L:%.1f, G:%.2f' % tuple(x))
         plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
         plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
@@ -155,6 +147,7 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
         plt.show()
     return x
 
+
 def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
     """Fit a mann model to the spectra
 
@@ -198,40 +191,59 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
     return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
 
 
-def var2ae(variance, spatial_resolution, N, L, G):
+def var2ae(variance, L, G, U, T=600, sample_frq=10, plt=False):
     """Fit alpha-epsilon to match variance of time series
 
     Parameters
     ----------
     variance : array-like
         variance of u vind component
-    spatial_resolution : int, float or array_like
-        Distance between samples in meters
-        - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points, 
-        Nx is number of points and Lx is box length in meters
-        - For time series: Sample frequency / U
-    N : int
     L : int or float
         Length scale of Mann model
     G : int or float
         Gamma of Mann model
+    U : int or float
+        Mean wind speed
+    T: int or float
+        Length [s] of signal, from which the variance is calculated
+    sample_frq: int or float
+        Sample frequency [Hz] of signal from which the variance is calculated
 
     Returns
     -------
     ae : float
-        Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
+        Alpha epsilon^(2/3) of Mann model that makes the energy of the model in the
+        frequency range [1/length, sample_frq] equal to the variance of u
     """
-    
-    k1 = np.logspace(1,10,1000)/100000000
+
+    k_low, k_high = 2 * np.pi / (U * np.array([T, 2 / sample_frq]))
+    k1 = 10 ** (np.linspace(np.log10(k_low), np.log10(k_high), 1000))
+
     def get_var(uu):
         return np.trapz(2 * uu[:], k1[:])
 
     v1 = get_var(get_mann_model_spectra(0.1, L, G, k1)[0])
     v2 = get_var(get_mann_model_spectra(0.2, L, G, k1)[0])
     ae = (variance - v1) / (v2 - v1) * .1 + .1
+    if plt is not False:
+        if not hasattr(plt, 'plot'):
+            import matplotlib.pyplot as plt
+        muu = get_mann_model_spectra(ae, L, G, k1)[0]
+        plt.semilogx(k1, k1 * muu, label='ae:%.3f, L:%.1f, G:%.2f' % (ae, L, G))
+        plt.legend()
+        plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
+        plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
     return ae
 
 
+def ae2ti(ae23, L, G, U, T=600, sample_frq=10):
+    k_low, k_high = 2 * np.pi / (U * np.array([T, 2 / sample_frq]))
+    k1 = 10 ** (np.linspace(np.log10(k_low), np.log10(k_high), 1000))
+
+    uu = get_mann_model_spectra(ae23, L, G, k1)[0]
+    var = np.trapz(2 * uu[:], k1[:])
+    return np.sqrt(var) / U
+
 
 def fit_ae(spatial_resolution, u, L, G, plt=False):
     """Fit alpha-epsilon to match variance of time series
@@ -239,8 +251,8 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
     Parameters
     ----------
     spatial_resolution : int, float or array_like
-        Distance between samples in meters
-        - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points, 
+        Number of points pr meterDistance between samples in meters
+        - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
         Nx is number of points and Lx is box length in meters
         - For time series: Sample frequency / U
     u : array-like
@@ -255,22 +267,22 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
     ae : float
         Alpha epsilon^(2/3) of Mann model that makes the energy of the model equal to the varians of u
     """
-    #if len(u.shape) == 1:
+    # if len(u.shape) == 1:
     #    u = u.reshape(len(u), 1)
 #     if min_bin_count is None:
 #         min_bin_count = max(2, 6 - u.shape[0] / 2)
 #     min_bin_count = 1
     def get_var(k1, uu):
-        l = 0  #128 // np.sqrt(u.shape[1])
-        return np.trapz(2 * uu[l:], k1[l:])
+        l = 0  # 128 // np.sqrt(u.shape[1])
+        return np.mean(np.trapz(2 * uu[l:], k1[l:], axis=0))
 
     k1, uu = spectra(spatial_resolution, u)[:2]
-    v = get_var(k1,uu)
+    v = get_var(k1, uu)
     v1 = get_var(k1, get_mann_model_spectra(0.1, L, G, k1)[0])
     v2 = get_var(k1, get_mann_model_spectra(0.2, L, G, k1)[0])
     ae = (v - v1) / (v2 - v1) * .1 + .1
 #     print (ae)
-#     
+#
 #     k1 = spectra(sf, u)[0]
 #     v1 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.1, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
 #     v2 = get_var(*logbin_spectra(k1, get_mann_model_spectra(0.2, L, G, k1)[0], min_bin_count=min_bin_count)[:2])
@@ -283,7 +295,7 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
         if not hasattr(plt, 'plot'):
             import matplotlib.pyplot as plt
         plt.semilogx(k1, k1 * uu, 'b-', label='uu')
-        k1_lb, uu_lb = logbin_spectra(*spectra(sf, u), min_bin_count=1)[:2]
+        k1_lb, uu_lb = logbin_spectra(*spectra(spatial_resolution, u), min_bin_count=1)[:2]
 
         plt.semilogx(k1_lb, k1_lb * uu_lb, 'r--', label='uu_logbin')
         muu = get_mann_model_spectra(ae, L, G, k1)[0]
@@ -296,24 +308,27 @@ def fit_ae(spatial_resolution, u, L, G, plt=False):
 
 
 def plot_fit(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
-#    if plt is None:
-#        import matplotlib.pyplot as plt
+    #    if plt is None:
+    #        import matplotlib.pyplot as plt
     plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, plt)
     plot_mann_spectra(ae, L, G, "-", mean_u, plt)
 
 
-
 def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
     if plt is None:
         import matplotlib.pyplot as plt
-    mf = 10 ** (np.linspace(-4, 1, 1000))
+    mf = 10 ** (np.linspace(-4, 3, 1000))
+    mf = 10 ** (np.linspace(-4, 3, 1000000))
     muu, mvv, mww, muw = get_mann_model_spectra(ae, L, G, mf)
-    plt.title("ae: %.3f, L: %.2f, G:%.2f"%(ae,L,G))
-    if 'uu' in spectra: plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
-    if 'vv' in spectra: plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
-    if 'ww' in spectra: plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
-    if 'uw' in spectra: plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
-
+    plt.title("ae: %.3f, L: %.2f, G:%.2f" % (ae, L, G))
+    if 'uu' in spectra:
+        plt.semilogx(mf, mf * muu * 10 ** 0 / u_ref ** 2, 'r' + style)
+    if 'vv' in spectra:
+        plt.semilogx(mf, mf * mvv * 10 ** 0 / u_ref ** 2, 'g' + style)
+    if 'ww' in spectra:
+        plt.semilogx(mf, mf * mww * 10 ** 0 / u_ref ** 2, 'b' + style)
+    if 'uw' in spectra:
+        plt.semilogx(mf, mf * muw * 10 ** 0 / u_ref ** 2, 'm' + style)
 
 
 if __name__ == "__main__":
@@ -323,30 +338,29 @@ if __name__ == "__main__":
     import matplotlib.pyplot as plt
 
     """Example of fitting Mann parameters to a time series"""
-    ds = gtsdf.Dataset(os.path.dirname(wind.__file__)+"/tests/test_files/WspDataset.hdf5")#'unit_test/test_files/wspdataset.hdf5')
+    ds = gtsdf.Dataset(os.path.dirname(wind.__file__) +
+                       "/tests/test_files/WspDataset.hdf5")  # 'unit_test/test_files/wspdataset.hdf5')
     f = 35
     u, v = wsp_dir2uv(ds.Vhub_85m, ds.Dir_hub_)
- 
+
     u_ref = np.mean(u)
     u -= u_ref
- 
+
     sf = f / u_ref
-    ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt = None)
-    print (ae, L, G)
-    print (u.shape)
-    print (ds.Time[-1])
+    ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v), plt=None)
+    print(ae, L, G)
+    print(u.shape)
+    print(ds.Time[-1])
     plt.plot(u)
     plt.plot(detrend_wsp(u)[0])
     plt.show()
-    print (fit_ae(sf, detrend_wsp(u)[0], L,  G, plt))
-    print (var2ae(detrend_wsp(u)[0].var(), L,  G,))
-    print ()
-    print (fit_ae(sf, u[:21000], L,  G))
-    print (var2ae(u[:21000].var(), L,  G,))
-    
+    print(fit_ae(sf, detrend_wsp(u)[0], L, G, plt))
+    print(var2ae(detrend_wsp(u)[0].var(), L, G,))
+    print()
+    print(fit_ae(sf, u[:21000], L, G))
+    print(var2ae(u[:21000].var(), L, G,))
 
 
- 
 #     """Example of fitting Mann parameters to a "series" of a turbulence box"""
 #     l = 16384
 #     nx = 8192
@@ -359,7 +373,6 @@ if __name__ == "__main__":
 #     print (ae, L, G)
 
 
-
 #     """Example of fitting Mann parameters to a "series" of a turbulence box"""
 #     l = 4778.3936
 #     nx = 8192
@@ -373,4 +386,4 @@ if __name__ == "__main__":
 #         #ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt)
 #         #print (fit_ae(sf, u, 73.0730383576,  2.01636095317))
 #         print (u.std())
-#     #print (ae, L, G)
\ No newline at end of file
+#     #print (ae, L, G)

wetb/wind/turbulence/mann_turbulence.py

@@ -8,38 +8,42 @@ import numpy as np
 from wetb.wind.turbulence.spectra import spectra, spectra_from_time_series
 name_format = "mann_l%.1f_ae%.4f_g%.1f_h%d_%dx%dx%d_%.3fx%.2fx%.2f_s%04d%c.turb"
 
-def load(filename, N=(32,32)):
+
+def load(filename, N=(32, 32)):
     """Load mann turbulence box
-    
+
     Parameters
     ----------
     filename : str
         Filename of turbulence box
     N : tuple, (ny,nz) or (nx,ny,nz)
         Number of grid points
-        
+
     Returns
     -------
     turbulence_box : nd_array
-    
+
     Examples
     --------
-    >>> u = load('turb_u.dat') 
+    >>> u = load('turb_u.dat')
     """
     data = np.fromfile(filename, np.dtype('<f'), -1)
-    if len(N)==2:
-        ny,nz = N
-        nx = len(data)/(ny*nz)
-        assert nx==int(nx), "Size of turbulence box (%d) does not match ny x nz (%d), nx=%.2f" % (len(data),ny*nz, nx)
-        nx=int(nx)
+    if len(N) == 2:
+        ny, nz = N
+        nx = len(data) / (ny * nz)
+        assert nx == int(nx), "Size of turbulence box (%d) does not match ny x nz (%d), nx=%.2f" % (
+            len(data), ny * nz, nx)
+        nx = int(nx)
     else:
-        nx,ny,nz = N
-        assert len(data) == nx*ny*nz, "Size of turbulence box (%d) does not match nx x ny x nz (%d)" % (len(data),nx*ny*nz)
-    return data.reshape(nx , ny * nz)
+        nx, ny, nz = N
+        assert len(data) == nx * ny * \
+            nz, "Size of turbulence box (%d) does not match nx x ny x nz (%d)" % (len(data), nx * ny * nz)
+    return data.reshape(nx, ny * nz)
 
-def load_uvw(filenames, N=(1024,32,32)):
+
+def load_uvw(filenames, N=(1024, 32, 32)):
     """Load u, v and w turbulence boxes
-    
+
     Parameters
     ----------
     filenames : list or str
@@ -47,34 +51,40 @@ def load_uvw(filenames, N=(1024,32,32)):
         if str: filename pattern where u,v,w are replaced with '%s'
     N : tuple
         Number of grid point in the x, y and z direction
-    
+
     Returns
     -------
     u,v,w : list of np_array
-    
+
     Examples
     --------
-    >>> u,v,w =load_uvw('turb_%s.dat') 
+    >>> u,v,w =load_uvw('turb_%s.dat')
     """
-    if isinstance(filenames,str):
-        return [load(filenames%uvw, N) for uvw in 'uvw']
+    if isinstance(filenames, str):
+        return [load(filenames % uvw, N) for uvw in 'uvw']
     else:
         return [load(f, N) for f in filenames]
-    
-def parameters2name(no_grid_points,box_dimension,ae23,L, Gamma,high_frq_compensation, seed, folder="./turb/"):
-    
+
+
+def save(turb, filename):
+    np.array(turb).astype('<f').tofile(filename)
+
+
+def parameters2name(no_grid_points, box_dimension, ae23, L, Gamma, high_frq_compensation, seed, folder="./turb/"):
+
     dxyz = tuple(np.array(box_dimension) / no_grid_points)
-    return ["./turb/" + name_format % ((L, ae23, Gamma, high_frq_compensation) + no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]    
+    return ["./turb/" + name_format % ((L, ae23, Gamma, high_frq_compensation) +
+                                       no_grid_points + dxyz + (seed, uvw)) for uvw in ['u', 'v', 'w']]
+
 
-    
 def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
     """Fit mann parameters, ae, L, G
-    
+
     Parameters
     ----------
     spatial_resolution : inf, float or array_like
         Distance between samples in meters
-        - For turbulence boxes: 1/dx = Lx/Nx where dx is distance between points, 
+        - For turbulence boxes: 1/dx = Lx/Nx where dx is distance between points,
         Nx is number of points and Lx is box length in meters
         - For time series: Sample frequency / U
     u : array_like
@@ -89,7 +99,7 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
         The w-wind component
         - if shape is (r,): One time series with *r* observations\n
         - if shape is (r,c): *c* different time series with *r* observations\n
-    
+
     Returns
     -------
     ae : int or float
@@ -105,4 +115,4 @@ def fit_mann_parameters(spatial_resolution, u, v, w=None, plt=None):
 
 def fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=None):
     from wetb.wind.turbulence.mann_parameters import fit_mann_model_spectra
-    return fit_mann_model_spectra(*spectra_from_time_series(sample_frq, Uvw_lst), plt=plt)
\ No newline at end of file
+    return fit_mann_model_spectra(*spectra_from_time_series(sample_frq, Uvw_lst), plt=plt)

wetb/wind/turbulence/spectra.py

@@ -7,6 +7,8 @@ Created on 27/11/2015
 
 import numpy as np
 import warnings
+
+
 def spectrum(x, y=None, k=1):
     """PSD or Cross spectrum (only positive half)
 
@@ -25,6 +27,7 @@ def spectrum(x, y=None, k=1):
     """
     if x is None:
         return None
+
     def fft(x):
         return np.fft.fft(x.T).T / len(x)
 
@@ -34,22 +37,23 @@ def spectrum(x, y=None, k=1):
     else:
         fftx = fft(x) * np.conj(fft(y))  # Cross spectra
 
-    fftx = fftx[:len(fftx) // 2 ] * 2  # positive half * 2
+    fftx = fftx[:len(fftx) // 2] * 2  # positive half * 2
 
 #    if len(fftx.shape) == 2:
 #        fftx = np.mean(fftx, 1)
     with warnings.catch_warnings():
-        warnings.simplefilter("ignore")  
+        warnings.simplefilter("ignore")
         return np.real(fftx * len(x) / (2 * k))[1:]
 
+
 def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
     """Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
 
     Parameters
     ----------
     spatial_resolution : int, float or array_like
-        Distance between samples in meters
-        - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points, 
+        Sample points per meter
+        - For turbulence boxes: 1/dx = Nx/Lx where dx is distance between points,
         Nx is number of points and Lx is box length in meters
         - For time series: Sample frequency / U
     u : array_like
@@ -81,24 +85,30 @@ def spectra(spatial_resolution, u, v=None, w=None, detrend=True):
     assert isinstance(spatial_resolution, (int, float))
 
     k = 2 * np.pi * spatial_resolution
-    if v is not None: assert u.shape==v.shape
-    if w is not None: assert u.shape==w.shape
-    
-    if 1 and len(u.shape)==2:
-        assert np.abs(np.mean(u,0)).max()<1
-        if v is not None: assert np.abs(np.mean(v,0)).max()<1
-        if w is not None: assert np.abs(np.mean(w,0)).max()<1
+    if v is not None:
+        assert u.shape == v.shape
+    if w is not None:
+        assert u.shape == w.shape
+
+    if 1 and len(u.shape) == 2:
+        #         assert np.abs(np.mean(u, 0)).max() < 2
+        #         if v is not None:
+        #             assert np.abs(np.mean(v, 0)).max() < 1
+        #         if w is not None:
+        #             assert np.abs(np.mean(w, 0)).max() < 1
         if isinstance(k, float):
-            k = np.repeat(k,u.shape[1])
+            k = np.repeat(k, u.shape[1])
         else:
-            assert u.shape[1]==k.shape[0]
-        k1_vec = np.array([np.linspace(0, k_ / 2, len(u) / 2)[1:] for k_ in k]).T
+            assert u.shape[1] == k.shape[0]
+        k1_vec = np.array([np.linspace(0, k_ / 2, len(u) // 2)[1:] for k_ in k]).T
     else:
-        assert np.abs(np.mean(u))<1
-        if v is not None: assert np.abs(np.mean(v))<1
-        if w is not None: assert np.abs(np.mean(w))<1
-        assert isinstance(k,float)
-        k1_vec = np.linspace(0, k / 2, len(u) / 2)[1:]
+        #assert np.abs(np.mean(u)) < 1
+        if v is not None:
+            assert np.abs(np.mean(v)) < 1
+        if w is not None:
+            assert np.abs(np.mean(w)) < 1
+        assert isinstance(k, float)
+        k1_vec = np.linspace(0, k / 2, int(len(u) / 2))[1:]
     if detrend:
         u, v, w = detrend_wsp(u, v, w)
 
@@ -114,7 +124,7 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
         Sample frequency
     Uvw_lst : array_like
         list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...], v and w are optional
-        
+
     Returns
     -------
     k1, uu[, vv[, ww, uw]] : 2-5 x array_like
@@ -129,34 +139,28 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
     assert isinstance(sample_frq, (int, float))
     Uvw_arr = np.array(Uvw_lst)
     Ncomp = Uvw_arr.shape[1]
-    U,v,w = [Uvw_arr[:,i,:].T for i in range(Ncomp)] + [None]*(3-Ncomp)
+    U, v, w = [Uvw_arr[:, i, :].T for i in range(Ncomp)] + [None] * (3 - Ncomp)
     k = 2 * np.pi * sample_frq / U.mean(0)
-    
-    if v is not None: assert np.abs(np.nanmean(v,0)).max()<1, "Max absolute mean of v is %f"%np.abs(np.nanmean(v,0)).max()
-    if w is not None: assert np.abs(np.nanmean(w,0)).max()<1
-    k1_vec = np.array([np.linspace(0, k_ / 2, U.shape[0] / 2)[1:] for k_ in k]).T
-    u = U-np.nanmean(U, 0)
+
+    if v is not None:
+        assert np.abs(np.nanmean(v, 0)).max() < 1, "Max absolute mean of v is %f" % np.abs(np.nanmean(v, 0)).max()
+    if w is not None:
+        assert np.abs(np.nanmean(w, 0)).max() < 1
+    k1_vec = np.array([np.linspace(0, k_ / 2, U.shape[0] // 2)[1:] for k_ in k]).T
+    u = U - np.nanmean(U, 0)
     u, v, w = detrend_wsp(u, v, w)
 
     return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]]
 
 
-{'ns':( 298, 0.0455033341592, 30.7642421893, 2.32693322102, 0.00538254487239),
-'nu':( 853, 0.0355868134195, 71.0190181051, 2.95576644546, 0.00268903502954),
-'s' :( 367, 0.0300454276184, 24.7375807569, 2.40165270878, 0.00253338624678),
-'vs':( 223, 0.0149281609201, 11.7467239752, 3.21842435078, 0.00143453350246),
-'n' :( 984, 0.0494970534618, 47.3846412548, 2.83466101838, 0.00513954022345),
-'u' :( 696, 0.0258456660845, 78.9511607824, 2.51128584334, 0.00250014140782),
-'vu':( 468, 0.0218274041889, 73.0730383576, 2.01636095317, 0.00196337613117)}
-
 def bin_spectrum(x, y, bin_size, min_bin_count=2):
     assert min_bin_count > 0
     x = x / bin_size
     low, high = np.floor(np.nanmin(x)), np.ceil(np.nanmax(x))
     bins = int(high - low)
     nbr_in_bins = np.histogram(x, bins, range=(low, high))[0]
-    if len(x.shape)==2:
-        min_bin_count*=x.shape[1]
+    if len(x.shape) == 2:
+        min_bin_count *= x.shape[1]
     mask = nbr_in_bins >= min_bin_count
     return np.histogram(x, bins, range=(low, high), weights=y)[0][mask] / nbr_in_bins[mask], nbr_in_bins
 
@@ -171,10 +175,11 @@ def logbin_spectrum(k1, xx, log10_bin_size=.2, min_bin_count=2):
 def logbin_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=0.2, min_bin_count=2):
     return tuple([logbin_spectrum(k1, xx, log10_bin_size, min_bin_count) for xx in [k1, uu, vv, ww, uw]])
 
-def plot_spectrum(spacial_frq, u, plt=None):
+
+def plot_spectrum(spatial_resolution, u, plt=None):
     if plt is None:
         import matplotlib.pyplot as plt
-    k1, uu = logbin_spectra(*spectra(spacial_frq, u))[:2]
+    k1, uu = logbin_spectra(*spectra(spatial_resolution, u))[:2]
 
     plt.semilogx(k1, k1 * uu, 'b-')
 
@@ -187,7 +192,8 @@ def detrend_wsp(u, v=None, w=None):
         t = np.arange(dwsp.shape[0])
         A = np.vstack([t, np.ones(len(t))]).T
         for i in range(dwsp.shape[1]):
-            trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0]
+            m = ~np.isnan(dwsp[:, i])
+            trend, offset = np.linalg.lstsq(A[:m.sum()], dwsp[:, i][m], rcond=None)[0]
             dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
         return dwsp.reshape(wsp.shape)
     return [_detrend(wsp) for wsp in [u, v, w]]
@@ -197,11 +203,12 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2,
     if plt is None:
         import matplotlib.pyplot as plt
     bk1, buu, bvv, bww, buw = logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)
+
     def plot(xx, label, color, plt):
-        plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , marker_style + color, label=label)
+        plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2, marker_style + color, label=label)
     plot(buu, 'uu', 'r', plt)
     plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
-    if mean_u ==1:
+    if mean_u == 1:
         plt.ylabel(r'Spectral density $k_{1} F(k_{1}) [m^2/s^2]$')
     else:
         plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [-]$')
@@ -209,4 +216,4 @@ def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2,
         plot(bvv, 'vv', 'g', plt)
     if bww is not None:
         plot(bww, 'ww', 'b', plt)
-        plot(buw, 'uw', 'm', plt)
\ No newline at end of file
+        plot(buw, 'uw', 'm', plt)

wetb/wind/turbulence/tests/test_man_turbulence.py

+'''
+Created on 20. jul. 2017
+
+@author: mmpe
+'''
+import unittest
+from wetb.utils.test_files import get_test_file
+from wetb.wind.turbulence import mann_turbulence
+from wetb.wind.turbulence.mann_turbulence import parameters2name
+from tests import npt
+
+
+class TestMannTurbulence(unittest.TestCase):
+
+    def testLoad(self):
+        fn = get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat')
+        self.assertRaises(AssertionError, mann_turbulence.load, fn, (31, 31))
+        self.assertRaises(AssertionError, mann_turbulence.load, fn, (8192, 32, 32))
+        u = mann_turbulence.load(fn, (8192, 8, 8))
+        self.assertEqual(u.shape, (8192, 8 * 8))
+
+    def test_loaduvw(self):
+        fn = 'h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat'
+        fn_lst = [get_test_file(fn % uvw) for uvw in 'uvw']
+
+        u, v, w = mann_turbulence.load_uvw(fn_lst, (8192, 8, 8))
+        self.assertEqual(u.shape, (8192, 8 * 8))
+        self.assertTrue(u.shape == v.shape == w.shape == (8192, 8 * 8))
+        u, v, w = mann_turbulence.load_uvw(fn_lst[0].replace("u.dat", "%s.dat"), (8192, 8, 8))
+        self.assertEqual(u.shape, (8192, 8 * 8))
+        self.assertTrue(u.shape == v.shape == w.shape == (8192, 8 * 8))
+
+    def test_parameters2name(self):
+        self.assertEqual(parameters2name((8192, 8, 8), (16384, 32, 32), 0.15, 10, 3.3, 1, 1)[
+                         0], './turb/mann_l10.0_ae0.1500_g3.3_h1_8192x8x8_2.000x4.00x4.00_s0001u.turb')
+
+    def test_save(self):
+        fn = get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat')
+        u_ref = mann_turbulence.load(fn, (8192, 8, 8))
+        mann_turbulence.save(u_ref, 'tmp_u.turb')
+        u = mann_turbulence.load('tmp_u.turb', (8192, 8, 8))
+        npt.assert_array_equal(u, u_ref)
+
+
+#     def test_name2parameters(self):
+#         self.assertEqual(name2parameters('./turb/mann_l10.0_ae0.15_g3.3_h1_8192x8x8_2.000x4.00x4.00_s0001u.turb'),
+#                          ((8192, 8, 8), (16384., 32., 32.), 0.15, 10, 3.3, 1, 1))
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()

wetb/wind/turbulence/tests/test_mann_parameters.py

+'''
+Created on 05/06/2012
+
+@author: Mads
+'''
+import os
+import sys
+import unittest
+
+import matplotlib.pyplot as plt
+import numpy as np
+from wetb.utils.test_files import get_test_file, move2test_files
+from wetb.wind.turbulence.mann_turbulence import load_uvw, fit_mann_parameters,\
+    fit_mann_parameters_from_time_series
+from wetb.utils.timing import print_time
+from wetb.wind.turbulence.spectra import spectra_from_time_series, spectra,\
+    plot_spectra, logbin_spectra
+import warnings
+from wetb.wind.turbulence import mann_turbulence
+from wetb.wind.turbulence.mann_parameters import var2ae, fit_ae, ae2ti
+from numpy import spacing
+
+
+tfp = os.path.join(os.path.dirname(__file__), 'test_files/')
+
+
+class TestMannTurbulence(unittest.TestCase):
+
+    def test_fit_mann_parameters_turbulence_box(self):
+        # for uvw in 'uvw':
+        #    move2test_files(tfp + 'h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat'%uvw)
+        fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
+        u, v, w = load_uvw(fn_lst, (8192, 8, 8))
+        dx = 16384 / 8192
+        fx = 1 / dx  # spatial resolution
+        plt = None
+        ae, L, G = fit_mann_parameters(fx, u, v, w, plt=plt)
+
+        self.assertAlmostEqual(ae, .15, delta=0.01)
+        self.assertAlmostEqual(L, 10, delta=0.3)
+        self.assertAlmostEqual(G, 3.3, delta=0.06)
+
+    def test_spectra_from_timeseries(self):
+        fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
+        u, v, w = load_uvw(fn_lst, (8192, 8, 8))
+        dx = 16384 / 8192
+        fx = 1 / dx  # spatial resolution
+        k1, uu, vv, ww, uw = logbin_spectra(*spectra(fx, u, v, w))
+
+        U = u + 4
+        sample_frq = 2
+        k12, uu2, vv2, ww2, uw2 = logbin_spectra(
+            *spectra_from_time_series(sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
+        np.testing.assert_allclose(uu, uu2, 0.02)
+
+        U = u + 8
+        sample_frq = 2
+        k13, uu3, vv3, ww3, uw3 = logbin_spectra(
+            *spectra_from_time_series(sample_frq, [(U_[::2], v_[::2], w_[::2]) for U_, v_, w_ in zip(U.T, v.T, w.T)]))
+        np.testing.assert_allclose(uu[:-3], uu3[:-2], 0.1)
+
+        # One set of time series with U=4
+        U = u + 4
+        Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
+        # Another set of time series with U=8 i.e. only every second point to have
+        # same sample_frq. (nan added to have same length)
+        U = u + 4
+        Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
+                         np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
+        sample_frq = 2
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore")
+            k14, uu4, vv4, ww4, uw4 = logbin_spectra(*spectra_from_time_series(sample_frq, Uvw_lst))
+        np.testing.assert_allclose(uu[:-3], uu3[:-2], rtol=0.1)
+        if 0:
+            import matplotlib.pyplot as plt
+            plt.semilogx(k1, uu * k1)
+            plt.semilogx(k12, uu2 * k12)
+            plt.semilogx(k13, uu3 * k13)
+            plt.semilogx(k14, uu4 * k14)
+            plt.show()
+
+    def test_fit_mann_parameters_from_timeseries(self):
+        fn_lst = [get_test_file('h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat' % uvw) for uvw in 'uvw']
+        u, v, w = load_uvw(fn_lst, (8192, 8, 8))
+        dx = 16384 / 8192
+        fx = 1 / dx  # spatial resolution
+        ae, L, G = fit_mann_parameters(fx, u, v, w)
+        self.assertAlmostEqual(ae, .15, delta=0.01)
+        self.assertAlmostEqual(L, 10, delta=0.3)
+        self.assertAlmostEqual(G, 3.3, delta=0.06)
+
+        #import matplotlib.pyplot as plt
+        plt = None
+        U = u + 4
+        sample_frq = 2
+        ae, L, G = fit_mann_parameters_from_time_series(
+            sample_frq, [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)], plt=plt)
+        self.assertAlmostEqual(ae, .15, delta=0.01)
+        self.assertAlmostEqual(L, 10, delta=0.3)
+        self.assertAlmostEqual(G, 3.3, delta=0.06)
+
+        # One set of time series with U=4
+        U = u + 4
+        Uvw_lst = [(U_, v_, w_) for U_, v_, w_ in zip(U.T, v.T, w.T)]
+        # Another set of time series with U=8 i.e. only every second point to have
+        # same sample_frq. (nan added to have same length)
+        U = u + 4
+        Uvw_lst.extend([(np.r_[U_[::2], U_[::2] + np.nan], np.r_[v_[::2], v_[::2] + np.nan],
+                         np.r_[w_[::2], w_[::2] + np.nan]) for U_, v_, w_ in zip(U.T, v.T, w.T)])
+        sample_frq = 2
+        ae, L, G = fit_mann_parameters_from_time_series(sample_frq, Uvw_lst, plt=plt)
+        self.assertAlmostEqual(ae, .15, delta=0.01)
+        self.assertAlmostEqual(L, 10, delta=0.3)
+        self.assertAlmostEqual(G, 3.3, delta=0.06)
+
+    def test_var2ae_U(self):
+        u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
+        dx = 2
+        for U in [1, 10, 100]:
+            # should be independent of U
+            dt = dx / U
+            T = 16384 / U
+            self.assertAlmostEqual(var2ae(variance=u.var(), L=10, G=3.3, U=U, T=T, sample_frq=1 / dt), .15, delta=.024)
+
+    def test_var2ae_T(self):
+        u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
+        dx = 2
+        U = 10
+        dt = dx / U
+        for i in np.arange(6):
+            # reshape to more and shorter series. Variance should decrease while ae should be roughly constant
+            n = 2**i
+            u_ = u.T.reshape((u.T.shape * np.array([n, 1 / n])).astype(int)).T
+            var = u_.var(0).mean()
+            ae = var2ae(variance=var, L=10, G=3.3, U=U, T=dx * u_.shape[0] / U, sample_frq=1 / dt)
+            self.assertAlmostEqual(ae, .15, delta=.027)
+
+    def test_var2ae_dt(self):
+        u = mann_turbulence.load(get_test_file("h2a16384_8_8_65536_32_32_0.15_40_4.0u.dat"), (16384, 8, 8))
+        dx = 4
+        U = 10
+        T = u.shape[0] * dx / U
+        for i in np.arange(9):
+            # average every neighbouring samples to decrease dt.
+            # Variance should decrease while ae should be roughly constant
+            n = 2**i
+
+            u_ = u.reshape(u.shape[0] // n, n, u.shape[1]).mean(1)
+            var = u_.var(0).mean()
+
+            ae = var2ae(variance=var, L=40, G=4, U=U, T=T, sample_frq=1 / (n * dx / U))
+            self.assertAlmostEqual(ae, .15, delta=.041 + i / 100)
+
+    def test_fit_ae2var(self):
+        u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
+        self.assertAlmostEqual(fit_ae(spatial_resolution=2, u=u, L=10, G=3.3), .15, delta=.02)
+
+    def test_ae2ti(self):
+        u = mann_turbulence.load(get_test_file("h2a8192_8_8_16384_32_32_0.15_10_3.3u.dat"), (8192, 8, 8))
+        dx = 2
+        U = 10
+        dt = dx / U
+        T = 16384 / U
+        ae23 = var2ae(variance=u.var(), L=10, G=3.3, U=U, T=T, sample_frq=1 / dt)
+        ti = u.std() / U
+        self.assertAlmostEqual(ae2ti(ae23, L=10, G=3.3, U=U, T=T, sample_frq=1 / dt), ti)
+
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()

wetb/wind/utils.py

@@ -38,14 +38,15 @@ def wsp_dir2uv(wsp, dir, dir_ref=None):
     v = -np.sin(rad(dir)) * wsp[:]
     return np.array([u, v])
 
+
 def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
-    """Convert horizontal wind speed and direction to u,v,w
+    r"""Convert horizontal wind speed and direction to u,v,w
 
     Parameters
     ----------
     wsp : array_like
-        - if wsp_horizontal is True: Horizontal wind speed, $\sqrt{u^2+v^2}\n
-        - if wsp_horizontal is False: Wind speed, $\sqrt{u^2+v^2+w^2}
+        - if wsp_horizontal is True: Horizontal wind speed, $\sqrt{u^2+v^2}$\n
+        - if wsp_horizontal is False: Wind speed, $\sqrt{u^2+v^2+w^2}$
     dir : array_like
         Wind direction
     tilt : array_like
@@ -66,6 +67,7 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
     w : array_like
         v wind component
     """
+
     wsp, dir, tilt = wsp[:], dir[:], tilt[:]
     if wsp_horizontal:
         w = tand(tilt) * wsp
@@ -76,7 +78,6 @@ def wsp_dir_tilt2uvw(wsp, dir, tilt, wsp_horizontal, dir_ref=None):
     return np.array([u, v, w])
 
 
-
 def xyz2uvw(x, y, z, left_handed=True):
     """Convert sonic x,y,z measurements to u,v,w wind components
 
@@ -108,7 +109,7 @@ def xyz2uvw(x, y, z, left_handed=True):
     SV = cosd(theta) * y - sind(theta) * x
 
 #     SUW = cosd(theta) * x + sind(theta) * y
-# 
+#
 #     #% rotation around y of tilt
 #     tilt = deg(np.arctan2(np.mean(z), np.mean(SUW)))
 #     SU = SUW * cosd(tilt) + z * sind(tilt);
@@ -120,8 +121,6 @@ def xyz2uvw(x, y, z, left_handed=True):
     return np.array([SU, SV, SW])
 
 
-
-
 def abvrel2xyz_old(alpha, beta, vrel):
     """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
 
@@ -143,27 +142,29 @@ def abvrel2xyz_old(alpha, beta, vrel):
     z : array_like
         Wind component in beta plane (positive for negative beta)
     """
-    alpha = np.array(alpha, dtype=np.float)
-    beta = np.array(beta, dtype=np.float)
-    vrel = np.array(vrel, dtype=np.float)
+    alpha = np.array(alpha, dtype=float)
+    beta = np.array(beta, dtype=float)
+    vrel = np.array(vrel, dtype=float)
 
     sign_vsx = -((np.abs(beta) > np.pi / 2) * 2 - 1)  # +1 for |beta| < 90, -1 for |beta|>90
-    sign_vsy = np.sign(alpha)  #+ for alpha > 0
-    sign_vsz = -np.sign(beta)  #- for beta>0
-
+    sign_vsy = np.sign(alpha)  # + for alpha > 0
+    sign_vsz = -np.sign(beta)  # - for beta>0
 
     x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
 
     m = alpha != 0
     y = np.zeros_like(alpha)
-    y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 + 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
+    y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 +
+                                                 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
 
     m = beta != 0
     z = np.zeros_like(alpha)
-    z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
+    z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 +
+                                                 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
 
     return x, y, z
 
+
 def abvrel2xyz(alpha, beta, vrel):
     """Convert pitot tube alpha, beta and relative velocity to local Cartesian wind speed velocities
 
@@ -197,41 +198,42 @@ def abvrel2xyz(alpha, beta, vrel):
     z : array_like
         Wind component in beta plane (positive for negative beta)
     """
-    alpha = np.array(alpha, dtype=np.float)
-    beta = np.array(beta, dtype=np.float)
-    vrel = np.array(vrel, dtype=np.float)
+    alpha = np.array(alpha, dtype=float)
+    beta = np.array(beta, dtype=float)
+    vrel = np.array(vrel, dtype=float)
 
     sign_vsx = ((np.abs(beta) > np.pi / 2) * 2 - 1)  # -1 for |beta| < 90, +1 for |beta|>90
-    sign_vsy = -np.sign(alpha)  #- for alpha > 0
+    sign_vsy = -np.sign(alpha)  # - for alpha > 0
     sign_vsz = np.sign(beta)  # for beta>0
 
-
     x = sign_vsx * np.sqrt(vrel ** 2 / (1 + np.tan(alpha) ** 2 + np.tan(beta) ** 2))
 
     m = alpha != 0
     y = np.zeros_like(alpha)
-    y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 + 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
+    y[m] = sign_vsy[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(alpha[m])) ** 2 +
+                                                 1 + (np.tan(beta[m]) / np.tan(alpha[m])) ** 2))
 
     m = beta != 0
     z = np.zeros_like(alpha)
-    z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 + 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
+    z[m] = sign_vsz[m] * np.sqrt(vrel[m] ** 2 / ((1 / np.tan(beta[m])) ** 2 +
+                                                 1 + (np.tan(alpha[m]) / np.tan(beta[m])) ** 2))
 
     return np.array([x, y, z]).T
 
 
 def detrend_uvw(u, v=None, w=None):
-#     def _detrend(wsp):
-#         if wsp is None:
-#             return None
-#         dwsp = np.atleast_2d(wsp.copy().T).T
-#         t = np.arange(dwsp.shape[0])
-#         A = np.vstack([t, np.ones(len(t))]).T
-#         for i in range(dwsp.shape[1]):
-#             trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0]
-#             dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
-#         return dwsp.reshape(wsp.shape)
+    #     def _detrend(wsp):
+    #         if wsp is None:
+    #             return None
+    #         dwsp = np.atleast_2d(wsp.copy().T).T
+    #         t = np.arange(dwsp.shape[0])
+    #         A = np.vstack([t, np.ones(len(t))]).T
+    #         for i in range(dwsp.shape[1]):
+    #             trend, offset = np.linalg.lstsq(A, dwsp[:, i])[0]
+    #             dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
+    #         return dwsp.reshape(wsp.shape)
     def _detrend(y):
         if y is None:
             return None
         return detrend(y)
-    return [_detrend(uvw) for uvw in [u, v, w]] 
\ No newline at end of file
+    return [_detrend(uvw) for uvw in [u, v, w]]

wetb/wind/weibull.py

@@ -103,7 +103,7 @@ def fit(wsp):
 if __name__ == "__main__":
     from wetb.wind import weibull
     from pylab import hist, show, plot
-    hist(weibull.random(10, 2, 10000), 20, normed=True)
+    hist(weibull.random(10, 2, 10000), 20, density=True)
     wsp = np.arange(0, 20, .5)
     plot(wsp, weibull.pdf(10, 2)(wsp))
     show()