From 2db6fde7090fd987fc7b9fd6e89d068f0460631e Mon Sep 17 00:00:00 2001
From: Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de>
Date: Sat, 26 Aug 2017 02:12:31 +0200
Subject: [PATCH] [add] check second maximum of HOS/AR-CF
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit
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 | 28 ++++++++++++++++++++++++++--
pylot/core/pick/autopick.py | 20 +++++++++++++-------
pylot/core/pick/picker.py | 14 ++++++++++----
pylot/core/pick/utils.py | 24 ++++++++++++++++++++++++
4 files changed, 73 insertions(+), 13 deletions(-)
diff --git a/pylot/core/io/default_parameters.py b/pylot/core/io/default_parameters.py
index 945e0497..f0150112 100644
--- a/pylot/core/io/default_parameters.py
+++ b/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',
diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py
index 0b60f4f7..4d39075d 100644
--- a/pylot/core/pick/autopick.py
+++ b/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,
- aicpick.getpick(), fig)
+ refPpick = PragPicker(cf2, tsnrz, pickwinP, iplot=iplot, aus=ausP, Tsmooth=tsmoothP,
+ 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,
- aictsmoothS, fig=fig)
+ aicarhpick = AICPicker(haiccf, tsnrh, pickwinS, checkwindow=checkwindowS,
+ 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,
- tsmoothS, aicarhpick.getpick(), fig)
+ refSpick = PragPicker(arhcf2, tsnrh, pickwinS, iplot=iplot, aus=ausS,
+ Tsmooth=tsmoothS, Pick1=aicarhpick.getpick(), fig=fig)
mpickS = refSpick.getpick()
#############################################################
if mpickS is not None:
diff --git a/pylot/core/pick/picker.py b/pylot/core/pick/picker.py
index 1715746c..59bbf143 100644
--- a/pylot/core/pick/picker.py
+++ b/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
- icfmax = np.argmax(self.Data[0].data)
+ # minimum in AIC function
+ 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))
diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py
index 0dce76fd..33c3bbd5 100644
--- a/pylot/core/pick/utils.py
+++ b/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