From 60b9f176f0608f5db1a2c43dff2aa77200a352a1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Thu, 3 Sep 2015 14:55:25 +0200 Subject: [PATCH 01/13] Cosmetics, changed print commands to keep compatibility to Python 3. --- pylot/core/pick/utils.py | 87 ++++++++++++++++++++-------------------- 1 file changed, 43 insertions(+), 44 deletions(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 161942cb..f891b6fe 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -7,7 +7,6 @@ :author: Ludger Kueperkoch / MAGS2 EP3 working group """ - import numpy as np import scipy as sc import matplotlib.pyplot as plt @@ -44,7 +43,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): LPick = None EPick = None PickError = None - print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...' + 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, @@ -60,8 +59,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): ilup, = np.where(x[isignal] > nlevel) ildown, = np.where(x[isignal] < -nlevel) if not ilup.size and not ildown.size: - print 'earllatepicker: Signal lower than noise level!' - print 'Skip this trace!' + print ("earllatepicker: Signal lower than noise level!") + print ("Skip this trace!") return LPick, EPick, PickError il = min(np.min(ilup) if ilup.size else float('inf'), np.min(ildown) if ildown.size else float('inf')) @@ -143,7 +142,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): FM = None if Pick is not None: - print 'fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...' + print ("fmpicker: Get first motion (polarity) of onset using unfiltered seismogram...") xraw = Xraw[0].data xfilt = Xfilt[0].data @@ -182,15 +181,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): else: li1 = index1[0] if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 0: - print 'fmpicker: Onset on unfiltered trace too emergent for first motion determination!' + print ("fmpicker: Onset on unfiltered trace too emergent for first motion determination!") P1 = None else: imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]])) if imax1 == 0: imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][index1[1]]])) if imax1 == 0: - print 'fmpicker: Zero crossings too close!' - print 'Skip first motion determination!' + print ("fmpicker: Zero crossings too close!") + print ("Skip first motion determination!") return FM islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1])) @@ -224,15 +223,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): else: li2 = index2[0] if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 0: - print 'fmpicker: Onset on filtered trace too emergent for first motion determination!' + print ("fmpicker: Onset on filtered trace too emergent for first motion determination!") P2 = None else: imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]])) if imax2 == 0: imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][index2[1]]])) if imax2 == 0: - print 'fmpicker: Zero crossings too close!' - print 'Skip first motion determination!' + print ("fmpicker: Zero crossings too close!") + print ("Skip first motion determination!") return FM islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2])) @@ -256,7 +255,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): elif P1[0] > 0 and P2[0] <= 0: FM = '+' - print 'fmpicker: Found polarity %s' % FM + print ("fmpicker: Found polarity %s" % FM) if iplot > 1: plt.figure(iplot) @@ -331,10 +330,10 @@ def getSNR(X, TSNR, t1): # get signal window isignal = getsignalwin(t, t1, TSNR[2]) if np.size(inoise) < 1: - print 'getSNR: Empty array inoise, check noise window!' + print ("getSNR: Empty array inoise, check noise window!") return elif np.size(isignal) < 1: - print 'getSNR: Empty array isignal, check signal window!' + print ("getSNR: Empty array isignal, check signal window!") return # demean over entire waveform @@ -372,7 +371,7 @@ def getnoisewin(t, t1, tnoise, tgap): inoise, = np.where((t <= max([t1 - tgap, 0])) \ & (t >= max([t1 - tnoise - tgap, 0]))) if np.size(inoise) < 1: - print 'getnoisewin: Empty array inoise, check noise window!' + print ("getnoisewin: Empty array inoise, check noise window!") return inoise @@ -396,7 +395,7 @@ def getsignalwin(t, t1, tsignal): isignal, = np.where((t <= min([t1 + tsignal, len(t)])) \ & (t >= t1)) if np.size(isignal) < 1: - print 'getsignalwin: Empty array isignal, check signal window!' + print ("getsignalwin: Empty array isignal, check signal window!") return isignal @@ -483,8 +482,8 @@ def wadaticheck(pickdic, dttolerance, iplot): # calculate vp/vs ratio before check vpvsr = p1[0] + 1 - print '###############################################' - print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr + print ("###############################################") + print ("wadaticheck: Average Vp/Vs ratio before check:", vpvsr) checkedPpicks = [] checkedSpicks = [] @@ -521,18 +520,18 @@ def wadaticheck(pickdic, dttolerance, iplot): # calculate vp/vs ratio after check cvpvsr = p2[0] + 1 - print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr - print 'wadatacheck: Skipped %d S pick(s).' % ibad + print ("wadaticheck: Average Vp/Vs ratio after check:", cvpvsr) + print ("wadatacheck: Skipped %d S pick(s)." % ibad) else: - print '###############################################' - print 'wadatacheck: Not enough checked S-P times available!' - print 'Skip Wadati check!' + print ("###############################################") + print ("wadatacheck: Not enough checked S-P times available!") + print ("Skip Wadati check!") checkedonsets = pickdic else: - print 'wadaticheck: Not enough S-P times available for reliable regression!' - print 'Skip wadati check!' + print ("wadaticheck: Not enough S-P times available for reliable regression!") + print ("Skip wadati check!") wfitflag = 1 # plot results @@ -592,7 +591,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): assert isinstance(X, Stream), "%s is not a stream object" % str(X) - print 'Checking signal length ...' + print ("Checking signal length ...") x = X[0].data t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, @@ -601,8 +600,8 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): # generate envelope function from Hilbert transform y = np.imag(sc.signal.hilbert(x)) e = np.sqrt(np.power(x, 2) + np.power(y, 2)) - # get noise window - inoise = getnoisewin(t, pick, TSNR[0], TSNR[1]) + # get noise window in front of pick plus saftey gap + inoise = getnoisewin(t, pick - 0.5, TSNR[0], TSNR[1]) # get signal window isignal = getsignalwin(t, pick, TSNR[2]) # calculate minimum adjusted signal level @@ -613,12 +612,12 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): numoverthr = len(np.where(e[isignal] >= minsiglevel)[0]) if numoverthr >= minnum: - print 'checksignallength: Signal reached required length.' + print ("checksignallength: Signal reached required length.") returnflag = 1 else: - print 'checksignallength: Signal shorter than required minimum signal length!' - print 'Presumably picked noise peak, pick is rejected!' - print '(min. signal length required:', minsiglength, 's)' + print ("checksignallength: Signal shorter than required minimum signal length!") + print ("Presumably picked noise peak, pick is rejected!") + print ("(min. signal length required:', minsiglength, 's)'") returnflag = 0 if iplot == 2: @@ -629,7 +628,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): p2, = plt.plot(t[inoise], e[inoise]) p3, = plt.plot(t[isignal],e[isignal], 'r') p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \ - [minsiglevel, minsiglevel], 'g') + [minsiglevel, minsiglevel], 'g', linewidth=2) p5, = plt.plot([pick, pick], [min(x), max(x)], 'b', linewidth=2) plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \ 'Envelope Signal Window', 'Minimum Signal Level', \ @@ -675,8 +674,8 @@ def checkPonsets(pickdic, dttolerance, iplot): stations.append(key) # apply jackknife bootstrapping on variance of P onsets - print '###############################################' - print 'checkPonsets: Apply jackknife bootstrapping on P-onset times ...' + print ("###############################################") + print ("checkPonsets: Apply jackknife bootstrapping on P-onset times ...") [xjack,PHI_pseudo,PHI_sub] = jackknife(Ppicks, 'VAR', 1) # get pseudo variances smaller than average variances # (times safety factor), these picks passed jackknife test @@ -684,7 +683,7 @@ def checkPonsets(pickdic, dttolerance, iplot): # these picks did not pass jackknife test badjk = np.where(PHI_pseudo > 2 * xjack) badjkstations = np.array(stations)[badjk] - print 'checkPonsets: %d pick(s) did not pass jackknife test!' % len(badjkstations) + print ("checkPonsets: %d pick(s) did not pass jackknife test!" % len(badjkstations)) # calculate median from these picks pmedian = np.median(np.array(Ppicks)[ij]) @@ -696,9 +695,9 @@ def checkPonsets(pickdic, dttolerance, iplot): goodstations = np.array(stations)[igood] badstations = np.array(stations)[ibad] - print 'checkPonsets: %d pick(s) deviate too much from median!' % len(ibad) - print 'checkPonsets: Skipped %d P pick(s) out of %d' % (len(badstations) \ - + len(badjkstations), len(stations)) + print ("checkPonsets: %d pick(s) deviate too much from median!" % len(ibad)) + print ("checkPonsets: Skipped %d P pick(s) out of %d" % (len(badstations) \ + + len(badjkstations), len(stations))) goodmarker = 'goodPonsetcheck' badmarker = 'badPonsetcheck' @@ -765,8 +764,8 @@ def jackknife(X, phi, h): g = len(X) / h if type(g) is not int: - print 'jackknife: Cannot divide quantity X in equal sized subgroups!' - print 'Choose another size for subgroups!' + print ("jackknife: Cannot divide quantity X in equal sized subgroups!") + print ("Choose another size for subgroups!") return PHI_jack, PHI_pseudo, PHI_sub else: # estimator of undisturbed spot check @@ -834,7 +833,7 @@ def checkZ4S(X, pick, zfac, checkwin, iplot): assert isinstance(X, Stream), "%s is not a stream object" % str(X) - print 'Check for spuriously picked S onset instead of P onset ...' + print ("Check for spuriously picked S onset instead of P onset ...") returnflag = 0 @@ -875,9 +874,9 @@ def checkZ4S(X, pick, zfac, checkwin, iplot): # vertical P-coda level must exceed horizontal P-coda level # zfac times encodalevel if zcodalevel < minsiglevel: - print 'checkZ4S: Maybe S onset? Skip this P pick!' + print ("checkZ4S: Maybe S onset? Skip this P pick!") else: - print 'checkZ4S: P onset passes checkZ4S test!' + print ("checkZ4S: P onset passes checkZ4S test!") returnflag = 1 if iplot > 1: From ab1d27747a7e5a7ce094dc2733958b1bef9ddd44 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Thu, 3 Sep 2015 15:42:20 +0200 Subject: [PATCH 02/13] Some reformating. --- pylot/core/pick/utils.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index f891b6fe..4476e6fc 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -7,6 +7,7 @@ :author: Ludger Kueperkoch / MAGS2 EP3 working group """ + import numpy as np import scipy as sc import matplotlib.pyplot as plt @@ -483,7 +484,7 @@ def wadaticheck(pickdic, dttolerance, iplot): # calculate vp/vs ratio before check vpvsr = p1[0] + 1 print ("###############################################") - print ("wadaticheck: Average Vp/Vs ratio before check:", vpvsr) + print ("wadaticheck: Average Vp/Vs ratio before check: %f" % vpvsr) checkedPpicks = [] checkedSpicks = [] @@ -520,8 +521,8 @@ def wadaticheck(pickdic, dttolerance, iplot): # calculate vp/vs ratio after check cvpvsr = p2[0] + 1 - print ("wadaticheck: Average Vp/Vs ratio after check:", cvpvsr) - print ("wadatacheck: Skipped %d S pick(s)." % ibad) + print ("wadaticheck: Average Vp/Vs ratio after check: %f" % cvpvsr) + print ("wadatacheck: Skipped %d S pick(s)" % ibad) else: print ("###############################################") print ("wadatacheck: Not enough checked S-P times available!") @@ -617,7 +618,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): else: print ("checksignallength: Signal shorter than required minimum signal length!") print ("Presumably picked noise peak, pick is rejected!") - print ("(min. signal length required:', minsiglength, 's)'") + print ("(min. signal length required: %s s)" % minsiglength) returnflag = 0 if iplot == 2: From 0dc10910784bd48378cac7aa5b612d107fa65846 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:12:25 +0200 Subject: [PATCH 03/13] restituteWFData: introduced return flag restflag to indicate whether restitution could be performed or not. --- pylot/core/read/data.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/pylot/core/read/data.py b/pylot/core/read/data.py index 30f82c36..ca50132d 100644 --- a/pylot/core/read/data.py +++ b/pylot/core/read/data.py @@ -228,6 +228,9 @@ class Data(object): :param streams: :return: """ + + restflag = 0 + if streams is None: st_raw = self.getWFData() st = st_raw.copy() @@ -236,7 +239,7 @@ class Data(object): for tr in st: # remove underscores - if tr.stats.station[3] == '_': + if len(tr.stats.station) > 3 and tr.stats.station[3] == '_': tr.stats.station = tr.stats.station[0:3] dlp = '%s/*.dless' % invdlpath invp = '%s/*.xml' % invdlpath @@ -273,6 +276,7 @@ class Data(object): 'date': st[ i].stats.starttime, 'units': "VEL"}) + restflag = 1 except ValueError as e: vmsg = '{0}'.format(e) print(vmsg) @@ -303,6 +307,7 @@ class Data(object): st[i].attach_response(inv) st[i].remove_response(output='VEL', pre_filt=prefilt) + restflag = 1 except ValueError as e: vmsg = '{0}'.format(e) print(vmsg) @@ -334,6 +339,7 @@ class Data(object): 'units': "VEL"} st[i].simulate(paz_remove=None, pre_filt=prefilt, seedresp=seedresp) + restflag = 1 except ValueError as e: vmsg = '{0}'.format(e) print(vmsg) @@ -346,7 +352,7 @@ class Data(object): print("Go on processing data without source parameter " "determination!") - return st + return st, restflag def getEvtData(self): """ From de608798b9b0aea80a8d8bc5f2bdf319390bdf98 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:13:52 +0200 Subject: [PATCH 04/13] Modified checking of signal length, uses RMS trace of all components now (if available). --- pylot/core/pick/autopick.py | 82 ++++++++++++++++++++++--------------- 1 file changed, 48 insertions(+), 34 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 6011f645..456023f5 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -41,9 +41,9 @@ def autopickevent(data, param): # quality control # median check and jackknife on P-onset times - jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, 2) + jk_checked_onsets = checkPonsets(all_onsets, mdttolerance, iplot) # check S-P times (Wadati) - return wadaticheck(jk_checked_onsets, wdttolerance, 2) + return wadaticheck(jk_checked_onsets, wdttolerance, iplot) def autopickstation(wfstream, pickparam): """ @@ -196,10 +196,36 @@ def autopickstation(wfstream, pickparam): ############################################################## if aicpick.getpick() is not None: # check signal length to detect spuriously picked noise peaks + # use all available components to avoid skipping correct picks + # on vertical traces with weak P coda z_copy[0].data = tr_filt.data - Pflag = checksignallength(z_copy, aicpick.getpick(), tsnrz, - minsiglength, \ - nfacsl, minpercent, iplot) + zne = z_copy + if len(ndat) == 0 or len(edat) == 0: + print ("One or more horizontal components missing!") + print ("Signal length only checked on vertical component!") + print ("Decreasing minsiglengh from %f to %f" \ + % (minsiglength, minsiglength / 2)) + Pflag = checksignallength(zne, aicpick.getpick(), tsnrz, + minsiglength / 2, \ + nfacsl, minpercent, iplot) + else: + # filter and taper horizontal traces + trH1_filt = edat.copy() + trH2_filt = ndat.copy() + trH1_filt.filter('bandpass', freqmin=bph1[0], + freqmax=bph1[1], \ + zerophase=False) + trH2_filt.filter('bandpass', freqmin=bph1[0], + freqmax=bph1[1], \ + zerophase=False) + trH1_filt.taper(max_percentage=0.05, type='hann') + trH2_filt.taper(max_percentage=0.05, type='hann') + zne += trH1_filt + zne += trH2_filt + Pflag = checksignallength(zne, aicpick.getpick(), tsnrz, + minsiglength, \ + nfacsl, minpercent, iplot) + if Pflag == 1: # check for spuriously picked S onset # both horizontal traces needed @@ -207,20 +233,6 @@ def autopickstation(wfstream, pickparam): print 'One or more horizontal components missing!' print 'Skipping control function checkZ4S.' else: - # filter and taper horizontal traces - trH1_filt = edat.copy() - trH2_filt = ndat.copy() - trH1_filt.filter('bandpass', freqmin=bph1[0], - freqmax=bph1[1], \ - zerophase=False) - trH2_filt.filter('bandpass', freqmin=bph1[0], - freqmax=bph1[1], \ - zerophase=False) - trH1_filt.taper(max_percentage=0.05, type='hann') - trH2_filt.taper(max_percentage=0.05, type='hann') - zne = z_copy - zne += trH1_filt - zne += trH2_filt Pflag = checkZ4S(zne, aicpick.getpick(), zfac, \ tsnrz[3], iplot) if Pflag == 0: @@ -515,18 +527,19 @@ def autopickstation(wfstream, pickparam): hdat = edat.copy() hdat += ndat h_copy = hdat.copy() - cordat = data.restituteWFData(invdir, h_copy) + [cordat, restflag] = data.restituteWFData(invdir, h_copy) # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude - if Sweight < 4: - wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) - else: - # use larger window for getting peak-to-peak amplitude - # as the S pick is quite unsure - wapp = WApp(cordat, mpickP, mpickP + sstop + \ - (0.5 * (mpickP + sstop)), iplot) + if restflag == 1: + if Sweight < 4: + wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) + else: + # use larger window for getting peak-to-peak amplitude + # as the S pick is quite unsure + wapp = WApp(cordat, mpickP, mpickP + sstop + \ + (0.5 * (mpickP + sstop)), iplot) - Ao = wapp.getwapp() + Ao = wapp.getwapp() else: print 'Bad initial (AIC) S-pick, skipping this onset!' @@ -544,12 +557,13 @@ def autopickstation(wfstream, pickparam): hdat = edat.copy() hdat += ndat h_copy = hdat.copy() - cordat = data.restituteWFData(invdir, h_copy) - # calculate WA-peak-to-peak amplitude - # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP \ - + sstop)), iplot) - Ao = wapp.getwapp() + [cordat, restflag] = data.restituteWFData(invdir, h_copy) + if restflag == 1: + # calculate WA-peak-to-peak amplitude + # using subclass WApp of superclass Magnitude + wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP \ + + sstop)), iplot) + Ao = wapp.getwapp() else: print 'autopickstation: No horizontal component data available or ' \ From 23430c9d90ff85932f4eebfa747f1a4180a71772 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:16:34 +0200 Subject: [PATCH 05/13] Modofied checksignallength: uses RMS trace of all components (if available) to check signal length. This avoids skipping of P pick, if P coda is very weak. If only vertical trace is available, rms of vertical trace is used instead with smaller required minimum signal length. --- pylot/core/pick/utils.py | 47 +++++++++++++++++++++++----------------- 1 file changed, 27 insertions(+), 20 deletions(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 4476e6fc..f94cac76 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -7,7 +7,7 @@ :author: Ludger Kueperkoch / MAGS2 EP3 working group """ - +import pdb import numpy as np import scipy as sc import matplotlib.pyplot as plt @@ -562,9 +562,9 @@ def wadaticheck(pickdic, dttolerance, iplot): def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): ''' Function to detect spuriously picked noise peaks. - Uses envelope to determine, how many samples [per cent] after - P onset are below certain threshold, calculated from noise - level times noise factor. + Uses RMS trace of all 3 components (if available) to determine, + how many samples [per cent] after P onset are below certain + threshold, calculated from noise level times noise factor. : param: X, time series (seismogram) : type: `~obspy.core.stream.Stream` @@ -594,23 +594,32 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): print ("Checking signal length ...") - x = X[0].data - t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, + if len(X) > 1: + # all three components available + # make sure, all components have equal lengths + ilen = min([len(X[0].data), len(X[1].data), len(X[2].data)]) + x1 = X[0][0:ilen] + x2 = X[1][0:ilen] + x3 = X[2][0:ilen] + # get RMS trace + rms = np.sqrt((np.power(x1, 2) + np.power(x2, 2) + np.power(x3, 2)) / 3) + else: + x1 = X[0].data + rms = np.sqrt(np.power(2, x1)) + + t = np.arange(0, ilen / X[0].stats.sampling_rate, X[0].stats.delta) - # generate envelope function from Hilbert transform - y = np.imag(sc.signal.hilbert(x)) - e = np.sqrt(np.power(x, 2) + np.power(y, 2)) # get noise window in front of pick plus saftey gap inoise = getnoisewin(t, pick - 0.5, TSNR[0], TSNR[1]) # get signal window - isignal = getsignalwin(t, pick, TSNR[2]) + isignal = getsignalwin(t, pick, minsiglength) # calculate minimum adjusted signal level - minsiglevel = max(e[inoise]) * nfac + minsiglevel = max(rms[inoise]) * nfac # minimum adjusted number of samples over minimum signal level minnum = len(isignal) * minpercent/100 # get number of samples above minimum adjusted signal level - numoverthr = len(np.where(e[isignal] >= minsiglevel)[0]) + numoverthr = len(np.where(rms[isignal] >= minsiglevel)[0]) if numoverthr >= minnum: print ("checksignallength: Signal reached required length.") @@ -623,16 +632,14 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): if iplot == 2: plt.figure(iplot) - p1, = plt.plot(t,x, 'k') - p2, = plt.plot(t[inoise], e[inoise], 'c') - p3, = plt.plot(t[isignal],e[isignal], 'r') - p2, = plt.plot(t[inoise], e[inoise]) - p3, = plt.plot(t[isignal],e[isignal], 'r') + p1, = plt.plot(t,rms, 'k') + p2, = plt.plot(t[inoise], rms[inoise], 'c') + p3, = plt.plot(t[isignal],rms[isignal], 'r') p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \ [minsiglevel, minsiglevel], 'g', linewidth=2) - p5, = plt.plot([pick, pick], [min(x), max(x)], 'b', linewidth=2) - plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \ - 'Envelope Signal Window', 'Minimum Signal Level', \ + p5, = plt.plot([pick, pick], [min(rms), max(rms)], 'b', linewidth=2) + plt.legend([p1, p2, p3, p4, p5], ['RMS Data', 'RMS Noise Window', \ + 'RMS Signal Window', 'Minimum Signal Level', \ 'Onset'], loc='best') plt.xlabel('Time [s] since %s' % X[0].stats.starttime) plt.ylabel('Counts') From 1001da728e85713d85b4cf6db6dae492516a3cd2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:17:46 +0200 Subject: [PATCH 06/13] Modified parameter minsiglength for modified routine checksignallength. --- autoPyLoT_local.in | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index 18849c2e..c2ee8508 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -7,8 +7,8 @@ #main settings# /DATA/Insheim #rootpath# %project path EVENT_DATA/LOCAL #datapath# %data path -2013.02_Insheim #database# %name of data base -e0019.048.13 #eventID# %event ID for single event processing +2015.08_Insheim #database# %name of data base +e0013.241.15 #eventID# %event ID for single event processing /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure# %choose data structure 0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything @@ -30,8 +30,8 @@ HYPOSAT #locrt# %location routine used ("HYPO 300 #Qp# %quality factor for P waves 100 #Qs# %quality factor for S waves #common settings picker# -20 #pstart# %start time [s] for calculating CF for P-picking -80 #pstop# %end time [s] for calculating CF for P-picking +15 #pstart# %start time [s] for calculating CF for P-picking +60 #pstop# %end time [s] for calculating CF for P-picking -1.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking 7 #sstop# %end time [s] after P-onset for calculating CF for S-picking 2 20 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] @@ -81,14 +81,14 @@ ARH #algoS# %choose algorithm for S-onset #inital AIC onset# 0.01 0.02 0.04 0.08 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P 0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S -50 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected +10 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected 1.2 #minAICPSNR# %below this SNR the initial P pick is rejected 6 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected #check duration of signal using envelope function# -2.5 #minsiglength# %minimum required length of signal [s] -3 #noisefactor# %noiselevel*noisefactor=threshold -70 #minpercent# %required percentage of samples higher than threshold +5 #minsiglength# %minimum required length of signal [s] +1.8 #noisefactor# %noiselevel*noisefactor=threshold +50 #minpercent# %required percentage of samples higher than threshold #check for spuriously picked S-onsets# 2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude #check statistics of P onsets# From 956e16fdcab32bd2a394ca9ec442d68d7cd23479 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:17:56 +0200 Subject: [PATCH 07/13] Modified parameter minsiglength for modified routine checksignallength. --- autoPyLoT_regional.in | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/autoPyLoT_regional.in b/autoPyLoT_regional.in index abcc0812..2231d007 100644 --- a/autoPyLoT_regional.in +++ b/autoPyLoT_regional.in @@ -9,6 +9,7 @@ EVENT_DATA/LOCAL #datapath# %data path 2006.01_Nisyros #database# %name of data base e1412.008.06 #eventID# %event ID for single event processing +/DATA/Egelados/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure# %choose data structure 0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything AUTOPHASES_AIC_HOS4_ARH #phasefile# %name of autoPILOT output phase file @@ -83,8 +84,8 @@ ARH #algoS# %choose algorithm for S-onset 8 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected #check duration of signal using envelope function# -6 #minsiglength# %minimum required length of signal [s] -1.5 #noisefactor# %noiselevel*noisefactor=threshold +30 #minsiglength# %minimum required length of signal [s] +2.5 #noisefactor# %noiselevel*noisefactor=threshold 60 #minpercent# %required percentage of samples higher than threshold #check for spuriously picked S-onsets# 1.5 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude From 254c745f2524898aff8fc75c1aa829e962c4475b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:19:57 +0200 Subject: [PATCH 08/13] Marginal changes. --- pylot/core/pick/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index f94cac76..c24ce1f8 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -7,7 +7,7 @@ :author: Ludger Kueperkoch / MAGS2 EP3 working group """ -import pdb + import numpy as np import scipy as sc import matplotlib.pyplot as plt From 0753071cfb8298c551842ea3de42664b373f1052 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 11:23:59 +0200 Subject: [PATCH 09/13] Removed import of scipy as this is no more necessary. --- pylot/core/pick/utils.py | 1 - 1 file changed, 1 deletion(-) diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index c24ce1f8..052c2870 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -9,7 +9,6 @@ """ import numpy as np -import scipy as sc import matplotlib.pyplot as plt from obspy.core import Stream, UTCDateTime import warnings From ac7d239b4062daadd3efb871e9d0326e4bb419c5 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 12:09:17 +0200 Subject: [PATCH 10/13] Optimized parameter pstop, as this is a cruicial parameter for stable AIC-CF calculation. --- autoPyLoT_regional.in | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/autoPyLoT_regional.in b/autoPyLoT_regional.in index 2231d007..23cdce87 100644 --- a/autoPyLoT_regional.in +++ b/autoPyLoT_regional.in @@ -31,7 +31,7 @@ HYPOSAT #locrt# %location routine used ("HYPO 100 #Qs# %quality factor for S waves #common settings picker# 20 #pstart# %start time [s] for calculating CF for P-picking -160 #pstop# %end time [s] for calculating CF for P-picking +100 #pstop# %end time [s] for calculating CF for P-picking 3.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking 100 #sstop# %end time [s] after P-onset for calculating CF for S-picking 3 10 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] @@ -79,9 +79,9 @@ ARH #algoS# %choose algorithm for S-onset #inital AIC onset# 0.04 0.08 0.16 0.32 #timeerrorsP# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P 0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S -5 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected +3 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected 1.2 #minAICPSNR# %below this SNR the initial P pick is rejected -8 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected +3 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected #check duration of signal using envelope function# 30 #minsiglength# %minimum required length of signal [s] From cfca52e576025c5cd5dfeb9f4468304eaccd555d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 15:28:37 +0200 Subject: [PATCH 11/13] Debuged slope determination [counts/s] within AICPicker. --- pylot/core/pick/Picker.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index d8a224e6..61eade5e 100644 --- a/pylot/core/pick/Picker.py +++ b/pylot/core/pick/Picker.py @@ -18,6 +18,7 @@ calculated after Diehl & Kissling (2009). :author: MAGS2 EP3 working group / Ludger Kueperkoch """ + import numpy as np import matplotlib.pyplot as plt from pylot.core.pick.utils import getnoisewin, getsignalwin @@ -245,8 +246,7 @@ class AICPicker(AutoPicking): if datafit[0] >= datafit[len(datafit) - 1]: print 'AICPicker: Negative slope, bad onset skipped!' return - - self.slope = 1 / tslope * datafit[len(dataslope) - 1] - datafit[0] + self.slope = 1 / tslope * (datafit[len(dataslope) - 1] - datafit[0]) else: self.SNR = None From 5d8346b1caebb4613ab3359e066305fb79ae89ac Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 15:28:57 +0200 Subject: [PATCH 12/13] Optimized some parameters. --- autoPyLoT_local.in | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index c2ee8508..84ae2137 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -7,8 +7,8 @@ #main settings# /DATA/Insheim #rootpath# %project path EVENT_DATA/LOCAL #datapath# %data path -2015.08_Insheim #database# %name of data base -e0013.241.15 #eventID# %event ID for single event processing +2013.02_Insheim #database# %name of data base +e0019.048.13 #eventID# %event ID for single event processing /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure# %choose data structure 0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything From 70b3f031f8026fed3544c109ec22c71202088948 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 4 Sep 2015 15:29:04 +0200 Subject: [PATCH 13/13] Optimized some parameters. --- autoPyLoT_regional.in | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/autoPyLoT_regional.in b/autoPyLoT_regional.in index 23cdce87..52f5c465 100644 --- a/autoPyLoT_regional.in +++ b/autoPyLoT_regional.in @@ -32,7 +32,7 @@ HYPOSAT #locrt# %location routine used ("HYPO #common settings picker# 20 #pstart# %start time [s] for calculating CF for P-picking 100 #pstop# %end time [s] for calculating CF for P-picking -3.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking +1.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking 100 #sstop# %end time [s] after P-onset for calculating CF for S-picking 3 10 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz] 3 12 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz] @@ -65,11 +65,11 @@ ARH #algoS# %choose algorithm for S-onset 0.3 #tpred2h# %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick 4 #Sarorder# %for AR-picker, order of AR process of H-components 10 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H) -6 #pickwinS# %for initial AIC and refined pick, length of S-pick window [s] +25 #pickwinS# %for initial AIC and refined pick, length of S-pick window [s] 5 0.2 3.0 3.0 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s] -3.0 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s] +3.5 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s] 1.0 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S) -0.4 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S) +0.2 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S) 1.5 #nfacS# %for AR-picker, noise factor for noise level determination (S) %first-motion picker% 1 #minfmweight# %minimum required p weight for first-motion determination @@ -81,16 +81,16 @@ ARH #algoS# %choose algorithm for S-onset 0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S 3 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected 1.2 #minAICPSNR# %below this SNR the initial P pick is rejected -3 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected -1.5 #minAICSSNR# %below this SNR the initial S pick is rejected +5 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected +2.5 #minAICSSNR# %below this SNR the initial S pick is rejected #check duration of signal using envelope function# 30 #minsiglength# %minimum required length of signal [s] 2.5 #noisefactor# %noiselevel*noisefactor=threshold 60 #minpercent# %required percentage of samples higher than threshold #check for spuriously picked S-onsets# -1.5 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude +1.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude #check statistics of P onsets# -35 #mdttolerance# %maximum allowed deviation of P picks from median [s] +45 #mdttolerance# %maximum allowed deviation of P picks from median [s] #wadati check# -2.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram +3.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram