added missing test file

wetb/gtsdf/__init__.py

@@ -44,6 +44,7 @@ from .gtsdf import compress2statistics
 
 class Dataset(object):
     def __init__(self, filename):
+        self.filename = filename
         self.time, self.data, self.info = load(filename)
     def __call__(self, id):
         if isinstance(id, str):

wetb/hawc2/at_time_file.py

@@ -49,8 +49,9 @@ class AtTimeFile(object):
         self.blade_radius = bladetip_radius
         with open(filename, encoding='utf-8') as fid:
             lines = fid.readlines()
-        atttribute_name_line = [l.startswith("# Radius_s") for l in lines].index(True)
-        self.attribute_names = lines[atttribute_name_line].lower().replace("#", "").split()
+        atttribute_name_line = [l.strip().startswith("# Radius_s") for l in lines].index(True)
+        #self.attribute_names = lines[atttribute_name_line].lower().replace("#", "").split()
+        self.attribute_names = [n.strip() for n in lines[atttribute_name_line].lower().split("#")[1:]]
         data = np.array([[float(l) for l in lines[i].split() ] for i in range(atttribute_name_line+1, len(lines))])
         self.data = data
         def func_factory(column):
@@ -97,10 +98,15 @@ class AtTimeFile(object):
 
 if __name__ == "__main__":
     at = AtTimeFile(r"tests/test_files/at.dat", 86.3655) # load file
+    at = AtTimeFile(r'U:\hama\HAWC2-AVATAR\res/avatar-7ntm-scaled-rad.dat')
+    
     at.attribute_names # Attribute names
     at[:3,1] # first 3 twist rows
-    at.twist()[:3] # Twist first 3 twist rows
-    print (at.twist(36, curved_length=True)) # Twist at curved_length = 36 (interpolated)
-    print (at.twist(36)) # Twist at 36 (interpolated)
+    print (len(at.attribute_names))
+    print ("\n".join(at.attribute_names))
+    print (at.data.shape)
+    #at.twist()[:3] # Twist first 3 twist rows
+    #print (at.twist(36, curved_length=True)) # Twist at curved_length = 36 (interpolated)
+    #print (at.twist(36)) # Twist at 36 (interpolated)
     
 

wetb/signal/filters/cy_filters.py

@@ -42,6 +42,33 @@ def cy_dynamic_low_pass_filter(inp, delta_t, tau, method=1):  #cpdef cy_dynamic_
             output[i] = output[i - 1] * np.exp(-delta_t / tau[i]) + inp[i] * (1 - np.exp(-delta_t / tau[i]))
     return output
 
+def cy_dynamic_low_pass_filter_2d(inp, delta_t, tau, method=1):  #cpdef cy_dynamic_low_pass_filter_2d(np.ndarray[double,ndim=2] inp, double delta_t, np.ndarray[double,ndim=1] tau, int method=1):
+    #cdef np.ndarray[double,ndim=2] output, alpha
+    #cdef int i
+
+    output = np.empty_like(inp, dtype=np.float)
+    output[0] = inp[0]
+
+    if method == 1:
+        alpha = delta_t / (tau + delta_t)
+        for i in range(1, inp.shape[0]):
+            output[i] = output[i - 1] + alpha[i] * (inp[i] - output[i - 1])  # Same as output[i] = alpha*inp[i]+(1-alpha)*output[i-1]
+    elif method == 2:
+        for i in range(1, inp.shape[0]):
+            output[i] = (delta_t * (inp[i] + inp[i - 1] - output[i - 1]) + 2 * tau[i] * output[i - 1]) / (delta_t + 2 * tau[i])
+    elif method == 3:
+        for i in range(1, inp.shape[0]):
+            output[i] = output[i - 1] * np.exp(-delta_t / tau[i]) + inp[i] * (1 - np.exp(-delta_t / tau[i]))
+    return output
+
+def cy_dynamic_low_pass_filter_test(inp):  #cpdef cy_dynamic_low_pass_filter_test(np.ndarray[double,ndim=2] inp):
+    #cdef np.ndarray[double,ndim=2] output, alpha
+    #cdef int i
+    output = np.empty_like(inp, dtype=np.float)
+    output[0] = inp[0]
+    for i in range(1, inp.shape[0]):
+        output[i] = inp[i]
+    return output
 
 @cython.ccall
 @cython.locals(alpha=cython.float, i=cython.int)

wetb/signal/filters/first_order.py

@@ -11,7 +11,10 @@ def low_pass(input, delta_t, tau, method=1):
     if isinstance(tau, (int, float)):
         return cy_filters.cy_low_pass_filter(input.astype(np.float64), delta_t, tau)
     else:
-        return cy_filters.cy_dynamic_low_pass_filter(input.astype(np.float64), delta_t, tau, method)
+        if len(input.shape)==2:
+            return cy_filters.cy_dynamic_low_pass_filter_2d(input.astype(np.float64), delta_t, tau, method)
+        else:
+            return cy_filters.cy_dynamic_low_pass_filter(input.astype(np.float64), delta_t, tau, method)
 
 def high_pass(input, delta_t, tau):
     from wetb.signal.filters import cy_filters

wetb/signal/subset_mean.py

@@ -137,7 +137,7 @@ def cycle_trigger(values, trigger_value=None, step=1, ascending=True, tolerance=
     else:
         return np.where((values[1:] < trigger_value - tolerance) & (values[:-1] >= trigger_value + tolerance))[0][::step]
 
-def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1):
+def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
     """Returns one index per revolution (minimum rotor position)
     
     Parameters
@@ -157,19 +157,25 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1):
         rotor_speed = np.ones_like(rotor_position)*rotor_speed
     deg_per_sample = rotor_speed*360/60/sample_frq
     thresshold = deg_per_sample.max()*3
-    drp = np.diff(rotor_position)%360
+    drp = (np.diff(rotor_position)+thresshold)%360-thresshold
     rp = rotor_position
 
     
 
     
-    rp = np.array(rotor_position).copy()
+    rp = np.array(rotor_position).copy()%360
     #filter degree increase > thresshold
     rp[np.r_[False, np.diff(rp)>thresshold]] = 180
     
     upper_indexes = np.where((rp[:-1]>(360-thresshold))&(rp[1:]<(360-thresshold)))[0]
     lower_indexes = np.where((rp[:-1]>thresshold)&(rp[1:]<thresshold))[0] +1 
     
+    if plt:
+        plt.plot(rp)
+        plt.plot(lower_indexes, rp[lower_indexes],'.')
+        plt.plot(upper_indexes, rp[upper_indexes],'.')
+        
+
     # Best lower is the first lower after upper
     best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
     upper2lower = best_lower - upper_indexes

wetb/wind/turbulence/mann_parameters.py

@@ -189,13 +189,13 @@ def residual(ae, L, G, k1, uu, vv=None, ww=None, uw=None, log10_bin_size=.2):
 
     Returns
     -------
-    residual : float
+    residual : array_like
+        rms of each spectrum
     """
-    _3to2list = list(np.array(logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size)))
-    bk1, sp_meas = _3to2list[:1] + [_3to2list[1:]]
-    sp_fit = np.array(logbin_spectra(k1, *get_mann_model_spectra(ae, L, G, k1)))[1:]
-
-    return np.sqrt(((bk1 * sp_meas - bk1 * sp_fit) ** 2).mean())
+    k1_sp = np.array([sp for sp in logbin_spectra(k1, uu, vv, ww, uw, log10_bin_size) if sp is not None])
+    bk1, sp_meas = k1_sp[0], k1_sp[1:]
+    sp_fit = np.array(get_mann_model_spectra(ae, L, G, bk1))[:sp_meas.shape[0]]
+    return np.sqrt(((bk1 * (sp_meas - sp_fit)) ** 2).mean(1))
 
 
 def fit_ae2var(variance, L, G):

wetb/wind/turbulence/spectra.py

@@ -113,7 +113,7 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
     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)...] 
+        list of U, v and w, [(U1,v1,w1),(U2,v2,w2)...], v and w are optional
         
     Returns
     -------
@@ -127,7 +127,9 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
         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)]
+    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)
     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()
@@ -138,6 +140,15 @@ def spectra_from_time_series(sample_frq, Uvw_lst):
 
     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