[add] check second maximum of HOS/AR-CF

Use the second maximum if it is larger than first maximum * minfactor. This helps when two onsets are close together (eg. S and SKS at 90°) and the second onset is stronger.

pylot/core/io/default_parameters.py

@@ -273,6 +273,26 @@ defaults = {'rootpath': {'type': str,
                       'value': 1.5,
                       'namestring': 'Noise factor S'},
 
+            'checkwindowP': {'type': float,
+                      'tooltip': 'time window before HOS/AR-maximum to check for smaller maxima [s]',
+                      'value': 10.0,
+                      'namestring': 'Check Window P'},
+
+            'minfactorP': {'type': float,
+                      'tooltip': 'Second maximum must be at least minfactor * first maximum [-]',
+                      'value': 0.7,
+                      'namestring': 'Minimum Factor P'},
+
+            'checkwindowS': {'type': float,
+                             'tooltip': 'time window before AR-maximum to check for smaller maxima [s]',
+                             'value': 10.0,
+                             'namestring': 'Check Window S'},
+
+            'minfactorS': {'type': float,
+                           'tooltip': 'Second maximum must be at least minfactor * first maximum [-]',
+                           'value': 0.7,
+                           'namestring': 'Minimum Factor S'},
+
             'minfmweight': {'type': int,
                             'tooltip': 'minimum required P weight for first-motion determination',
                             'value': 1,
@@ -450,7 +470,9 @@ settings_special_pick = {
         'aictsmooth',
         'tsmoothP',
         'ausP',
-        'nfacP'],
+        'nfacP',
+        'checkwindowP',
+        'minfactorP'],
     'h': [
         'algoS',
         'tdet1h',
@@ -464,7 +486,9 @@ settings_special_pick = {
         'aictsmoothS',
         'tsmoothS',
         'ausS',
-        'nfacS'],
+        'nfacS',
+        'checkwindowS',
+        'minfactorS'],
     'fm': [
         'minfmweight',
         'minFMSNR',

pylot/core/pick/autopick.py

@@ -175,6 +175,10 @@ def autopickstation(wfstream, pickparam, verbose=False,
     # parameter to check for spuriously picked S onset
     zfac = pickparam.get('zfac')
     # path to inventory-, dataless- or resp-files
+    checkwindowP = pickparam.get('checkwindowP')
+    minfactorP = pickparam.get('minfactorP')
+    checkwindowS = pickparam.get('checkwindowS')
+    minfactorS = pickparam.get('minfactorS')
 
     # initialize output
     Pweight = 4  # weight for P onset
@@ -322,7 +326,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
             fig = fig_dict[key]
         else:
             fig = None
-        aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP, fig=fig)
+        aicpick = AICPicker(aiccf, tsnrz, pickwinP, checkwindow=checkwindowP, minfactor=minfactorP,
+                            iplot=iplot, Tsmooth=tsmoothP, fig=fig)
         # add pstart and pstop to aic plot
         if fig:
             for ax in fig.axes:
@@ -442,8 +447,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
                 fig = fig_dict['refPpick']
             else:
                 fig = None
-            refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot, ausP, tsmoothP,
+            refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot=iplot, aus=ausP, Tsmooth=tsmoothP,
-                                  aicpick.getpick(), fig)
+                                  Pick1 = aicpick.getpick(), fig=fig)
             mpickP = refPpick.getpick()
             #############################################################
             if mpickP is not None:
@@ -624,8 +629,9 @@ def autopickstation(wfstream, pickparam, verbose=False,
             fig = fig_dict['aicARHfig']
         else:
             fig = None
-        aicarhpick = AICPicker(haiccf, tsnrh, pickwinS, iplot, None,
+        aicarhpick = AICPicker(haiccf, tsnrh, pickwinS, checkwindow=checkwindowS,
-                               aictsmoothS, fig=fig)
+                               minfactor=minfactorS, iplot=iplot, Tsmooth=aictsmoothS,
+                               fig=fig)
         ###############################################################
         # go on with processing if AIC onset passes quality control
         slope = aicarhpick.getSlope()
@@ -686,8 +692,8 @@ def autopickstation(wfstream, pickparam, verbose=False,
                 fig = fig_dict['refSpick']
             else:
                 fig = None
-            refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot, ausS,
+            refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot=iplot, aus=ausS,
-                                  tsmoothS, aicarhpick.getpick(), fig)
+                                  Tsmooth=tsmoothS, Pick1=aicarhpick.getpick(), fig=fig)
             mpickS = refSpick.getpick()
             #############################################################
             if mpickS is not None:

pylot/core/pick/picker.py

@@ -24,7 +24,7 @@ import warnings
 import matplotlib.pyplot as plt
 import numpy as np
 from pylot.core.pick.charfuns import CharacteristicFunction
-from pylot.core.pick.utils import getnoisewin, getsignalwin
+from pylot.core.pick.utils import getnoisewin, getsignalwin, get_maximum_index
 
 
 class AutoPicker(object):
@@ -35,7 +35,7 @@ class AutoPicker(object):
 
     warnings.simplefilter('ignore')
 
-    def __init__(self, cf, TSNR, PickWindow, iplot=0, aus=None, Tsmooth=None, Pick1=None, fig=None):
+    def __init__(self, cf, TSNR, PickWindow, checkwindow=None, minfactor=None, iplot=0, aus=None, Tsmooth=None, Pick1=None, fig=None):
         '''
         :param: cf, characteristic function, on which the picking algorithm is applied
         :type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
@@ -74,6 +74,8 @@ class AutoPicker(object):
         self.setTsmooth(Tsmooth)
         self.setpick1(Pick1)
         self.fig = fig
+        self.setCheckWindow(checkwindow)
+        self.minfactor = minfactor
         self.calcPick()
 
     def __str__(self):
@@ -89,6 +91,10 @@ class AutoPicker(object):
                    aus=self.getaus(),
                    Tsmooth=self.getTsmooth(),
                    Pick1=self.getpick1())
+    def setCheckWindow(self, checkwindow):
+        '''convert checkwindow to samples'''
+        if checkwindow:
+            self.checkwindow = int(checkwindow / self.Data[0].stats.delta)
 
     def getTSNR(self):
         return self.TSNR
@@ -187,8 +193,8 @@ class AICPicker(AutoPicker):
         offset = abs(min(aic) - min(aicsmooth))
         aicsmooth = aicsmooth - offset
         # get maximum of HOS/AR-CF as startimg point for searching
-        # minimum in AIC function 
+        # minimum in AIC function
-        icfmax = np.argmax(self.Data[0].data)
+        icfmax = get_maximum_index(self.Data[0].data, self.checkwindow, self.minfactor)
 
         # find minimum in AIC-CF front of maximum of HOS/AR-CF
         lpickwindow = int(round(self.PickWindow / self.dt))

pylot/core/pick/utils.py

@@ -13,6 +13,7 @@ import warnings
 import matplotlib.pyplot as plt
 import numpy as np
 from obspy.core import Stream, UTCDateTime
+from scipy.signal import argrelextrema
 
 
 def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None):
@@ -1171,6 +1172,29 @@ def removePicksAbove(pickDic, minWeight):
                     newdic[eventKey][station].update({phasename: phaseinfo})
     return newdic
 
+
+def get_maximum_index(data, checkwindow, minfactor):
+    '''get maximum of CF as starting point, then check for highest local maximum
+    in front of it.
+    return second maximum if its larger than first maximum * minfactor, else
+    return first maximum'''
+    icfmax1 = np.argmax(data)
+    imax_local = argrelextrema(data[icfmax1 - checkwindow:icfmax1], np.greater)  # indices of local maxima
+    if (len(imax_local[0]) > 0):
+        imax_local = imax_local[0] + icfmax1 - checkwindow
+        local_maxima = (imax_local, data[imax_local])
+        icfmax2 = local_maxima[0][np.where(local_maxima[1] == max(local_maxima[1]))][0]
+        if data[icfmax2] > data[icfmax1] * minfactor:
+            print("Found valid local maximum in front of first maximum")
+            return icfmax2
+        else:
+            print("First maximum is the largest: {}>{}".format(data[icfmax1],
+                                                               data[icfmax2]))
+            return icfmax1
+    else:
+        print("No local maxima found in check window")
+        return icfmax1
+
 if __name__ == '__main__':
     import doctest