Merge branch 'develop' of 134.147.164.251:/data/git/pylot into develop

pylot/core/pick/CharFuns.py

@@ -117,8 +117,8 @@ class CharacteristicFunction(object):
 
     def getTimeArray(self):
         if self.getTime1():
-           incr = self.getARdetStep()
+           incr = self.getARdetStep()[0]
-           self.TimeArray = np.arange(0, len(self.getCF()) * incr[0], incr[0]) + self.getCut()[0] \
+           self.TimeArray = np.arange(0, len(self.getCF()) * incr, incr) + self.getCut()[0] \
                            + self.getTime1() + self.getTime2()  
         else:
            incr = self.getIncrement()
@@ -143,18 +143,31 @@ class CharacteristicFunction(object):
         cutting window
         '''
         if cut is not None:
-            if self.cut[0] == 0:
+            if len(self.orig_data) == 1:
+                if self.cut[0] == 0 and self.cut[1] == 0:
                    start = 0
+                   stop = len(self.orig_data[0])
+                elif self.cut[0] == 0 and self.cut[1] is not 0:
+                   start = 0
+                   stop = self.cut[1] / self.dt
                 else:
                    start = self.cut[0] / self.dt
                    stop = self.cut[1] / self.dt
-            if len(self.orig_data) == 1:
                 zz = self.orig_data.copy()
                 z1 = zz[0].copy()
                 zz[0].data = z1.data[start:stop]
                 data = zz
                 return data
             elif len(self.orig_data) == 2:
+                if self.cut[0] == 0 and self.cut[1] == 0:
+                   start = 0
+                   stop = min([len(self.orig_data[0]), len(self.orig_data[1])])
+                elif self.cut[0] == 0 and self.cut[1] is not 0:
+                   start = 0
+                   stop = self.cut[1] / self.dt
+                else:
+                   start = self.cut[0] / self.dt
+                   stop = self.cut[1] / self.dt
                 hh = self.orig_data.copy()
                 h1 = hh[0].copy()
                 h2 = hh[1].copy()
@@ -163,6 +176,15 @@ class CharacteristicFunction(object):
                 data = hh 
                 return data
             elif len(self.orig_data) == 3:
+                if self.cut[0] == 0 and self.cut[1] == 0:
+                   start = 0
+                   stop = min([len(self.orig_data[0]), len(self.orig_data[1]), len(self.orig_data[2])])
+                elif self.cut[0] == 0 and self.cut[1] is not 0:
+                   start = 0
+                   stop = self.cut[1] / self.dt
+                else:
+                   start = self.cut[0] / self.dt
+                   stop = self.cut[1] / self.dt
                 hh = self.orig_data.copy()
                 h1 = hh[0].copy()
                 h2 = hh[1].copy()
@@ -195,6 +217,9 @@ class AICcf(CharacteristicFunction):
         print 'Calculating AIC ...'
         x = self.getDataArray()
         xnp = x[0].data
+        nn = np.isnan(xnp)
+        if len(nn) > 1:
+           xnp[nn] = 0
         datlen = len(xnp)
         k = np.arange(1, datlen)
         cf = np.zeros(datlen)
@@ -224,6 +249,9 @@ class HOScf(CharacteristicFunction):
 
         x = self.getDataArray(self.getCut())
         xnp =x[0].data
+        nn = np.isnan(xnp)
+        if len(nn) > 1:
+           xnp[nn] = 0
         if self.getOrder() == 3:  # this is skewness
             print 'Calculating skewness ...'
             y = np.power(xnp, 3)
@@ -256,6 +284,9 @@ class HOScf(CharacteristicFunction):
             elif self.getOrder() == 4:
                 LTA[j] = lta / np.power(lta1, 2)
         
+        nn = np.isnan(LTA)
+        if len(nn) > 1:
+           LTA[nn] = 0
         self.cf = LTA
 
 
@@ -266,6 +297,9 @@ class ARZcf(CharacteristicFunction):
         print 'Calculating AR-prediction error from single trace ...'
         x = self.getDataArray(self.getCut())
         xnp = x[0].data
+        nn = np.isnan(xnp)
+        if len(nn) > 1:
+           xnp[nn] = 0
         #some parameters needed
         #add noise to time series
         xnoise = xnp + np.random.normal(0.0, 1.0, len(xnp)) * self.getFnoise() * max(abs(xnp))
@@ -288,6 +322,9 @@ class ARZcf(CharacteristicFunction):
         
         #convert list to numpy array
         cf = np.asarray(cf)
+        nn = np.isnan(cf)
+        if len(nn) > 1:
+           cf[nn] = 0
         self.cf = cf
 
 
@@ -376,6 +413,12 @@ class ARHcf(CharacteristicFunction):
         print 'Calculating AR-prediction error from both horizontal traces ...'
         
         xnp = self.getDataArray(self.getCut())
+        n0 = np.isnan(xnp[0].data)
+        if len(n0) > 1:
+           xnp[0].data[n0] = 0
+        n1 = np.isnan(xnp[1].data)
+        if len(n1) > 1:
+           xnp[1].data[n1] = 0
         
         #some parameters needed
         #add noise to time series
@@ -390,7 +433,7 @@ class ARHcf(CharacteristicFunction):
               
         cf = []
         loopstep = self.getARdetStep()
-        for i in range(ldet + self.getOrder() - 3, tend - lpred + 1, loopstep[1]):
+        for i in range(ldet + self.getOrder() - 1, tend - lpred + 1, loopstep[1]):
             self.arDetH(Xnoise, self.getOrder(), i-ldet, i)
             #AR prediction of waveform using calculated AR coefficients
             self.arPredH(xnp, self.arpara, i + 1, lpred)
@@ -401,6 +444,9 @@ class ARHcf(CharacteristicFunction):
 
         #convert list to numpy array
         cf = np.asarray(cf)
+        nn = np.isnan(cf)
+        if len(nn) > 1:
+           cf[nn] = 0
         self.cf = cf
 
     def arDetH(self, data, order, rind, ldet):
@@ -493,6 +539,15 @@ class AR3Ccf(CharacteristicFunction):
         print 'Calculating AR-prediction error from all 3 components ...'
         
         xnp = self.getDataArray(self.getCut())
+        n0 = np.isnan(xnp[0].data)
+        if len(n0) > 1:
+           xnp[0].data[n0] = 0
+        n1 = np.isnan(xnp[1].data)
+        if len(n1) > 1:
+           xnp[1].data[n1] = 0
+        n2 = np.isnan(xnp[2].data)
+        if len(n2) > 1:
+           xnp[2].data[n2] = 0
         
         #some parameters needed
         #add noise to time series
@@ -508,7 +563,7 @@ class AR3Ccf(CharacteristicFunction):
               
         cf = []
         loopstep = self.getARdetStep()
-        for i in range(ldet + self.getOrder() - 3, tend - lpred + 1, loopstep[1]):
+        for i in range(ldet + self.getOrder() - 1, tend - lpred + 1, loopstep[1]):
             self.arDet3C(Xnoise, self.getOrder(), i-ldet, i)
             #AR prediction of waveform using calculated AR coefficients
             self.arPred3C(xnp, self.arpara, i + 1, lpred)
@@ -520,6 +575,9 @@ class AR3Ccf(CharacteristicFunction):
 
         #convert list to numpy array
         cf = np.asarray(cf)
+        nn = np.isnan(cf)
+        if len(nn) > 1:
+           cf[nn] = 0
         self.cf = cf
 
     def arDet3C(self, data, order, rind, ldet):

pylot/core/pick/Picker.py

@@ -3,16 +3,16 @@
 Created Dec 2014
 Implementation of the picking algorithms published and described in:
 
-Küperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010:
+Kueperkoch, L., Meier, T., Lee, J., Friederich, W., & Egelados Working Group, 2010:
 Automated determination of P-phase arrival times at regional and local distances
 using higher order statistics, Geophys. J. Int., 181, 1159-1170
 
-Küperkoch, L., Meier, T., Brüstle, A., Lee, J., Friederich, W., & Egelados
+Kueperkoch, L., Meier, T., Bruestle, A., Lee, J., Friederich, W., & Egelados
 Working Group, 2012: Automated determination of S-phase arrival times using
 autoregressive prediction: application ot local and regional distances, Geophys. J. Int.,
 188, 687-702.
 
-:author: MAGS2 EP3 working group / Ludger Küperkoch	
+:author: MAGS2 EP3 working group / Ludger Kueperkoch	
 """
 import numpy as np
 import matplotlib.pyplot as plt
@@ -23,7 +23,7 @@ class AutoPicking(object):
     Superclass of different, automated picking algorithms applied on a CF determined
     using AIC, HOS, or AR prediction.
     ''' 
-    def __init__(self, cf, Tslope, aerr, TSNR, PickWindow, peps=None, Tsmooth=None):
+    def __init__(self, cf, Tslope, aerr, TSNR, PickWindow, aus=None, Tsmooth=None, Pick1=None):
         '''
         :param: cf, characteristic function, on which the picking algorithm is applied
         :type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
@@ -41,11 +41,14 @@ class AutoPicking(object):
         :param: PickWindow, length of pick window [s]
         :type: float
 
-        :param: peps, find local minimum at i if aic(i-1)*(1+peps) >= aic(i)
+        :param: aus ("artificial uplift of samples"), find local minimum at i if aic(i-1)*(1+aus) >= aic(i)
         :type: float
 
         :param: Tsmooth, length of moving smoothing window to calculate smoothed CF [s]
         :type: float
+
+        :param: Pick1, initial (prelimenary) onset time, starting point for PragPicker
+        :type: float
         '''
 
         #assert isinstance(cf, CharFuns), "%s is not a CharacteristicFunction object" % str(cf)
@@ -58,8 +61,9 @@ class AutoPicking(object):
         self.setaerr(aerr)
         self.setTSNR(TSNR)
         self.setPickWindow(PickWindow)
-        self.setpeps(peps)
+        self.setaus(aus)
         self.setTsmooth(Tsmooth)
+        self.setpick1(Pick1)
         self.calcPick()
 
     def __str__(self):
@@ -68,15 +72,17 @@ class AutoPicking(object):
         aerr:\t{aerr}\n
         TSNR:\t\t\t{TSNR}\n
         PickWindow:\t{PickWindow}\n
-        peps:\t{peps}\n
+        aus:\t{aus}\n
         Tsmooth:\t{Tsmooth}\n
+        Pick1:\t{Pick1}\n
         '''.format(name=type(self).__name__,
                    Tslope=self.getTslope(),
                    aerr=self.getaerr(),
                    TSNR=self.getTSNR(),
                    PickWindow=self.getPickWindow(),
-                   peps=self.getpeps(),
+                   aus=self.getaus(),
-                   Tsmooth=self.getTsmooth())     
+                   Tsmooth=self.getTsmooth(),
+                   Pick1=self.getpick1())     
 
     def getTslope(self):
         return self.Tslope
@@ -102,11 +108,11 @@ class AutoPicking(object):
     def setPickWindow(self, PickWindow):
         self.PickWindow = PickWindow
 
-    def getpeps(self):
+    def getaus(self):
-        return self.peps
+        return self.aus
 
-    def setpeps(self, peps):
+    def setaus(self, aus):
-        self.peps = peps
+        self.aus = aus
 
     def setTsmooth(self, Tsmooth):
         self.Tsmooth = Tsmooth
@@ -117,6 +123,12 @@ class AutoPicking(object):
     def getpick(self):
         return self.Pick
 
+    def getpick1(self):
+        return self.Pick1
+
+    def setpick1(self, Pick1):
+        self.Pick1 = Pick1
+
     def calcPick(self):
         self.Pick = None
  
@@ -128,7 +140,7 @@ class AICPicker(AutoPicking):
 
     def calcPick(self):
  
-        print 'Get onset (pick) from AIC-CF ...'
+        print 'Get onset time (pick) from AIC-CF ...'
 
         self.Pick = -1
         #taper AIC-CF to get rid off side maxima
@@ -155,4 +167,78 @@ class PragPicker(AutoPicking):
 
     def calcPick(self):
 
-        print 'Get onset (pick) from HOS- or AR-CF using pragmatic picking algorithm ...'
+        if self.getpick1() is not None:
+           print 'Get onset time (pick) from HOS- or AR-CF using pragmatic picking algorithm ...'
+
+           self.Pick = -1
+           #smooth CF
+           ismooth = round(self.Tsmooth / self.dt);
+           cfsmooth = np.zeros(len(self.cf))
+           if len(self.cf) < ismooth:
+              print 'PragPicker: Tsmooth larger than CF!'
+              return
+           else:
+              for i in range(1, len(self.cf)):
+                 if i > ismooth:
+                    ii1 = i - ismooth;
+                    cfsmooth[i] = cfsmooth[i - 1] + (self.cf[i] - self.cf[ii1]) / ismooth
+                 else:
+                    cfsmooth[i] = np.mean(self.cf[1 : i])
+
+           #select picking window
+           #which is centered around tpick1
+           ipick = np.where((self.Tcf >= self.getpick1() - self.PickWindow / 2) \
+                 & (self.Tcf <= self.getpick1() + self.PickWindow / 2))
+           cfipick = self.cf[ipick]
+           Tcfpick = self.Tcf[ipick]
+           cfsmoothipick = cfsmooth[ipick]
+           ipick1 = np.argmin(abs(self.Tcf - self.getpick1()))
+           cfpick1 = 2 * self.cf[ipick1]
+
+           #check trend of CF, i.e. differences of CF and adjust aus regarding this trend
+           #prominent trend: decrease aus
+           #flat: use given aus
+           cfdiff = np.diff(cfipick);
+           i0diff = np.where(cfdiff > 0)
+           cfdiff = cfdiff[i0diff]
+           minaus = min(cfdiff * (1 + self.aus));
+           aus1 = max([minaus, self.aus]);
+
+           #at first we look to the right until the end of the pick window is reached
+           flagpick_r = 0
+           flagpick_l = 0
+           flagpick = 0
+           lpickwindow = int(round(self.PickWindow / self.dt))
+           for i in range(max(np.insert(ipick, 0, 2)), min([ipick1 + lpickwindow + 1, len(self.cf) - 1])):
+              if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]:
+                 if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
+                    if cfpick1 >= self.cf[i]:     
+                       pick_r = self.Tcf[i]
+                       self.Pick = pick_r
+                       flagpick_l = 1
+                       cfpick_r = self.cf[i]
+                       break
+
+           #now we look to the left
+           for i in range(ipick1, max([ipick1 - lpickwindow + 1, 2]), -1):
+              if self.cf[i + 1] > self.cf[i] and self.cf[i - 1] >= self.cf[i]:
+                 if cfsmooth[i - 1] * (1 + aus1) >= cfsmooth[i]:
+                    if cfpick1 >= self.cf[i]:     
+                       pick_l = self.Tcf[i]
+                       self.Pick = pick_l
+                       flagpick_r = 1
+                       cfpick_l = self.cf[i]
+                       break
+
+           #now decide which pick: left or right?
+           if flagpick_l > 0 and flagpick_r > 0:
+              if cfpick_l <= cfpick_r:
+                 self.Pick = pick_l
+              else:
+                 self.Pick = pick_r
+
+        else: 
+           self.Pick = -1
+           print 'PragPicker: No initial onset time given! Check input!'
+           return
+

pylot/core/pick/run_makeCF.py

@@ -19,7 +19,7 @@ def run_makeCF(project, database, event, iplot, station=None):
     #parameters for CF calculation
     t2 = 7              #length of moving window for HOS calculation [sec]
     p = 4               #order of statistics
-    cuttimes = [10, 40] #start and end time vor CF calculation
+    cuttimes = [5, 40] #start and end time for CF calculation
     bpz = [2, 30]       #corner frequencies of bandpass filter, vertical component
     bph = [2, 15]       #corner frequencies of bandpass filter, horizontal components
     tdetz= 1.2          #length of AR-determination window [sec], vertical component
@@ -59,9 +59,6 @@ def run_makeCF(project, database, event, iplot, station=None):
           #calculate HOS-CF using subclass HOScf of class CharacteristicFunction
           hoscf = HOScf(st_copy, cuttimes, t2, p) #instance of HOScf
           ##############################################################
-          #get onset time from HOS-CF using class Picker
-          #hospick = PragPicker(hoscf, 2, 70, [1, 0.5, 0.2], 2, 0.001, 0.2)
-          ##############################################################
           #calculate AIC-HOS-CF using subclass AICcf of class CharacteristicFunction
           #class needs stream object => build it
           tr_aic = tr_filt.copy()
@@ -69,9 +66,12 @@ def run_makeCF(project, database, event, iplot, station=None):
           st_copy[0].data = tr_aic.data
           aiccf = AICcf(st_copy, cuttimes, t2) #instance of AICcf
           ##############################################################
-          #get onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking
+          #get prelimenary onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking
           aicpick = AICPicker(aiccf, 2, 70, [1, 0.5, 0.2], 3)
           ##############################################################
+          #get refined onset time from HOS-CF using class Picker
+          hospick = PragPicker(hoscf, 2, 70, [1, 0.5, 0.2], 2, 0.001, 0.2, aicpick.getpick())
+          ##############################################################
           #calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction
           #get stream object of filtered data
           st_copy[0].data = tr_filt.data
@@ -86,6 +86,9 @@ def run_makeCF(project, database, event, iplot, station=None):
           ##############################################################
           #get onset time from AIC-ARZ-CF using subclass AICPicker of class AutoPicking
           aicarzpick = AICPicker(araiccf,  2, 70, [1, 0.5, 0.2], 2)
+          ##############################################################
+          #get refined onset time from ARZ-CF using class Picker
+          arzpick = PragPicker(arzcf, 2, 70, [1, 0.5, 0.2], 2, 0, 0.2, aicarzpick.getpick())
     elif not wfzfiles:
        print 'No vertical component data found!'
 
@@ -156,9 +159,15 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([aicpick.getpick(), aicpick.getpick()], [-1, 1], 'b--')
              plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [1, 1], 'b')
              plt.plot([aicpick.getpick()-0.5, aicpick.getpick()+0.5], [-1, -1], 'b')
+             plt.plot([hospick.getpick(), hospick.getpick()], [-1.3, 1.3], 'r--')
+             plt.plot([hospick.getpick()-0.5, hospick.getpick()+0.5], [1.3, 1.3], 'r')
+             plt.plot([hospick.getpick()-0.5, hospick.getpick()+0.5], [-1.3, -1.3], 'r')
              plt.plot([aicarzpick.getpick(), aicarzpick.getpick()], [-1.2, 1.2], 'y--')
              plt.plot([aicarzpick.getpick()-0.5, aicarzpick.getpick()+0.5], [1.2, 1.2], 'y')
              plt.plot([aicarzpick.getpick()-0.5, aicarzpick.getpick()+0.5], [-1.2, -1.2], 'y')
+             plt.plot([arzpick.getpick(), arzpick.getpick()], [-1.4, 1.4], 'g--')
+             plt.plot([arzpick.getpick()-0.5, arzpick.getpick()+0.5], [1.4, 1.4], 'g')
+             plt.plot([arzpick.getpick()-0.5, arzpick.getpick()+0.5], [-1.4, -1.4], 'g')
              plt.yticks([])
              plt.xlabel('Time [s]')
              plt.ylabel('Normalized Counts')
@@ -173,7 +182,7 @@ def run_makeCF(project, database, event, iplot, station=None):
              th2data = np.arange(0, trH2_filt.stats.npts / trH2_filt.stats.sampling_rate, trH2_filt.stats.delta)
              tarhcf = np.arange(0, len(arhcf.getCF()) * tsteph, tsteph) + cuttimes[0] + tdeth +tpredh
              p21 = plt.plot(th1data, trH1_filt.data/max(trH1_filt.data), 'k')
-             p22 = plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r') 
+             p22 = plt.plot(arhcf.getTimeArray(), arhcf.getCF()/max(arhcf.getCF()), 'r') 
              p23 = plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF()))
              plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b--')
              plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b')
@@ -185,7 +194,6 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.legend([p21, p22, p23], ['Data', 'ARH-CF', 'ARHAIC-CF']) 
              plt.subplot(212)
              plt.plot(th2data, trH2_filt.data/max(trH2_filt.data), 'k')
-             plt.plot(tarhcf, arhcf.getCF()/max(arhcf.getCF()), 'r') 
              plt.plot(arhaiccf.getTimeArray(), arhaiccf.getCF()/max(arhaiccf.getCF()))
              plt.plot([aicarhpick.getpick(), aicarhpick.getpick()], [-1, 1], 'b--')
              plt.plot([aicarhpick.getpick()-0.5, aicarhpick.getpick()+0.5], [1, 1], 'b')