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

pylot/core/pick/Picker.py

@@ -27,15 +27,11 @@ class AutoPicking(object):
     Superclass of different, automated picking algorithms applied on a CF determined
     using AIC, HOS, or AR prediction.
     ''' 
-    def __init__(self, cf, nfac, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, Pick1=None):
+    def __init__(self, cf, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, Pick1=None):
         '''
         :param: cf, characteristic function, on which the picking algorithm is applied
         :type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
 
-        :param: nfac (noise factor), nfac times noise level to calculate latest possible pick
-         in EarlLatePicker
-        :type: int
-
         :param: TSNR, length of time windows around pick used to determine SNR [s]
         :type: tuple (T_noise, T_gap, T_signal)
 
@@ -63,7 +59,6 @@ class AutoPicking(object):
         self.Tcf = cf.getTimeArray()
         self.Data = cf.getXCF()
         self.dt = cf.getIncrement()
-        self.setnfac(nfac)
         self.setTSNR(TSNR)
         self.setPickWindow(PickWindow)
         self.setiplot(iplot)
@@ -74,25 +69,18 @@ class AutoPicking(object):
 
     def __str__(self):
         return '''\n\t{name} object:\n
-        nfac:\t{nfac}\n
         TSNR:\t\t\t{TSNR}\n
         PickWindow:\t{PickWindow}\n
         aus:\t{aus}\n
         Tsmooth:\t{Tsmooth}\n
         Pick1:\t{Pick1}\n
         '''.format(name=type(self).__name__,
-                   nfac=self.getnfac(),
                    TSNR=self.getTSNR(),
                    PickWindow=self.getPickWindow(),
                    aus=self.getaus(),
                    Tsmooth=self.getTsmooth(),
                    Pick1=self.getpick1())     
 
-    def getnfac(self):
-        return self.nfac 
-              
-    def setnfac(self, nfac):
-        self.nfac = nfac
  
     def getTSNR(self):
         return self.TSNR
@@ -127,15 +115,6 @@ class AutoPicking(object):
     def getSlope(self):
         return self.slope
 
-    def getLpick(self):
-        return self.LPick
-
-    def getEpick(self):
-        return self.EPick
-
-    def getPickError(self):
-        return self.PickError
-
     def getiplot(self):
         return self.iplot
 
@@ -165,7 +144,6 @@ class AICPicker(AutoPicking):
         print 'AICPicker: Get initial onset time (pick) from AIC-CF ...'
 
         self.Pick = None
-        self.PickError = None
         #find NaN's
         nn = np.isnan(self.cf)
         if len(nn) > 1:
@@ -263,7 +241,7 @@ class AICPicker(AutoPicking):
            p1, = plt.plot(self.Tcf, x / max(x), 'k')
            p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r')
            if self.Pick is not None:
-              p3, = plt.plot([self.Pick, self.Pick], [-1 , 1], 'b', linewidth=2)
+              p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2)
               plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
            else:
               plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
@@ -281,7 +259,7 @@ class AICPicker(AutoPicking):
               p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2)
               plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], \
                          loc='best')
-              plt.title('SNR and Slope, Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \
+              plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \
                         self.SNR, self.slope))
               plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
               plt.ylabel('Counts')
@@ -307,7 +285,6 @@ class PragPicker(AutoPicking):
            print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...'
 
            self.Pick = None 
-           self.PickError = None
            self.SNR = None
            self.slope = None
            #smooth CF
@@ -392,313 +369,3 @@ class PragPicker(AutoPicking):
            self.Pick = None
            print 'PragPicker: No initial onset time given! Check input!'
            return
-
-
-class EarlLatePicker(AutoPicking):
-    '''
-    Method to derive earliest and latest possible pick after Diehl & Kissling (2009) 
-    as reasonable uncertainties. Latest possible pick is based on noise level,
-    earliest possible pick is half a signal wavelength in front of most likely 
-    pick given by PragPicker. Most likely pick (initial pick) must be given. 
-    '''
-       
-    def calcPick(self):
-
-        self.LPick = None 
-        self.EPick = None
-        self.PickError = None
-        self.SNR = None
-        self.slope = None
-        if self.getpick1() is not None:
-           print 'EarlLatePicker: Get earliest and latest possible pick relative to most likely pick ...'
-
-           ti = self.getpick1()
-           x = self.Data
-           t = self.Tcf
-           #some parameters needed:
-           tnoise = self.TSNR[0]  #noise window length for calculating noise level
-           tsignal = self.TSNR[2] #signal window length
-           tsafety = self.TSNR[1] #safety gap between signal onset and noise window
-
-           #get latest possible pick
-           #get noise window
-           inoise = np.where((self.Tcf <= max([ti - tsafety, 0])) \
-                    & (self.Tcf >= max([ti - tnoise - tsafety, 0])))    
-           #get signal window
-           isignal = np.where((self.Tcf <= min([ti + tsignal + tsafety, len(x[0].data)])) \
-                     & (self.Tcf >= ti))
-           #calculate noise level
-           if len(x) == 1:
-              nlevel = max(abs(x[0].data[inoise])) * self.nfac
-              #get time where signal exceeds nlevel
-              ilup = np.where(x[0].data[isignal] > nlevel)
-              ildown = np.where(x[0].data[isignal] < -nlevel)
-              if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
-                 print 'EarlLatePicker: Signal lower than noise level, misspick?'
-                 return
-              il = min([ilup[0][0], ildown[0][0]])
-              self.LPick = t[isignal][il]
-           elif len(x) == 2:
-              nlevel = max(np.sqrt(np.power(x[0].data[inoise], 2) + np.power(x[1].data[inoise], 2)))
-              #get earliest time where signal exceeds nlevel
-              ilup1 = np.where(x[0].data[isignal] > nlevel)
-              ilup2 = np.where(x[1].data[isignal] > nlevel)
-              ildown1 = np.where(x[0].data[isignal] < -nlevel)
-              ildown2 = np.where(x[1].data[isignal] < -nlevel)
-              if np.size(ilup1) < 1 and np.size(ilup2) > 1:
-                 ilup = ilup2
-              elif np.size(ilup1) > 1 and np.size(ilup2) < 1:
-                 ilup = ilup1
-              elif np.size(ilup1) < 1 and np.size(ilup2) < 1:
-                 ilup = None
-              else:
-                 ilup = min([ilup1[0][0], ilup2[0][0]])
- 
-              if np.size(ildown1) < 1 and np.size(ildown2) > 1:
-                 ildown = ildown2
-              elif np.size(ildown1) > 1 and np.size(ildown2) < 1:
-                 ildown = ildown1
-              elif np.size(ildown1) < 1 and np.size(ildown2) < 1:
-                 ildown = None
-              else:
-                 ildown = min([ildown1[0][0], ildown2[0][0]])
-              if ilup == None  and ildown == None:
-                 print 'EarlLatePicker: Signal lower than noise level, misspick?'
-                 return
-              il = min([ilup, ildown])
-              self.LPick = t[isignal][il]
-           elif len(x) == 3:
-              nlevel = max(np.sqrt(np.power(x[0].data[inoise], 2) + np.power(x[1].data[inoise], 2) + \
-                       np.power(x[2].data[inoise], 2)))
-              #get earliest time where signal exceeds nlevel
-              ilup1 = np.where(x[0].data[isignal] > nlevel)
-              ilup2 = np.where(x[1].data[isignal] > nlevel)
-              ilup3 = np.where(x[2].data[isignal] > nlevel)
-              ildown1 = np.where(x[0].data[isignal] < -nlevel)
-              ildown2 = np.where(x[1].data[isignal] < -nlevel)
-              ildown3 = np.where(x[2].data[isignal] < -nlevel)
-              if np.size(ilup1) > 1 and np.size(ilup2) < 1 and np.size(ilup3) < 1:
-                 ilup = ilup1
-              elif np.size(ilup1) > 1 and np.size(ilup2) > 1 and np.size(ilup3) < 1:
-                 ilup = min([ilup1[0][0], ilup2[0][0]])
-              elif np.size(ilup1) > 1 and np.size(ilup2) > 1 and np.size(ilup3) > 1: 
-                 ilup = min([ilup1[0][0], ilup2[0][0], ilup3[0][0]])
-              elif np.size(ilup1) < 1 and np.size(ilup2) > 1 and np.size(ilup3) > 1:
-                 ilup = min([ilup2[0][0], ilup3[0][0]])
-              elif np.size(ilup1) > 1 and np.size(ilup2) < 1 and np.size(ilup3) > 1:
-                 ilup = min([ilup1[0][0], ilup3[0][0]])
-              elif np.size(ilup1) < 1 and np.size(ilup2) < 1 and np.size(ilup3) < 1:
-                 ilup = None
-              else:
-                 ilup = min([ilup1[0][0], ilup2[0][0], ilup3[0][0]])
-
-              if np.size(ildown1) > 1 and np.size(ildown2) < 1 and np.size(ildown3) < 1:
-                 ildown = ildown1
-              elif np.size(ildown1) > 1 and np.size(ildown2) > 1 and np.size(ildown3) < 1:
-                 ildown = min([ildown1[0][0], ildown2[0][0]])
-              elif np.size(ildown1) > 1 and np.size(ildown2) > 1 and np.size(ildown3) > 1: 
-                 ildown = min([ildown1[0][0], ildown2[0][0], ildown3[0][0]])
-              elif np.size(ildown1) < 1 and np.size(ildown2) > 1 and np.size(ildown3) > 1:
-                 ildown = min([ildown2[0][0], ildown3[0][0]])
-              elif np.size(ildown1) > 1 and np.size(ildown2) < 1 and np.size(ildown3) > 1:
-                 ildown = min([ildown1[0][0], ildown3[0][0]])
-              elif np.size(ildown1) < 1 and np.size(ildown2) < 1 and np.size(ildown3) < 1:
-                 ildown = None
-              else:
-                 ildown = min([ildown1[0][0], ildown2[0][0], ildown3[0][0]])
-              if ilup == None and ildown == None:
-                 print 'EarlLatePicker: Signal lower than noise level, misspick?'
-                 return
-              il = min([ilup, ildown])
-              self.LPick = t[isignal][il]
-
-           #get earliest possible pick
-           #get next 2 zero crossings after most likely pick
-           #if there is one trace in stream
-           if len(x) == 1:
-              zc = []
-              zc.append(ti)
-              i = 0
-              for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
-                  i = i+ 1
-                  if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
-                     zc.append(t[isignal][i])
-                  elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
-                     zc.append(t[isignal][i])
-                  if len(zc) == 3:
-                     break
-              #calculate maximum period of signal out of zero crossings
-              Ts = max(np.diff(zc))
-           #if there are two traces in stream
-           #get maximum of two signal periods
-           if len(x) == 2:
-              zc1 = []
-              zc2 = []
-              zc1.append(ti)
-              zc2.append(ti)
-              i = 0
-              for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
-                  i = i+ 1
-                  if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
-                     zc1.append(t[isignal][i])
-                  elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
-                     zc1.append(t[isignal][i])
-                  if x[1].data[j-1] <= 0 and x[1].data[j] >= 0:
-                     zc2.append(t[isignal][i])
-                  elif x[1].data[j-1] > 0 and x[1].data[j] <= 0:
-                     zc2.append(t[isignal][i])
-                  if len(zc1) >= 3 and len(zc2) >= 3:
-                     break
-              Ts = max([max(np.diff(zc1)), max(np.diff(zc2))])
-           #if there are three traces in stream
-           #get maximum of three signal periods
-           if len(x) == 3:
-              zc1 = []
-              zc2 = []
-              zc3 = []
-              zc1.append(ti)
-              zc2.append(ti)
-              zc3.append(ti)
-              i = 0
-              for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
-                  i = i+ 1
-                  if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
-                     zc1.append(t[isignal][i])
-                  elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
-                     zc1.append(t[isignal][i])
-                  if x[1].data[j-1] <= 0 and x[1].data[j] >= 0:
-                     zc2.append(t[isignal][i])
-                  elif x[1].data[j-1] > 0 and x[1].data[j] <= 0:
-                     zc2.append(t[isignal][i])
-                  if x[2].data[j-1] <= 0 and x[2].data[j] >= 0:
-                     zc3.append(t[isignal][i])
-                  elif x[2].data[j-1] > 0 and x[2].data[j] <= 0:
-                     zc3.append(t[isignal][i])
-                  if len(zc1) >= 3 and len(zc2) >= 3 and len(zc3) >= 3:
-                     break
-              Ts = max([max(np.diff(zc1)), max(np.diff(zc2)), max(np.diff(zc3))])
-
-           #Ts/4 is assumed as time difference between most likely and earliest possible pick!
-           self.EPick = ti - Ts/4
-
-           #get symmetric pick error as mean from earliest and latest possible pick
-           #by weighting latest possible pick tow times earliest possible pick
-           diffti_tl = self.LPick - ti
-           diffti_te = ti - self.EPick
-           self.PickError = (diffti_te + 2 * diffti_tl)  / 3
-
-           if self.iplot is not None:
-              plt.figure(self.iplot)
-              if len(x) == 1:
-                 p1, = plt.plot(t, x[0].data, 'k')
-                 p2, = plt.plot(t[inoise], x[0].data[inoise])
-                 p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
-                 p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14)
-                 plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
-                             loc='best')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
-                 plt.yticks([])
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-                 plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
-              elif len(x) == 2:
-                 plt.subplot(2,1,1)
-                 p1, = plt.plot(t, x[0].data, 'k')
-                 p2, = plt.plot(t[inoise], x[0].data[inoise])
-                 p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
-                 p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 p5, = plt.plot(zc1[0:3], [0, 0, 0], '*g', markersize=14)
-                 plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
-                             loc='best')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot(zc1[0:3], [0, 0, 0], '*g')
-                 plt.yticks([])
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-                 plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
-                 plt.subplot(2,1,2)
-                 plt.plot(t, x[1].data, 'k')
-                 plt.plot(t[inoise], x[1].data[inoise])
-                 plt.plot(t[isignal], x[1].data[isignal], 'r')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[1].data), -max(x[1].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot(zc2[0:3], [0, 0, 0], '*g', markersize=14)
-                 plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-                 plt.yticks([])
-              elif len(x) == 3:
-                 plt.subplot(3,1,1)
-                 p1, = plt.plot(t, x[0].data, 'k')
-                 p2, = plt.plot(t[inoise], x[0].data[inoise])
-                 p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
-                 p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 p5, = plt.plot(zc1[0:3], [0, 0, 0], '*g', markersize=14)
-                 plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
-                             loc='best')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.yticks([])
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-                 plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
-                 plt.subplot(3,1,2)
-                 plt.plot(t, x[1].data, 'k')
-                 plt.plot(t[inoise], x[1].data[inoise])
-                 plt.plot(t[isignal], x[1].data[isignal], 'r')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[1].data), -max(x[1].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot(zc2[0:3], [0, 0, 0], '*g', markersize=14)
-                 plt.yticks([])
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-                 plt.subplot(3,1,3)
-                 plt.plot(t, x[2].data, 'k')
-                 plt.plot(t[inoise], x[2].data[inoise])
-                 plt.plot(t[isignal], x[2].data[isignal], 'r')
-                 plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
-                 plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
-                 plt.plot([ti, ti], [max(x[2].data), -max(x[2].data)], 'b', linewidth=2)
-                 plt.plot([self.LPick, self.LPick], [max(x[2].data)/2, -max(x[2].data)/2], '--k')
-                 plt.plot([self.EPick, self.EPick], [max(x[2].data)/2, -max(x[2].data)/2], '--k')
-                 plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
-                 plt.plot(zc3[0:3], [0, 0, 0], '*g', markersize=14)
-                 plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
-                 plt.yticks([])
-                 ax = plt.gca()
-                 ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
-              plt.show()
-              raw_input()
-              plt.close(self.iplot)
-
-        elif self.getpick1() == None: 
-           print 'EarlLatePicker: No initial onset time given! Check input!'
-           return
-

pylot/core/pick/earllatepicker.py

@@ -0,0 +1,138 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+   Created Mar 2015
+   Transcription of the rezipe of Diehl et al. (2009) for consistent phase
+   picking. For a given inital (the most likely) pick, the corresponding earliest
+   and latest possible picks are calculated based on noise measurements in front of
+   the most likely pick and signal wavelength derived from zero crossings.
+
+   :author: MAGS2 EP3 working group / Ludger Kueperkoch	
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from obspy.core import Stream
+import argparse
+
+def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
+    '''
+    Function to derive earliest and latest possible pick after Diehl & Kissling (2009) 
+    as reasonable uncertainties. Latest possible pick is based on noise level,
+    earliest possible pick is half a signal wavelength in front of most likely 
+    pick given by PragPicker. Most likely pick (initial pick) must be given. 
+    
+    :param: x, time series (seismogram)
+    :type:  `~obspy.core.stream.Stream`
+   
+    :param: nfac (noise factor), nfac times noise level to calculate latest possible pick
+            in EarlLatePicker
+    :type: int
+
+    :param: TSNR, length of time windows around pick used to determine SNR [s]
+    :type: tuple (T_noise, T_gap, T_signal)
+ 
+    :param: Pick1, initial (prelimenary) onset time, starting point for EarlLatePicker
+    :type: float
+
+    '''
+
+    assert isinstance(X, Stream), "%s is not a stream object" % str(X)
+
+    LPick = None 
+    EPick = None
+    PickError = None
+    if Pick1 is not None:
+       print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...'
+
+       x =X[0].data
+       t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta)
+       #some parameters needed:
+       tnoise = TSNR[0]  #noise window length for calculating noise level
+       tsignal = TSNR[2] #signal window length
+       tsafety = TSNR[1] #safety gap between signal onset and noise window
+
+       #get latest possible pick
+       #get noise window
+       inoise = np.where((t <= max([Pick1 - tsafety, 0])) \
+                & (t >= max([Pick1 - tnoise - tsafety, 0])))    
+       #get signal window
+       isignal = np.where((t <= min([Pick1 + tsignal + tsafety, len(x)])) \
+                 & (t >= Pick1))
+       #calculate noise level
+       nlevel = max(abs(x[inoise])) * nfac
+       #get time where signal exceeds nlevel
+       ilup = np.where(x[isignal] > nlevel)
+       ildown = np.where(x[isignal] < -nlevel)
+       if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
+          print 'earllatepicker: Signal lower than noise level, misspick?'
+          return
+       il = min([ilup[0][0], ildown[0][0]])
+       LPick = t[isignal][il]
+
+       #get earliest possible pick
+       #get next 2 zero crossings after most likely pick
+       #if there is one trace in stream
+       zc = []
+       zc.append(Pick1)
+       i = 0
+       for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
+           i = i+ 1
+           if x[j-1] <= 0 and x[j] >= 0:
+              zc.append(t[isignal][i])
+           elif x[j-1] > 0 and x[j] <= 0:
+                zc.append(t[isignal][i])
+           if len(zc) == 3:
+               break
+       #calculate maximum period of signal out of zero crossings
+       Ts = max(np.diff(zc))
+       #Ts/4 is assumed as time difference between most likely and earliest possible pick!
+       EPick = Pick1 - Ts/4
+
+       #get symmetric pick error as mean from earliest and latest possible pick
+       #by weighting latest possible pick tow times earliest possible pick
+       diffti_tl = LPick -Pick1 
+       diffti_te = Pick1 - EPick
+       PickError = (diffti_te + 2 * diffti_tl)  / 3
+
+       if iplot is not None:
+          plt.figure(iplot)
+          p1, = plt.plot(t, x, 'k')
+          p2, = plt.plot(t[inoise], x[inoise])
+          p3, = plt.plot(t[isignal], x[isignal], 'r')
+          p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') 
+          p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14)
+          plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
+                      loc='best')
+          plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') 
+          plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2)
+          plt.plot([LPick, LPick], [max(x)/2, -max(x)/2], '--k')
+          plt.plot([EPick, EPick], [max(x)/2, -max(x)/2], '--k')
+          plt.plot([Pick1 + PickError, Pick1 + PickError], [max(x)/2, -max(x)/2], 'r--')
+          plt.plot([Pick1 - PickError, Pick1 - PickError], [max(x)/2, -max(x)/2], 'r--')
+          plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
+          plt.yticks([])
+          ax = plt.gca()
+          ax.set_xlim([t[inoise[0][0]] - 2, t[isignal[0][len(isignal) - 1]] + 3])
+          plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % X[0].stats.station)
+          plt.show()
+          raw_input()
+          plt.close(iplot)
+
+    elif Pick1 == None: 
+         print 'earllatepicker: No initial onset time given! Check input!'
+         return
+
+    return EPick, LPick, PickError
+
+if __name__ == "__main__":
+   parser = argparse.ArgumentParser()
+   parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read')
+   parser.add_argument('--nfac', type=int, help='(noise factor), nfac times noise level to calculate latest possible pick')
+   parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \
+                       [s] (Tnoise, Tgap, Tsignal)')
+   parser.add_argument('--Pick1', type=float, help='Onset time of most likely pick')
+   parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
+   args = parser.parse_args()
+   earllatepicker(args.X, args.nfac, args.TSNR, args.Pick1, args.iplot)
+
+#earllatepicker(X, nfac, TSNR, Pick1, iplot)

pylot/core/pick/fmpicker.py

@@ -0,0 +1,183 @@
+#!/usr/bin/python
+# -*- coding: utf-8 -*-
+"""
+   Created Mar 2015
+   Function to derive first motion (polarity) for given phase onset based on zero crossings.
+
+   :author: MAGS2 EP3 working group / Ludger Kueperkoch	
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from obspy.core import Stream
+import argparse
+
+def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
+    '''
+    Function to derive first motion (polarity) of given phase onset Pick.
+    Calculation is based on zero crossings determined within time window pickwin
+    after given onset time. 
+    
+    :param: Xraw, unfiltered time series (seismogram)
+    :type:  `~obspy.core.stream.Stream`
+   
+    :param: Xfilt, filtered time series (seismogram)
+    :type:  `~obspy.core.stream.Stream`
+
+    :param: pickwin, time window after onset Pick within zero crossings are calculated
+    :type: float
+
+    :param: Pick, initial (most likely) onset time, starting point for fmpicker
+    :type: float
+
+    :param: iplot, if given, results are plotted in figure(iplot)
+    :type: int
+    '''
+
+    assert isinstance(Xraw, Stream), "%s is not a stream object" % str(Xraw)
+    assert isinstance(Xfilt, Stream), "%s is not a stream object" % str(Xfilt)
+
+    FM = None
+    if Pick is not None:
+       print 'fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...'
+
+       xraw = Xraw[0].data
+       xfilt = Xfilt[0].data
+       t = np.arange(0, Xraw[0].stats.npts / Xraw[0].stats.sampling_rate, Xraw[0].stats.delta)
+       #get pick window
+       ipick = np.where((t <= min([Pick + pickwin, len(Xraw[0])])) & (t >= Pick))
+       #remove mean
+       xraw[ipick] = xraw[ipick] - np.mean(xraw[ipick])
+       xfilt[ipick] = xfilt[ipick] - np.mean(xfilt[ipick])
+
+       #get next zero crossing after most likely pick
+       #initial onset is assumed to be the first zero crossing
+       #first from unfiltered trace
+       zc1 = []
+       zc1.append(Pick)
+       index1 = []
+       i = 0
+       for j in range(ipick[0][1],ipick[0][len(t[ipick]) - 1]):
+           i = i+ 1
+           if xraw[j-1] <= 0 and xraw[j] >= 0:
+              zc1.append(t[ipick][i])
+              index1.append(i)
+           elif xraw[j-1] > 0 and xraw[j] <= 0:
+              zc1.append(t[ipick][i])
+              index1.append(i)
+           if len(zc1) == 3:
+               break
+
+       #if time difference betweeen 1st and 2cnd zero crossing
+       #is too short, get time difference between 1st and 3rd
+       #to derive maximum 
+       if zc1[1] - zc1[0] <= Xraw[0].stats.delta:
+          li1 = index1[1]
+       else:
+          li1 = index1[0]
+       if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0:
+          print 'earllatepicker: Onset on unfiltered trace too emergent for first motion determination!'
+          P1 = None
+       else: 
+          imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]]))
+          islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1]))
+          #calculate slope as polynomal fit of order 1
+          xslope1 = np.arange(0, len(xraw[islope1]), 1)
+          P1 = np.polyfit(xslope1, xraw[islope1], 1)
+          datafit1 = np.polyval(P1, xslope1)
+
+       #now using filterd trace
+       #next zero crossing after most likely pick
+       zc2 = []
+       zc2.append(Pick)
+       index2 = []
+       i = 0
+       for j in range(ipick[0][1],ipick[0][len(t[ipick]) - 1]):
+           i = i+ 1
+           if xfilt[j-1] <= 0 and xfilt[j] >= 0:
+              zc2.append(t[ipick][i])
+              index2.append(i)
+           elif xfilt[j-1] > 0 and xfilt[j] <= 0:
+              zc2.append(t[ipick][i])
+              index2.append(i)
+           if len(zc2) == 3:
+               break
+
+       #if time difference betweeen 1st and 2cnd zero crossing
+       #is too short, get time difference between 1st and 3rd
+       #to derive maximum 
+       if zc2[1] - zc2[0] <= Xfilt[0].stats.delta:
+          li2 = index2[1]
+       else:
+          li2 = index2[0]
+       if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0:
+          print 'earllatepicker: Onset on filtered trace too emergent for first motion determination!'
+          P2 = None
+       else: 
+          imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]]))
+          islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2]))
+          #calculate slope as polynomal fit of order 1
+          xslope2 = np.arange(0, len(xfilt[islope2]), 1)
+          P2 = np.polyfit(xslope2, xfilt[islope2], 1)
+          datafit2 = np.polyval(P2, xslope2)
+       
+       #compare results
+       if P1 is not None and P2 is not None:
+          if P1[0] < 0 and P2[0] < 0:
+             FM = 'D'
+          elif P1[0] >= 0 and P2[0] < 0:
+             FM = '-'
+          elif P1[0] < 0 and P2[0]>= 0:
+             FM = '-'
+          elif P1[0] > 0 and P2[0] > 0:
+             FM = 'U'
+          elif P1[0] <= 0 and P2[0] > 0:
+             FM = '+'
+          elif P1[0] > 0 and P2[0] <= 0:
+             FM = '+'
+
+    if iplot is not None:
+       plt.figure(iplot)
+       plt.subplot(2,1,1)
+       plt.plot(t, xraw, 'k')
+       p1, = plt.plot([Pick, Pick], [max(xraw), -max(xraw)], 'b', linewidth=2)
+       if P1 is not None:
+          p2, = plt.plot(t[islope1], xraw[islope1])
+          p3, = plt.plot(zc1, np.zeros(len(zc1)), '*g', markersize=14)
+          p4, = plt.plot(t[islope1], datafit1, '--g', linewidth=2) 
+          plt.legend([p1, p2, p3, p4], ['Pick', 'Slope Window', 'Zero Crossings', 'Slope'], \
+                      loc='best')
+          plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14) 
+          ax = plt.gca()
+          ax.set_xlim([t[islope1[0][0]] - 0.1, t[islope1[0][len(islope1) - 1]] + 0.3])
+       plt.yticks([])
+       plt.title('First-Motion Determination, %s, Unfiltered Data' % Xraw[0].stats.station)
+
+       plt.subplot(2,1,2)
+       plt.title('First-Motion Determination, Filtered Data')
+       plt.plot(t, xfilt, 'k')
+       p1, = plt.plot([Pick, Pick], [max(xfilt), -max(xfilt)], 'b', linewidth=2)
+       if P2 is not None:
+          p2, = plt.plot(t[islope2], xfilt[islope2])
+          p3, = plt.plot(zc2, np.zeros(len(zc2)), '*g', markersize=14)
+          p4, = plt.plot(t[islope2], datafit2, '--g', linewidth=2) 
+          plt.text(Pick + 0.02, max(xraw) / 2, '%s' % FM, fontsize=14) 
+          ax = plt.gca()
+          ax.set_xlim([t[islope2[0][0]] - 0.1, t[islope2[0][len(islope2) - 1]] + 0.3])
+       plt.xlabel('Time [s] since %s' % Xraw[0].stats.starttime)
+       plt.yticks([])
+       plt.show()
+       raw_input()
+       plt.close(iplot)
+
+    return FM
+
+if __name__ == "__main__":
+   parser = argparse.ArgumentParser()
+   parser.add_argument('--Xraw', type=~obspy.core.stream.Stream, help='unfiltered time series (seismogram) read with obspy module read')
+   parser.add_argument('--Xfilt', type=~obspy.core.stream.Stream, help='filtered time series (seismogram) read with obspy module read')
+   parser.add_argument('--pickwin', type=float, help='length of pick window [s] for first motion determination')
+   parser.add_argument('--Pick', type=float, help='Onset time of most likely pick')
+   parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs')
+   args = parser.parse_args()
+   earllatepicker(args.Xraw, args.Xfilt, args.pickwin, args.Pick, args.iplot)
+

pylot/core/pick/run_makeCF.py

@@ -11,6 +11,8 @@ import matplotlib.pyplot as plt
 import numpy as np
 from pylot.core.pick.CharFuns import CharacteristicFunction
 from pylot.core.pick.Picker import AutoPicking
+from earllatepicker import earllatepicker
+from fmpicker import fmpicker
 import glob
 import argparse
 
@@ -28,9 +30,9 @@ def run_makeCF(project, database, event, iplot, station=None):
     addnoise = 0.001         #add noise to seismogram for stable AR prediction
     arzorder = 2             #chosen order of AR process, vertical component
     arhorder = 4             #chosen order of AR process, horizontal components
-    TSNRhos = [5, 0.5, 1, 0.1] #window lengths [s] for calculating SNR for earliest/latest pick and quality assessment
+    TSNRhos = [5, 0.5, 1, .6] #window lengths [s] for calculating SNR for earliest/latest pick and quality assessment
                              #from HOS-CF [noise window, safety gap, signal window, slope determination window]
-    TSNRarz = [5, 0.5, 1, 0.5] #window lengths [s] for calculating SNR for earliest/lates pick and quality assessment
+    TSNRarz = [5, 0.5, 1, 1.0] #window lengths [s] for calculating SNR for earliest/lates pick and quality assessment
                              #from ARZ-CF
     #get waveform data
     if station:
@@ -70,17 +72,20 @@ def run_makeCF(project, database, event, iplot, station=None):
           aiccf = AICcf(st_copy, cuttimes) #instance of AICcf
           ##############################################################
           #get prelimenary onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking
-          aicpick = AICPicker(aiccf, None, TSNRhos, 3, 10, None, 0.1)
+          aicpick = AICPicker(aiccf, TSNRhos, 3, 10, None, 0.1)
           ##############################################################
           #get refined onset time from HOS-CF using class Picker
-          hospick = PragPicker(hoscf, None, TSNRhos, 2, 10, 0.001, 0.2, aicpick.getpick())
+          hospick = PragPicker(hoscf, TSNRhos, 2, 10, 0.001, 0.2, aicpick.getpick())
+          #############################################################
           #get earliest and latest possible picks
-          hosELpick = EarlLatePicker(hoscf, 1.5, TSNRhos, None, 10, None, None, hospick.getpick())
+          st_copy[0].data = tr_filt.data
+          [lpickhos, epickhos, pickerrhos] = earllatepicker(st_copy, 1.5, TSNRhos, hospick.getpick(), 10)
+          #############################################################
+          #get first motion of onset
+          hosfm = fmpicker(st, st_copy, 0.2, hospick.getpick(), 11)
           ##############################################################
           #calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction
-          #get stream object of filtered data
+          arzcf = ARZcf(st, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf
-          st_copy[0].data = tr_filt.data
-          arzcf = ARZcf(st_copy, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf
           ##############################################################
           #calculate AIC-ARZ-CF using subclass AICcf of class CharacteristicFunction
           #class needs stream object => build it
@@ -90,12 +95,13 @@ def run_makeCF(project, database, event, iplot, station=None):
           araiccf = AICcf(st_copy, cuttimes, tpredz, 0, tdetz) #instance of AICcf
           ##############################################################
           #get onset time from AIC-ARZ-CF using subclass AICPicker of class AutoPicking
-          aicarzpick = AICPicker(araiccf, 1.5, TSNRarz, 2, 10, None, 0.1)
+          aicarzpick = AICPicker(araiccf, TSNRarz, 2, 10, None, 0.1)
           ##############################################################
           #get refined onset time from ARZ-CF using class Picker
-          arzpick = PragPicker(arzcf, 1.5, TSNRarz, 2.0, 10, 0.1, 0.05, aicarzpick.getpick())
+          arzpick = PragPicker(arzcf, TSNRarz, 2.0, 10, 0.1, 0.05, aicarzpick.getpick())
           #get earliest and latest possible picks
-          arzELpick = EarlLatePicker(arzcf, 1.5, TSNRarz, None, 10, None, None, arzpick.getpick())
+          st_copy[0].data = tr_filt.data
+          [lpickarz, epickarz, pickerrarz] = earllatepicker(st_copy, 1.5, TSNRarz, arzpick.getpick(), 10)
     elif not wfzfiles:
        print 'No vertical component data found!'
 
@@ -131,12 +137,23 @@ def run_makeCF(project, database, event, iplot, station=None):
           arhaiccf = AICcf(H_copy, cuttimes, tpredh, 0, tdeth) #instance of AICcf
           ##############################################################
           #get onset time from AIC-ARH-CF using subclass AICPicker of class AutoPicking
-          aicarhpick = AICPicker(arhaiccf, 1.5, TSNRarz, 4, 10, None, 0.1)
+          aicarhpick = AICPicker(arhaiccf, TSNRarz, 4, 10, None, 0.1)
           ###############################################################
           #get refined onset time from ARH-CF using class Picker
-          arhpick = PragPicker(arhcf, 1.5, TSNRarz, 2.5, 10, 0.1, 0.05, aicarhpick.getpick())
+          arhpick = PragPicker(arhcf, TSNRarz, 2.5, 10, 0.1, 0.05, aicarhpick.getpick())
           #get earliest and latest possible picks
-          arhELpick = EarlLatePicker(arhcf, 1.5, TSNRarz, None, 10, None, None, arhpick.getpick())
+          H_copy[0].data = trH1_filt.data
+          [lpickarh1, epickarh1, pickerrarh1] = earllatepicker(H_copy, 1.5, TSNRarz, arhpick.getpick(), 10)
+          H_copy[0].data = trH2_filt.data
+          [lpickarh2, epickarh2, pickerrarh2] = earllatepicker(H_copy, 1.5, TSNRarz, arhpick.getpick(), 10)
+          #get earliest pick of both earliest possible picks
+          epick = [epickarh1, epickarh2]
+          lpick = [lpickarh1, lpickarh2]
+          pickerr = [pickerrarh1, pickerrarh2]
+          ipick =np.argmin([epickarh1, epickarh2])
+          epickarh = epick[ipick]
+          lpickarh = lpick[ipick]
+          pickerrarh = pickerr[ipick]
 
           #create stream with 3 traces
           #merge streams
@@ -158,8 +175,6 @@ def run_makeCF(project, database, event, iplot, station=None):
           AllC[2].data = All3_filt.data
           #calculate AR3C-CF using subclass AR3Ccf of class CharacteristicFunction
           ar3ccf = AR3Ccf(AllC, cuttimes, tpredz, arhorder, tdetz, addnoise) #instance of AR3Ccf
-          #get earliest and latest possible pick from initial ARH-pick
-          ar3cELpick = EarlLatePicker(ar3ccf, 1.5, TSNRarz, None, 10, None, None, arhpick.getpick())
           ##############################################################
           if iplot:
              #plot vertical trace
@@ -177,16 +192,16 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([hospick.getpick(), hospick.getpick()], [-1.3, 1.3], 'r', linewidth=2)
              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([hosELpick.getLpick(), hosELpick.getLpick()], [-1.1, 1.1], 'r--')
+             plt.plot([lpickhos, lpickhos], [-1.1, 1.1], 'r--')
-             plt.plot([hosELpick.getEpick(), hosELpick.getEpick()], [-1.1, 1.1], 'r--')
+             plt.plot([epickhos, epickhos], [-1.1, 1.1], 'r--')
              plt.plot([aicarzpick.getpick(), aicarzpick.getpick()], [-1.2, 1.2], 'y', linewidth=2)
              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', linewidth=2)
              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.plot([arzELpick.getLpick(), arzELpick.getLpick()], [-1.2, 1.2], 'g--')
+             plt.plot([lpickarz, lpickarz], [-1.2, 1.2], 'g--')
-             plt.plot([arzELpick.getEpick(), arzELpick.getEpick()], [-1.2, 1.2], 'g--')
+             plt.plot([epickarz, epickarz], [-1.2, 1.2], 'g--')
              plt.yticks([])
              plt.ylim([-1.5, 1.5])
              plt.xlabel('Time [s]')
@@ -211,12 +226,10 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([arhpick.getpick(), arhpick.getpick()], [-1, 1], 'r')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [1, 1], 'r')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [-1, -1], 'r')
-             plt.plot([arhELpick.getLpick(), arhELpick.getLpick()], [-0.8, 0.8], 'r--')
+             plt.plot([lpickarh, lpickarh], [-0.8, 0.8], 'r--')
-             plt.plot([arhELpick.getEpick(), arhELpick.getEpick()], [-0.8, 0.8], 'r--')
+             plt.plot([epickarh, epickarh], [-0.8, 0.8], 'r--')
-             plt.plot([arhpick.getpick() + arhELpick.getPickError(), arhpick.getpick() + arhELpick.getPickError()], \
+             plt.plot([arhpick.getpick() + pickerrarh, arhpick.getpick() + pickerrarh], [-0.2, 0.2], 'r--')
-                       [-0.2, 0.2], 'r--')
+             plt.plot([arhpick.getpick() - pickerrarh, arhpick.getpick() - pickerrarh], [-0.2, 0.2], 'r--')
-             plt.plot([arhpick.getpick() - arhELpick.getPickError(), arhpick.getpick() - arhELpick.getPickError()], \
-                       [-0.2, 0.2], 'r--')
              plt.yticks([])
              plt.ylim([-1.5, 1.5])
              plt.ylabel('Normalized Counts')
@@ -233,12 +246,10 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([arhpick.getpick(), arhpick.getpick()], [-1, 1], 'r')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [1, 1], 'r')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [-1, -1], 'r')
-             plt.plot([arhELpick.getLpick(), arhELpick.getLpick()], [-0.8, 0.8], 'r--')
+             plt.plot([lpickarh, lpickarh], [-0.8, 0.8], 'r--')
-             plt.plot([arhELpick.getEpick(), arhELpick.getEpick()], [-0.8, 0.8], 'r--')
+             plt.plot([epickarh, epickarh], [-0.8, 0.8], 'r--')
-             plt.plot([arhpick.getpick() + arhELpick.getPickError(), arhpick.getpick() + arhELpick.getPickError()], \
+             plt.plot([arhpick.getpick() + pickerrarh, arhpick.getpick() + pickerrarh], [-0.2, 0.2], 'r--')
-                       [-0.2, 0.2], 'r--')
+             plt.plot([arhpick.getpick() - pickerrarh, arhpick.getpick() - pickerrarh], [-0.2, 0.2], 'r--')
-             plt.plot([arhpick.getpick() - arhELpick.getPickError(), arhpick.getpick() - arhELpick.getPickError()], \
-                       [-0.2, 0.2], 'r--')
              plt.title([trH2_filt.stats.station, trH2_filt.stats.channel])
              plt.yticks([])
              plt.ylim([-1.5, 1.5])
@@ -252,8 +263,6 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([arhpick.getpick(), arhpick.getpick()], [-1, 1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [-1, -1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [1, 1], 'b')
-             plt.plot([ar3cELpick.getLpick(), ar3cELpick.getLpick()], [-0.8, 0.8], 'b--')
-             plt.plot([ar3cELpick.getEpick(), ar3cELpick.getEpick()], [-0.8, 0.8], 'b--')
              plt.yticks([])
              plt.xticks([])
              plt.ylabel('Normalized Counts')
@@ -266,8 +275,6 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([arhpick.getpick(), arhpick.getpick()], [-1, 1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [-1, -1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [1, 1], 'b')
-             plt.plot([ar3cELpick.getLpick(), ar3cELpick.getLpick()], [-0.8, 0.8], 'b--')
-             plt.plot([ar3cELpick.getEpick(), ar3cELpick.getEpick()], [-0.8, 0.8], 'b--')
              plt.yticks([])
              plt.xticks([])
              plt.ylabel('Normalized Counts')
@@ -278,8 +285,6 @@ def run_makeCF(project, database, event, iplot, station=None):
              plt.plot([arhpick.getpick(), arhpick.getpick()], [-1, 1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [-1, -1], 'b')
              plt.plot([arhpick.getpick()-0.5, arhpick.getpick()+0.5], [1, 1], 'b')
-             plt.plot([ar3cELpick.getLpick(), ar3cELpick.getLpick()], [-0.8, 0.8], 'b--')
-             plt.plot([ar3cELpick.getEpick(), ar3cELpick.getEpick()], [-0.8, 0.8], 'b--')
              plt.yticks([])
              plt.ylabel('Normalized Counts')
              plt.title([trH2_filt.stats.station, trH2_filt.stats.channel])