From 8ba34db05c6fbdf4e4d1573464b5040020cf55b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 18 Mar 2015 14:44:08 +0100 Subject: [PATCH 1/5] New function to calculate earliest and latest possible pick from a given initial (most likely) pick. Replaces class EarlLatePicker of object Picker. --- pylot/core/pick/earllatepicker.py | 138 ++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100755 pylot/core/pick/earllatepicker.py diff --git a/pylot/core/pick/earllatepicker.py b/pylot/core/pick/earllatepicker.py new file mode 100755 index 00000000..1a2663f7 --- /dev/null +++ b/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) From 16ae4bdfe952f163fd5bc3cdf2d9b3c2c5e0fa02 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 18 Mar 2015 14:45:08 +0100 Subject: [PATCH 2/5] Modified for using new function earllatepicker instead of removed class EarlLatePicker of object Picker. --- pylot/core/pick/run_makeCF.py | 75 +++++++++++++++++------------------ 1 file changed, 37 insertions(+), 38 deletions(-) diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index 5b801715..de081e23 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/pylot/core/pick/run_makeCF.py @@ -11,6 +11,7 @@ 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 import glob import argparse @@ -28,9 +29,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 +71,16 @@ 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()) + [lpickhos, epickhos, pickerrhos] = earllatepicker(st, 1.5, TSNRhos, hospick.getpick(), 10) ############################################################## #calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction - #get stream object of filtered data - st_copy[0].data = tr_filt.data - arzcf = ARZcf(st_copy, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf + arzcf = ARZcf(st, 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 +90,12 @@ 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()) + [lpickarz, epickarz, pickerrarz] = earllatepicker(st, 1.5, TSNRarz, arzpick.getpick(), 10) elif not wfzfiles: print 'No vertical component data found!' @@ -131,12 +131,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 +169,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 +186,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([hosELpick.getEpick(), hosELpick.getEpick()], [-1.1, 1.1], 'r--') + plt.plot([lpickhos, lpickhos], [-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([arzELpick.getEpick(), arzELpick.getEpick()], [-1.2, 1.2], 'g--') + plt.plot([lpickarz, lpickarz], [-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 +220,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([arhELpick.getEpick(), arhELpick.getEpick()], [-0.8, 0.8], 'r--') - plt.plot([arhpick.getpick() + arhELpick.getPickError(), arhpick.getpick() + arhELpick.getPickError()], \ - [-0.2, 0.2], 'r--') - plt.plot([arhpick.getpick() - arhELpick.getPickError(), arhpick.getpick() - arhELpick.getPickError()], \ - [-0.2, 0.2], 'r--') + plt.plot([lpickarh, lpickarh], [-0.8, 0.8], 'r--') + plt.plot([epickarh, epickarh], [-0.8, 0.8], 'r--') + plt.plot([arhpick.getpick() + pickerrarh, arhpick.getpick() + pickerrarh], [-0.2, 0.2], 'r--') + plt.plot([arhpick.getpick() - pickerrarh, arhpick.getpick() - pickerrarh], [-0.2, 0.2], 'r--') plt.yticks([]) plt.ylim([-1.5, 1.5]) plt.ylabel('Normalized Counts') @@ -233,12 +240,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([arhELpick.getEpick(), arhELpick.getEpick()], [-0.8, 0.8], 'r--') - plt.plot([arhpick.getpick() + arhELpick.getPickError(), arhpick.getpick() + arhELpick.getPickError()], \ - [-0.2, 0.2], 'r--') - plt.plot([arhpick.getpick() - arhELpick.getPickError(), arhpick.getpick() - arhELpick.getPickError()], \ - [-0.2, 0.2], 'r--') + plt.plot([lpickarh, lpickarh], [-0.8, 0.8], 'r--') + plt.plot([epickarh, epickarh], [-0.8, 0.8], 'r--') + plt.plot([arhpick.getpick() + pickerrarh, arhpick.getpick() + pickerrarh], [-0.2, 0.2], 'r--') + plt.plot([arhpick.getpick() - pickerrarh, arhpick.getpick() - pickerrarh], [-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 +257,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 +269,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 +279,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]) From 787cac7d686214f66a09eb293a5f69cb68f212bb Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 18 Mar 2015 14:45:49 +0100 Subject: [PATCH 3/5] Removed class EarlLatePicker, replaced by new function earllatepicker. --- pylot/core/pick/Picker.py | 339 +------------------------------------- 1 file changed, 3 insertions(+), 336 deletions(-) diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index dc928a85..b3579663 100644 --- a/pylot/core/pick/Picker.py +++ b/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 - From a606b030e2f620b557f92975f3a2e0764b7aaf65 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Thu, 19 Mar 2015 14:32:50 +0100 Subject: [PATCH 4/5] New function to derive automatically first motion (polarity) of phase onset based on zero crossings and slope determination. --- pylot/core/pick/fmpicker.py | 183 ++++++++++++++++++++++++++++++++++++ 1 file changed, 183 insertions(+) create mode 100755 pylot/core/pick/fmpicker.py diff --git a/pylot/core/pick/fmpicker.py b/pylot/core/pick/fmpicker.py new file mode 100755 index 00000000..6aee8eb3 --- /dev/null +++ b/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) + From dc78abed09bd5f05322dbfe1c0c967e9c745f70e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Thu, 19 Mar 2015 14:36:56 +0100 Subject: [PATCH 5/5] Modified to handle new function fmpicker. --- pylot/core/pick/run_makeCF.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index de081e23..690fdcff 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/pylot/core/pick/run_makeCF.py @@ -12,6 +12,7 @@ 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 @@ -77,7 +78,11 @@ def run_makeCF(project, database, event, iplot, station=None): hospick = PragPicker(hoscf, TSNRhos, 2, 10, 0.001, 0.2, aicpick.getpick()) ############################################################# #get earliest and latest possible picks - [lpickhos, epickhos, pickerrhos] = earllatepicker(st, 1.5, TSNRhos, hospick.getpick(), 10) + 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 arzcf = ARZcf(st, cuttimes, tpredz, arzorder, tdetz, addnoise) #instance of ARZcf @@ -95,7 +100,8 @@ def run_makeCF(project, database, event, iplot, station=None): #get refined onset time from ARZ-CF using class Picker arzpick = PragPicker(arzcf, TSNRarz, 2.0, 10, 0.1, 0.05, aicarzpick.getpick()) #get earliest and latest possible picks - [lpickarz, epickarz, pickerrarz] = earllatepicker(st, 1.5, TSNRarz, arzpick.getpick(), 10) + 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!'