bugfix + improved revolution trigger

wetb/hawc2/htc_contents.py

@@ -120,7 +120,7 @@ class HTCContents(object):
         self._add_contents(section)
         return section
 
-    def add_line(self, name, values, comments):
+    def add_line(self, name, values, comments=""):
         line = HTCLine(name, values, comments)
         self._add_contents(line)
         return line

wetb/signal/fix/_rotor_position.py

@@ -5,6 +5,8 @@ Created on 30. mar. 2017
 '''
 import numpy as np
 from wetb.signal.fit._spline_fit import spline_fit
+from wetb.signal.filters._differentiation import differentiation
+
 def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt=None):
     """Rotor position fitted with spline
     
@@ -33,7 +35,7 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
     from wetb.signal.subset_mean import revolution_trigger
     
     t = np.arange(len(rotor_position))
-    indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_no_round_diff=4)
+    indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=4)
     
     rp = rotor_position[:].copy()
     
@@ -69,7 +71,9 @@ def fix_rotor_position(rotor_position, sample_frq, rotor_speed, fix_dt=None, plt
         plt.plot(t/sample_frq, spline_fit(x,y)(t)-t*a, label='Spline (detrended)')
         plt.legend()
         plt.show()
-    return spline_fit(x,y)(t)%360
+    fit = spline_fit(x,y)(t)
+    fit[differentiation(fit, "left")<0]=np.nan
+    return fit%360
 
 def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
     """Find the optimal fix_dt parameter for fix_rotor_position (function above).
@@ -95,7 +99,8 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
     """
     from wetb.signal.filters import differentiation
     def err(i):
-        rpm_pos = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))%180 / 360 * sample_frq * 60
+        drp = differentiation(fix_rotor_position(rotor_position, sample_frq, rotor_speed, i))
+        rpm_pos = drp%180 / 360 * sample_frq * 60
         return np.sum((rpm_pos - rotor_speed)**2)
     
     best = 27
@@ -110,4 +115,69 @@ def find_fix_dt(rotor_position, sample_frq, rotor_speed, plt=None):
         plt.show()
     
     return best 
-    
\ No newline at end of file
+
+def check_rotor_position(rotor_position, sample_frq, rotor_speed, max_rev_diff=1, plt=None):
+    """Rotor position fitted with spline
+    
+    Parameters
+    ----------
+    rotor_position : array_like
+        Rotor position [deg] (0-360)
+    sample_frq : int or float
+        Sample frequency [Hz]
+    rotor_speed : array_like
+        Rotor speed [RPM]
+    fix_dt : int, float or None, optional
+        Time distance [s] between spline fix points\n
+        If None (default) a range of seconds is tested and the result that minimize the RMS 
+        between differentiated rotor position fit and rotor speed is used.\n 
+        Note that a significant speed up is achievable by specifying the parameter
+    plt : PyPlot or None
+        If PyPlot a visual interpretation is plotted
+            
+    Returns
+    -------
+    y : nd_array
+        Fitted rotor position
+    """
+    
+    from wetb.signal.subset_mean import revolution_trigger
+    
+    t = np.arange(len(rotor_position))/sample_frq
+    indexes = revolution_trigger(rotor_position[:], sample_frq, rotor_speed, max_rev_diff=max_rev_diff)
+    
+    rotor_position = rotor_position[:]%360
+        
+    rp_fit = fix_rotor_position(rotor_position, sample_frq, rotor_speed, None)
+        
+    if plt:
+        #a = (rp.max()-rp.min())/t.max()
+        #plt.plot(t/sample_frq,rp-t*a,label='Continus rotor position (detrended)')
+        f, axarr = plt.subplots(2)
+        print (rp_fit)
+        axarr[0].plot(t, rotor_position)
+        axarr[0].plot(indexes/sample_frq, rotor_position[indexes],'.')
+        axarr[0].plot(t, rp_fit)
+        axarr[1].plot(t, rotor_speed)
+        axarr[1].plot(indexes/sample_frq,60/( differentiation(indexes)/sample_frq))
+        print (t[170:200])
+        drp = differentiation(rp_fit)
+        #drp[(drp<0)]= 0
+        axarr[1].plot(t, drp%180 / 360 * sample_frq * 60, label='fix')
+        plt.legend()
+    
+    i1,i2 = indexes[0], indexes[-1]
+    print ("Rev from rotor position", np.sum(np.diff(rotor_position[i1:i2])%180)/360)
+    print ("Rev from rotor speed", np.mean(rotor_speed[i1:i2])*(i2-i1)/60/sample_frq)
+    print ("Rev from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360)
+    
+    print ("Mean RPM from rotor speed", np.mean(rotor_speed))
+    print ("Mean RPM from fitted rotor position", np.sum(np.diff(rp_fit[i1:i2])%180)/360 / ((i2-i1)/60/sample_frq))
+    
+    spr1 = (np.diff(indexes)/sample_frq)
+    
+    #rs1 = 60/( np.diff(indexes)/sample_frq)
+    spr2 = np.array([60/rotor_speed[i1:i2].mean() for i1,i2 in zip(indexes[:-1], indexes[1:])])
+    
+    err = spr1-spr2
+    print (err.max())
\ No newline at end of file

wetb/signal/subset_mean.py

@@ -9,6 +9,7 @@ import unittest
 from wetb.utils.geometry import rpm2rads
 from _collections import deque
 from tables.tests.test_index_backcompat import Indexes2_0TestCase
+from wetb.signal.filters._differentiation import differentiation
 
 
 def power_mean(power, trigger_indexes, I, rotor_speed, time, air_density=1.225, rotor_speed_mean_samples=1) :
@@ -136,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_no_round_diff=1):
+def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_rev_diff=1):
     """Returns one index per revolution (minimum rotor position)
     
     Parameters
@@ -155,16 +156,12 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
     if isinstance(rotor_speed, (float, int)):
         rotor_speed = np.ones_like(rotor_position)*rotor_speed
     deg_per_sample = rotor_speed*360/60/sample_frq
-    sample_per_round = 1/(rotor_speed/60/sample_frq)
-    thresshold = deg_per_sample.max()*2
-    
-    nround_rotor_speed = np.nansum(rotor_speed/60/sample_frq)
-    
-    mod = [v for v in [5,10,30,60,90] if v>thresshold][0]
+    thresshold = deg_per_sample.max()*3
+    drp = np.diff(rotor_position)%360
+    rp = rotor_position
+
     
-    nround_rotor_position = np.nansum(np.diff(rotor_position)%mod)/360
-    #assert abs(nround_rotor_position-nround_rotor_speed)<max_no_round_diff, "No of rounds from rotor_position (%.2f) mismatch with no_rounds from rotor_speed (%.2f)"%(nround_rotor_position, nround_rotor_speed)
-    #print (nround_rotor_position, nround_rotor_speed)
+
     
     rp = np.array(rotor_position).copy()
     #filter degree increase > thresshold
@@ -176,10 +173,47 @@ def revolution_trigger(rotor_position, sample_frq, rotor_speed, max_no_round_dif
     # Best lower is the first lower after upper
     best_lower = lower_indexes[np.searchsorted(lower_indexes, upper_indexes)]
     upper2lower = best_lower - upper_indexes
-    best_lower = best_lower[upper2lower<upper2lower.mean()*2]
-    #max_dist_error = max([np.abs((i2-i1)- np.mean(sample_per_round[i1:i2])) for i1,i2 in zip(best_lower[:-1], best_lower[1:])])
-    #assert max_dist_error < sample_frq/5, max_dist_error
-    return best_lower
+    trigger_indexes = best_lower[upper2lower<upper2lower.mean()*2].tolist()
+
+    if len(trigger_indexes)>1:
+        
+    
+        
+        rpm_rs = np.array([rev.mean() for rev in np.split(rotor_speed, trigger_indexes)[1:-1]])
+        rpm_i = 1/np.diff(trigger_indexes)*60*sample_frq
+        spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+        spr_i = np.diff(trigger_indexes)
+        
+
+        while np.any(spr_rs>spr_i*1.9):
+            i = np.where(spr_rs>spr_i*1.9)[0][0]
+            if np.abs(spr_i[i-1] + spr_i[i] - spr_rs[i])< np.abs(spr_i[i] + spr_i[i+1] - spr_rs[i]):
+                del trigger_indexes[i]
+            else:
+                del trigger_indexes[i+1]
+            spr_i = np.diff(trigger_indexes)
+            spr_rs = np.array([rev.mean() for rev in np.split(1/rotor_speed*60*sample_frq, trigger_indexes)[1:-1]])
+        
+            
+        
+        # if a revolution is too long add trigger
+        if np.any(rpm_rs>rpm_i*2.1):
+            # >1 missing triggers
+            raise NotImplementedError
+        trigger_indexes.extend([np.mean(trigger_indexes[i:i+2]) for i in np.where(rpm_rs>rpm_i*1.9)[0]])
+        trigger_indexes = np.sort(trigger_indexes).astype(np.int)
+    
+    
+        i1,i2 = trigger_indexes[0], trigger_indexes[-1]
+        nround_rotor_speed = np.nansum(rotor_speed[i1:i2]/60/sample_frq)
+        #mod = [v for v in [5,10,30,60,90] if v>thresshold][0]    
+        nround_rotor_position = len(trigger_indexes)-1 #np.nansum(np.diff(rotor_position[i1:i2])%mod)/360
+        if max_rev_diff is not None:
+            diff_pct = abs(nround_rotor_position-nround_rotor_speed)/nround_rotor_position*100
+            assert diff_pct<max_rev_diff, "No of revolution mismatch: rotor_position (%d), rotor_speed (%.1f), diff %.1f%%"%(nround_rotor_position, nround_rotor_speed, diff_pct)
+         
+    
+    return trigger_indexes
     
 
 def revolution_trigger_old(values, rpm_dt=None, dmin=5, dmax=10, ):

wetb/signal/tests/test_fix.py

@@ -42,6 +42,8 @@ class TestFix(unittest.TestCase):
          
         self.assertEqual(find_fix_dt(ds.azi, sample_frq, ds.Rot_cor), 4)
         
+        
+        
 
 if __name__ == "__main__":
     #import sys;sys.argv = ['', 'Test.testRotorPositionFix']

wetb/signal/tests/test_subset_mean.py

@@ -92,21 +92,32 @@ class TestSubsetMean(unittest.TestCase):
             x2 = np.random.randint(0, len(rotor_position),10)
             print (list(x2))
         else:
-            x1 = [447, 854, 595, 804, 847, 488, 412, 199, 675, 766]
+            x1 = [447, 854,   847, 488, 412, 199, 675]
             x2 = [92, 647, 821, 422, 33, 159, 369, 99, 157, 464]
             
         rotor_position[x1] += 360
         rotor_position[x2] -= 360
         rotor_position[90] = 180
+        rotor_position[270:276] = 15
+        rotor_position[598:602] = [360,360,0,0]
+        rotor_position[748:752] = [360,360,0,0]
         
         indexes = revolution_trigger(rotor_position, 20,1/(90/20/60))
-        np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
         if 0:
             import matplotlib.pyplot as plt
             plt.plot(rotor_position)
             plt.plot(indexes, np.zeros_like(indexes),'.')
             plt.show()
+        np.testing.assert_array_equal(indexes, [ 91, 180, 270, 360, 450, 540, 630, 720, 810])
+        
+    def test_revolution_trigger_max_rev_diff(self):
+        rotor_position = np.arange(0.,360*10,4)
+        rotor_position += np.random.random(len(rotor_position))
+        rotor_position = rotor_position % 360
+        rotor_speed = 1/(90/20/60)
+        revolution_trigger(rotor_position, 20, rotor_speed*1.02, max_rev_diff=3)
+        self.assertRaisesRegex(AssertionError, "No of revolution mismatch", revolution_trigger, rotor_position, 20, rotor_speed*1.02, max_rev_diff=1)
+        
         
-
 if __name__ == "__main__":
     unittest.main()

wetb/wind/__init__.py

+from .turbulence import *
\ No newline at end of file

wetb/wind/turbulence/mann_parameters.py

@@ -254,20 +254,18 @@ def fit_ae(sf, u, L, G, min_bin_count=None, plt=False):
     return ae
 
 
-def plot_spectra(ae, L, G, k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, show=True, plt=None):
+def plot_spectra(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_mann_spectra(ae, L, G, "-", mean_u, plt)
-    if show:
-        plt.show()
 
-def _plot_spectra(k1, uu, vv, ww=None, uw=None, mean_u=1, log10_bin_size=.2, plt=None):
+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)
+        plt.semilogx(bk1, bk1 * xx * 10 ** 0 / mean_u ** 2 , '.' + color, label=label)
     plot(buu, 'uu', 'r', plt)
     if (bvv) is not None:
         plot(bvv, 'vv', 'g', plt)
@@ -280,8 +278,9 @@ def plot_mann_spectra(ae, L, G, style='-', u_ref=1, plt=None, spectra=['uu', 'vv
         import matplotlib.pyplot as plt
     mf = 10 ** (np.linspace(-4, 1, 1000))
     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 '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)