first test_file implementation

wetb/utils/test_files.py

+'''
+Created on 20. jul. 2017
+
+@author: mmpe
+'''
+import os
+import wetb
+import urllib.request
+from fileinput import filename
+import inspect
+wetb_rep_path = os.path.join(os.path.dirname(wetb.__file__), "../")                                   
+default_TestFile_rep_path=os.path.join(os.path.dirname(wetb.__file__) + "/../../TestFiles/")
+
+def get_test_file(filename):
+    if not os.path.isabs(filename):
+        index = [os.path.realpath(s[1]) for s in inspect.stack()].index(__file__) + 1
+        tfp = os.path.dirname(inspect.stack()[index][1]) + "/test_files/"
+        filename = tfp + filename
+    
+    if os.path.exists(filename):
+        return filename
+    else:
+        filename2 = os.path.realpath(os.path.join(wetb_rep_path, 'downloaded_test_files', os.path.relpath(filename, wetb_rep_path)))
+        if not os.path.isfile(filename2):
+            #url = 'https://gitlab.windenergy.dtu.dk/toolbox/TestFiles/%s'%os.path.relpath(filename, wetb_rep_path)
+            url = 'http://tools.windenergy.dtu.dk/TestFiles/%s.txt'%os.path.relpath(filename, wetb_rep_path).replace("\\","/")
+            print ("download %s\nfrom %s"%(filename, url))
+            if not os.path.exists(os.path.dirname(filename2)):
+                os.makedirs(os.path.dirname(filename2))
+            urllib.request.urlretrieve(url, filename2)
+        return filename2
+
+
+
+def move2test_files(filename,TestFile_rep_path=default_TestFile_rep_path):
+    wetb_rep_path = os.path.join(os.path.dirname(wetb.__file__), "../")
+    folder = os.path.dirname(TestFile_rep_path + os.path.relpath(filename, wetb_rep_path))
+    if not os.path.exists(folder):
+        os.makedirs(folder)
+    os.rename(filename, os.path.join(folder, os.path.basename(filename)+'.txt'))
+    
+    
+    
+    
+    
\ No newline at end of file

wetb/utils/tests/test_test_files.py

+'''
+Created on 20. jul. 2017
+
+@author: mmpe
+'''
+import unittest
+from wetb.utils.test_files import move2test_files, get_test_file
+import os
+from wetb.utils import test_files
+import wetb
+
+
+tfp = os.path.join(os.path.dirname(__file__) + "/test_files/")
+class Test_test_files(unittest.TestCase):
+    def test_move2test_files(self):
+        dst = test_files.default_TestFile_rep_path+ "wetb/utils/tests/test_files/tmp_test_file.txt"
+        src = tfp+'tmp_test_file.txt'
+        if os.path.isdir(test_files.default_TestFile_rep_path):
+            if os.path.isfile(dst):
+                os.remove(dst)
+            if not os.path.isfile(src):
+                with open(src,'w') as fid:
+                    fid.write("This is a test file")
+            move2test_files(src)
+            self.assertTrue(os.path.isfile(dst))
+
+    def test_test_files(self):
+        fn = os.path.realpath(os.path.dirname(wetb.__file__) + '/../downloaded_test_files/wetb/utils/tests/test_files/test_file.txt')
+        if os.path.isfile(fn):
+            os.remove(fn)
+        fn1 = get_test_file(tfp+'test_file.txt')
+        self.assertTrue(fn1)
+        os.remove(fn1)
+        fn2 = get_test_file('test_file.txt')
+        self.assertEqual(fn2, fn)
+        self.assertTrue(os.path.isfile(fn2))
+        
+
+if __name__ == "__main__":
+    #import sys;sys.argv = ['', 'Test.testName']
+    unittest.main()
\ No newline at end of file

wetb/wind/turbulence/mann_parameters.py

@@ -11,7 +11,8 @@ from scipy.interpolate import RectBivariateSpline
 
 import numpy as np
 from wetb.wind.turbulence.spectra import spectra, logbin_spectra, \
-    plot_spectrum
+    plot_spectrum, plot_spectra
+from wetb.utils.test_files import move2test_files
 
 
 
@@ -24,12 +25,7 @@ RBS2 = RectBivariateSpline(xp, yp, sp2)
 RBS3 = RectBivariateSpline(xp, yp, sp3)
 RBS4 = RectBivariateSpline(xp, yp, sp4)
 
-def estimate_mann_parameters(sf, u, v, w=None):
-    if isinstance(u, (list, tuple)):
-        #return fit_mann_model_spectra(*mean_spectra(sf, u, v, w))
-        raise NotImplementedError
-    else:
-        return fit_mann_model_spectra(*spectra(sf, u, v, w))
+
 
 #def mean_spectra(fs, u_ref_lst, u_lst, v_lst=None, w_lst=None):
 #    if isinstance(fs, (int, float)):
@@ -145,11 +141,13 @@ def fit_mann_model_spectra(k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2,
     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)
 
-    if plt is not False:
+    if plt:
         if not hasattr(plt, 'plot'):
             import matplotlib.pyplot as plt
-        _plot_spectra(k1, uu, vv, ww, uw, plt=plt)
-        plot_mann_spectra(*x, plt=plt)
+#         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)
         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]$')
@@ -254,26 +252,13 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
     return ae
 
 
-def plot_spectra(ae, L, G, k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+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
-    _plot_spectra(k1, uu, vv, ww, uw, mean_u, log10_bin_size, 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_spectra(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
-    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 , '.' + color, label=label)
-    plot(buu, 'uu', 'r', plt)
-    plt.xlabel('Wavenumber $k_{1}$ [$m^{-1}$]')
-    plt.ylabel(r'Spectral density $k_{1} F(k_{1})/U^{2} [m^2/s^2]$')
-    if (bvv) is not None:
-        plot(bvv, 'vv', 'g', plt)
-    if bww is not None:
-        plot(bww, 'ww', 'b', plt)
-        plot(buw, 'uw', 'm', plt)
+
 
 def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv', 'ww', 'uw']):
     if plt is None:
@@ -312,7 +297,8 @@ if __name__ == "__main__":
     nx = 8192
     ny, nz = 8, 8
     sf = (nx / l)
-    fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
+    #fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
+    #fn = os.path.dirname(wind.__file__)+"/tests/test_files/turb/h2a8192_8_8_16384_32_32_0.15_10_3.3%s.dat"
     u, v, w = [np.fromfile(fn % uvw, np.dtype('<f'), -1).reshape(nx , ny * nz) for uvw in ['u', 'v', 'w']]
     ae, L, G = fit_mann_model_spectra(*spectra(sf, u, v, w), plt=plt)
     print (ae, L, G)

wetb/wind/turbulence/spectra.py

@@ -35,20 +35,21 @@ def spectrum(x, y=None, k=1):
 
     fftx = fftx[:len(fftx) // 2 ] * 2  # positive half * 2
 
-    if len(fftx.shape) == 2:
-        fftx = np.mean(fftx, 1)
+#    if len(fftx.shape) == 2:
+#        fftx = np.mean(fftx, 1)
 
     return np.real(fftx * len(x) / (2 * k))[1:]
 
-def spectra(spacial_frq, u, v=None, w=None, detrend=True):
+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
     ----------
-    spacial_frq : inf or float
-        - For time series: Sample frequency, see Notes\n
-        - For turbulence boxes: 1/dx where
-        dx is nx/Lx (number of grid points in length direction / box length in meters)
+    spatial_resolution : int, float or array_like
+        Distance between samples in meters
+        - 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
         The u-wind component\n
         - if shape is (r,): One time series with *r* observations\n
@@ -74,36 +75,75 @@ def spectra(spacial_frq, u, v=None, w=None, detrend=True):
         - uw: The uw cross spectrum, only if w is provided\n
 
         For two dimensional input, the mean spectra are returned
-
-    Notes
-    -----
-    If the 'mean wind speed' or 'shear compensated reference wind speed' are subtracted from u in advance,
-    then the spacial_frq must be divided by the subtracted wind speed
     """
-    assert isinstance(spacial_frq, (int, float))
-
-    k = 2 * np.pi * spacial_frq
-
-    if np.mean(u) > 1:
-        k /= np.mean(u)
-        u = u - np.mean(u, 0)
-    if v is not None:
-        v = v - np.mean(v, 0)
-    if w is not None:
-        w = w - np.mean(w, 0)
+    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 isinstance(k, float):
+            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
+    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:]
     if detrend:
         u, v, w = detrend_wsp(u, v, w)
 
-    k1_vec = np.linspace(0, k / 2, len(u) / 2)[1:]
     return [k1_vec] + [spectrum(x1, x2, k=k) for x1, x2 in [(u, u), (v, v), (w, w), (w, u)]]
 
 
+def spectra_from_time_series(sample_frq, Uvw_lst):
+    """Return the wave number, the uu, vv, ww autospectra and the uw cross spectra
+
+    Parameters
+    ----------
+    sample_frq : int, float or array_like
+        Sample frequency
+    Uvw_lst : array_like
+        list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...] 
+        
+    Returns
+    -------
+    k1, uu[, vv[, ww, uw]] : 2-5 x array_like
+        - k1: wave number
+        - uu: The uu auto spectrum\n
+        - vv: The vv auto spectrum, only if v is provided\n
+        - ww: The ww auto spectrum, only if w is provided\n
+        - uw: The uw cross spectrum, only if w is provided\n
+
+        For two dimensional input, the mean spectra are returned
+    """
+    assert isinstance(sample_frq, (int, float))
+    U,v,w = [np.array(Uvw_lst)[:,i,:].T for i in range(3)]
+    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)
+    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)]]
+
 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.min(x)), np.ceil(np.max(x))
+    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]
     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
 
@@ -138,3 +178,22 @@ def detrend_wsp(u, v=None, w=None):
             dwsp[:, i] = dwsp[:, i] - t * trend + t[-1] / 2 * trend
         return dwsp.reshape(wsp.shape)
     return [_detrend(wsp) for wsp in [u, v, w]]
+
+
+def plot_spectra(k1, uu, vv=None, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None, marker_style='.'):
+    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)
+    plot(buu, 'uu', 'r', plt)
+    plt.xlabel('Wavenumber $k_{1}$ [$m^{-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} [-]$')
+    if (bvv) is not None:
+        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