From 59930c32381ee254978a048d9e5f8fd5db1627dd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Mon, 15 Dec 2014 15:03:41 +0100 Subject: [PATCH 1/5] Implemented pragmatic picking algorithm developed by TM, JL, and LK --- pylot/core/pick/Picker.py | 116 +++++++++++++++++++++++++++++++++----- 1 file changed, 101 insertions(+), 15 deletions(-) diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index e5857e49..b7cb7465 100644 --- a/pylot/core/pick/Picker.py +++ b/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(), - Tsmooth=self.getTsmooth()) + aus=self.getaus(), + 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): - return self.peps + def getaus(self): + return self.aus - def setpeps(self, peps): - self.peps = peps + def setaus(self, aus): + 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 + From fa58ec2aee26bab1fa8370f4133fe3a7cddb6a66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Mon, 15 Dec 2014 15:04:48 +0100 Subject: [PATCH 2/5] Modified for applying pragmatic picking algorithm, new class PragPicker in Picker.py --- pylot/core/pick/run_makeCF.py | 17 +++++++++++++---- 1 file changed, 13 insertions(+), 4 deletions(-) diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index 182404d7..efbd84d0 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/pylot/core/pick/run_makeCF.py @@ -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') From 13b8a9daec738090e9d8b867b38509e7378f9299 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 16 Dec 2014 16:13:52 +0100 Subject: [PATCH 3/5] Debugged --- pylot/core/pick/run_makeCF.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index efbd84d0..37bfeee6 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/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 @@ -182,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') @@ -194,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') From 2fcf325a6e01119545da7df17c88a1ef2b5df1ab Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 16 Dec 2014 16:15:53 +0100 Subject: [PATCH 4/5] Debugged getDataArray, same data lengths are now guaranteed --- pylot/core/pick/CharFuns.py | 76 ++++++++++++++++++++++++++++++++----- 1 file changed, 67 insertions(+), 9 deletions(-) diff --git a/pylot/core/pick/CharFuns.py b/pylot/core/pick/CharFuns.py index fed1bc29..64c4f06e 100644 --- a/pylot/core/pick/CharFuns.py +++ b/pylot/core/pick/CharFuns.py @@ -117,8 +117,8 @@ class CharacteristicFunction(object): def getTimeArray(self): if self.getTime1(): - incr = self.getARdetStep() - self.TimeArray = np.arange(0, len(self.getCF()) * incr[0], incr[0]) + self.getCut()[0] \ + incr = self.getARdetStep()[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: - start = 0 - else: - start = self.cut[0] / self.dt - stop = self.cut[1] / self.dt 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 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): From f0d60de7454a4c2137bed6f88169b2d3aaf39943 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 17 Dec 2014 12:16:32 +0100 Subject: [PATCH 5/5] add save data method --- QtPyLoT.py | 13 ++++++++++--- pylot/core/read/data.py | 4 +++- pylot/core/util/errors.py | 3 +++ 3 files changed, 16 insertions(+), 4 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 22579058..9176326a 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -143,18 +143,25 @@ class MainWindow(QMainWindow): filt = """Supported event formats (*.mat *.qml *.xml *.kor *.evt)""" caption = 'Select event to open' - fname = QFileDialog().getOpenFileName(self, caption=caption, + self.fname = QFileDialog().getOpenFileName(self, caption=caption, filter=filt) else: - fname = unicode(action.data().toString()) + self.fname = unicode(action.data().toString()) if not self.okToContinue(): return else: return if fname: - self.data = Data(evtdata=fname) + self.fname = fname + self.data = Data(evtdata=self.fname) def saveData(self): + settings = QSettings() + exform = settings.value('data/exportFormat', 'None') + try: + self.data.exportEvent(self.fname, exform) + except FormatError: + return False return True def getComponent(self): diff --git a/pylot/core/read/data.py b/pylot/core/read/data.py index 82c75db8..3be72605 100644 --- a/pylot/core/read/data.py +++ b/pylot/core/read/data.py @@ -7,6 +7,7 @@ from obspy.core import (read, Stream) from obspy import readEvents from obspy.core.event import (Event, Catalog) from pylot.core.util import fnConstructor +from pylot.core.util.errors import FormatError class Data(object): @@ -71,7 +72,8 @@ class Data(object): from pylot.core.util.defaults import OUTPUTFORMATS if evtformat.strip() not in OUTPUTFORMATS.values(): - evtformat = OUTPUTFORMATS.values()[0] + errmsg = 'selected format {0} not available'.format(evtformat) + raise FormatError(errmsg) if fnout is None: ID = self.evtdata.getEventID() diff --git a/pylot/core/util/errors.py b/pylot/core/util/errors.py index 68125c2d..0b2a4780 100644 --- a/pylot/core/util/errors.py +++ b/pylot/core/util/errors.py @@ -8,3 +8,6 @@ Created on Thu Mar 20 09:47:04 2014 class OptionsError(Exception): pass + +class FormatError(Exception): + pass