marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

Merge branch 'develop' of ariadne:/data/git/pylot into develop

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2015-09-05 09:41:52 +02:00
commit e9c4987ca0
6 changed files with 147 additions and 120 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
autoPyLoT_local.in
View File

@ -30,8 +30,8 @@ HYPOSAT #locrt# %location routine used ("HYPO
300 #Qp# %quality factor for P waves 300 #Qp# %quality factor for P waves
100 #Qs# %quality factor for S waves 100 #Qs# %quality factor for S waves
#common settings picker# #common settings picker#
20 #pstart# %start time [s] for calculating CF for P-picking 15 #pstart# %start time [s] for calculating CF for P-picking
80 #pstop# %end 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 -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 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] 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# #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.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 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 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 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 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected
#check duration of signal using envelope function# #check duration of signal using envelope function#
2.5 #minsiglength# %minimum required length of signal [s] 5 #minsiglength# %minimum required length of signal [s]
3 #noisefactor# %noiselevel*noisefactor=threshold 1.8 #noisefactor# %noiselevel*noisefactor=threshold
70 #minpercent# %required percentage of samples higher than threshold 50 #minpercent# %required percentage of samples higher than threshold
#check for spuriously picked S-onsets# #check for spuriously picked S-onsets#
2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude 2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude
#check statistics of P onsets# #check statistics of P onsets#

27
autoPyLoT_regional.in
View File

@ -9,6 +9,7 @@
EVENT_DATA/LOCAL #datapath# %data path EVENT_DATA/LOCAL #datapath# %data path
2006.01_Nisyros #database# %name of data base 2006.01_Nisyros #database# %name of data base
e1412.008.06 #eventID# %event ID for single event processing 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 PILOT #datastructure# %choose data structure
0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything 0 #iplot# %flag for plotting: 0 none, 1, partly, >1 everything
AUTOPHASES_AIC_HOS4_ARH #phasefile# %name of autoPILOT output phase file AUTOPHASES_AIC_HOS4_ARH #phasefile# %name of autoPILOT output phase file
@ -30,8 +31,8 @@ HYPOSAT #locrt# %location routine used ("HYPO
100 #Qs# %quality factor for S waves 100 #Qs# %quality factor for S waves
#common settings picker# #common settings picker#
20 #pstart# %start time [s] for calculating CF for P-picking 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 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 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 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] 3 12 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
@ -64,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 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 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) 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] 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) 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) 1.5 #nfacS# %for AR-picker, noise factor for noise level determination (S)
%first-motion picker% %first-motion picker%
1 #minfmweight# %minimum required p weight for first-motion determination 1 #minfmweight# %minimum required p weight for first-motion determination
@ -78,18 +79,18 @@ ARH #algoS# %choose algorithm for S-onset
#inital AIC 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 #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 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 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 5 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected
1.5 #minAICSSNR# %below this SNR 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# #check duration of signal using envelope function#
6 #minsiglength# %minimum required length of signal [s] 30 #minsiglength# %minimum required length of signal [s]
1.5 #noisefactor# %noiselevel*noisefactor=threshold 2.5 #noisefactor# %noiselevel*noisefactor=threshold
60 #minpercent# %required percentage of samples higher than threshold 60 #minpercent# %required percentage of samples higher than threshold
#check for spuriously picked S-onsets# #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# #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# #wadati check#
2.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram 3.0 #wdttolerance# %maximum allowed deviation from Wadati-diagram

4
pylot/core/pick/Picker.py
View File

@ -18,6 +18,7 @@ calculated after Diehl & Kissling (2009).
:author: MAGS2 EP3 working group / Ludger Kueperkoch :author: MAGS2 EP3 working group / Ludger Kueperkoch
""" """
import numpy as np import numpy as np
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from pylot.core.pick.utils import getnoisewin, getsignalwin from pylot.core.pick.utils import getnoisewin, getsignalwin
@ -245,8 +246,7 @@ class AICPicker(AutoPicking):
if datafit[0] >= datafit[len(datafit) - 1]: if datafit[0] >= datafit[len(datafit) - 1]:
print 'AICPicker: Negative slope, bad onset skipped!' print 'AICPicker: Negative slope, bad onset skipped!'
return return
self.slope = 1 / tslope * (datafit[len(dataslope) - 1] - datafit[0])
self.slope = 1 / tslope * datafit[len(dataslope) - 1] - datafit[0]
else: else:
self.SNR = None self.SNR = None

82
pylot/core/pick/autopick.py
View File

@ -41,9 +41,9 @@ def autopickevent(data, param):
# quality control # quality control
# median check and jackknife on P-onset times # 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) # check S-P times (Wadati)
return wadaticheck(jk_checked_onsets, wdttolerance, 2) return wadaticheck(jk_checked_onsets, wdttolerance, iplot)
def autopickstation(wfstream, pickparam): def autopickstation(wfstream, pickparam):
""" """
@ -196,10 +196,36 @@ def autopickstation(wfstream, pickparam):
############################################################## ##############################################################
if aicpick.getpick() is not None: if aicpick.getpick() is not None:
# check signal length to detect spuriously picked noise peaks # 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 z_copy[0].data = tr_filt.data
Pflag = checksignallength(z_copy, aicpick.getpick(), tsnrz, zne = z_copy
minsiglength, \ if len(ndat) == 0 or len(edat) == 0:
nfacsl, minpercent, iplot) 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: if Pflag == 1:
# check for spuriously picked S onset # check for spuriously picked S onset
# both horizontal traces needed # both horizontal traces needed
@ -207,20 +233,6 @@ def autopickstation(wfstream, pickparam):
print 'One or more horizontal components missing!' print 'One or more horizontal components missing!'
print 'Skipping control function checkZ4S.' print 'Skipping control function checkZ4S.'
else: 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, \ Pflag = checkZ4S(zne, aicpick.getpick(), zfac, \
tsnrz[3], iplot) tsnrz[3], iplot)
if Pflag == 0: if Pflag == 0:
@ -515,18 +527,19 @@ def autopickstation(wfstream, pickparam):
hdat = edat.copy() hdat = edat.copy()
hdat += ndat hdat += ndat
h_copy = hdat.copy() h_copy = hdat.copy()
cordat = data.restituteWFData(invdir, h_copy) [cordat, restflag] = data.restituteWFData(invdir, h_copy)
# calculate WA-peak-to-peak amplitude # calculate WA-peak-to-peak amplitude
# using subclass WApp of superclass Magnitude # using subclass WApp of superclass Magnitude
if Sweight < 4: if restflag == 1:
wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) if Sweight < 4:
else: wapp = WApp(cordat, mpickS, mpickP + sstop, iplot)
# use larger window for getting peak-to-peak amplitude else:
# as the S pick is quite unsure # use larger window for getting peak-to-peak amplitude
wapp = WApp(cordat, mpickP, mpickP + sstop + \ # as the S pick is quite unsure
(0.5 * (mpickP + sstop)), iplot) wapp = WApp(cordat, mpickP, mpickP + sstop + \
(0.5 * (mpickP + sstop)), iplot)
Ao = wapp.getwapp() Ao = wapp.getwapp()
else: else:
print 'Bad initial (AIC) S-pick, skipping this onset!' print 'Bad initial (AIC) S-pick, skipping this onset!'
@ -544,12 +557,13 @@ def autopickstation(wfstream, pickparam):
hdat = edat.copy() hdat = edat.copy()
hdat += ndat hdat += ndat
h_copy = hdat.copy() h_copy = hdat.copy()
cordat = data.restituteWFData(invdir, h_copy) [cordat, restflag] = data.restituteWFData(invdir, h_copy)
# calculate WA-peak-to-peak amplitude if restflag == 1:
# using subclass WApp of superclass Magnitude # calculate WA-peak-to-peak amplitude
wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP \ # using subclass WApp of superclass Magnitude
+ sstop)), iplot) wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP \
Ao = wapp.getwapp() + sstop)), iplot)
Ao = wapp.getwapp()
else: else:
print 'autopickstation: No horizontal component data available or ' \ print 'autopickstation: No horizontal component data available or ' \

132
pylot/core/pick/utils.py
View File

@ -9,7 +9,6 @@
""" """
import numpy as np import numpy as np
import scipy as sc
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
import warnings import warnings
@ -44,7 +43,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
LPick = None LPick = None
EPick = None EPick = None
PickError = 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 x = X[0].data
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, 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) ilup, = np.where(x[isignal] > nlevel)
ildown, = np.where(x[isignal] < -nlevel) ildown, = np.where(x[isignal] < -nlevel)
if not ilup.size and not ildown.size: if not ilup.size and not ildown.size:
print 'earllatepicker: Signal lower than noise level!' print ("earllatepicker: Signal lower than noise level!")
print 'Skip this trace!' print ("Skip this trace!")
return LPick, EPick, PickError return LPick, EPick, PickError
il = min(np.min(ilup) if ilup.size else float('inf'), il = min(np.min(ilup) if ilup.size else float('inf'),
np.min(ildown) if ildown.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 FM = None
if Pick is not 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 xraw = Xraw[0].data
xfilt = Xfilt[0].data xfilt = Xfilt[0].data
@ -182,15 +181,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
else: else:
li1 = index1[0] li1 = index1[0]
if np.size(xraw[ipick[0][1]:ipick[0][li1]]) == 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 P1 = None
else: else:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]])) imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][li1]]))
if imax1 == 0: if imax1 == 0:
imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][index1[1]]])) imax1 = np.argmax(abs(xraw[ipick[0][1]:ipick[0][index1[1]]]))
if imax1 == 0: if imax1 == 0:
print 'fmpicker: Zero crossings too close!' print ("fmpicker: Zero crossings too close!")
print 'Skip first motion determination!' print ("Skip first motion determination!")
return FM return FM
islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1])) islope1 = np.where((t >= Pick) & (t <= Pick + t[imax1]))
@ -224,15 +223,15 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None):
else: else:
li2 = index2[0] li2 = index2[0]
if np.size(xfilt[ipick[0][1]:ipick[0][li2]]) == 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 P2 = None
else: else:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]])) imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][li2]]))
if imax2 == 0: if imax2 == 0:
imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][index2[1]]])) imax2 = np.argmax(abs(xfilt[ipick[0][1]:ipick[0][index2[1]]]))
if imax2 == 0: if imax2 == 0:
print 'fmpicker: Zero crossings too close!' print ("fmpicker: Zero crossings too close!")
print 'Skip first motion determination!' print ("Skip first motion determination!")
return FM return FM
islope2 = np.where((t >= Pick) & (t <= Pick + t[imax2])) 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: elif P1[0] > 0 and P2[0] <= 0:
FM = '+' FM = '+'
print 'fmpicker: Found polarity %s' % FM print ("fmpicker: Found polarity %s" % FM)
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)
@ -331,10 +330,10 @@ def getSNR(X, TSNR, t1):
# get signal window # get signal window
isignal = getsignalwin(t, t1, TSNR[2]) isignal = getsignalwin(t, t1, TSNR[2])
if np.size(inoise) < 1: if np.size(inoise) < 1:
print 'getSNR: Empty array inoise, check noise window!' print ("getSNR: Empty array inoise, check noise window!")
return return
elif np.size(isignal) < 1: elif np.size(isignal) < 1:
print 'getSNR: Empty array isignal, check signal window!' print ("getSNR: Empty array isignal, check signal window!")
return return
# demean over entire waveform # demean over entire waveform
@ -372,7 +371,7 @@ def getnoisewin(t, t1, tnoise, tgap):
inoise, = np.where((t <= max([t1 - tgap, 0])) \ inoise, = np.where((t <= max([t1 - tgap, 0])) \
& (t >= max([t1 - tnoise - tgap, 0]))) & (t >= max([t1 - tnoise - tgap, 0])))
if np.size(inoise) < 1: if np.size(inoise) < 1:
print 'getnoisewin: Empty array inoise, check noise window!' print ("getnoisewin: Empty array inoise, check noise window!")
return inoise return inoise
@ -396,7 +395,7 @@ def getsignalwin(t, t1, tsignal):
isignal, = np.where((t <= min([t1 + tsignal, len(t)])) \ isignal, = np.where((t <= min([t1 + tsignal, len(t)])) \
& (t >= t1)) & (t >= t1))
if np.size(isignal) < 1: if np.size(isignal) < 1:
print 'getsignalwin: Empty array isignal, check signal window!' print ("getsignalwin: Empty array isignal, check signal window!")
return isignal return isignal
@ -483,8 +482,8 @@ def wadaticheck(pickdic, dttolerance, iplot):
# calculate vp/vs ratio before check # calculate vp/vs ratio before check
vpvsr = p1[0] + 1 vpvsr = p1[0] + 1
print '###############################################' print ("###############################################")
print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr print ("wadaticheck: Average Vp/Vs ratio before check: %f" % vpvsr)
checkedPpicks = [] checkedPpicks = []
checkedSpicks = [] checkedSpicks = []
@ -521,18 +520,18 @@ def wadaticheck(pickdic, dttolerance, iplot):
# calculate vp/vs ratio after check # calculate vp/vs ratio after check
cvpvsr = p2[0] + 1 cvpvsr = p2[0] + 1
print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr print ("wadaticheck: Average Vp/Vs ratio after check: %f" % cvpvsr)
print 'wadatacheck: Skipped %d S pick(s).' % ibad print ("wadatacheck: Skipped %d S pick(s)" % ibad)
else: else:
print '###############################################' print ("###############################################")
print 'wadatacheck: Not enough checked S-P times available!' print ("wadatacheck: Not enough checked S-P times available!")
print 'Skip Wadati check!' print ("Skip Wadati check!")
checkedonsets = pickdic checkedonsets = pickdic
else: else:
print 'wadaticheck: Not enough S-P times available for reliable regression!' print ("wadaticheck: Not enough S-P times available for reliable regression!")
print 'Skip wadati check!' print ("Skip wadati check!")
wfitflag = 1 wfitflag = 1
# plot results # plot results
@ -562,9 +561,9 @@ def wadaticheck(pickdic, dttolerance, iplot):
def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot): def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
''' '''
Function to detect spuriously picked noise peaks. Function to detect spuriously picked noise peaks.
Uses envelope to determine, how many samples [per cent] after Uses RMS trace of all 3 components (if available) to determine,
P onset are below certain threshold, calculated from noise how many samples [per cent] after P onset are below certain
level times noise factor. threshold, calculated from noise level times noise factor.
: param: X, time series (seismogram) : param: X, time series (seismogram)
: type: `~obspy.core.stream.Stream` : type: `~obspy.core.stream.Stream`
@ -592,47 +591,54 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
assert isinstance(X, Stream), "%s is not a stream object" % str(X) assert isinstance(X, Stream), "%s is not a stream object" % str(X)
print 'Checking signal length ...' print ("Checking signal length ...")
x = X[0].data if len(X) > 1:
t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, # 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) X[0].stats.delta)
# generate envelope function from Hilbert transform # get noise window in front of pick plus saftey gap
y = np.imag(sc.signal.hilbert(x)) inoise = getnoisewin(t, pick - 0.5, TSNR[0], TSNR[1])
e = np.sqrt(np.power(x, 2) + np.power(y, 2))
# get noise window
inoise = getnoisewin(t, pick, TSNR[0], TSNR[1])
# get signal window # get signal window
isignal = getsignalwin(t, pick, TSNR[2]) isignal = getsignalwin(t, pick, minsiglength)
# calculate minimum adjusted signal level # calculate minimum adjusted signal level
minsiglevel = max(e[inoise]) * nfac minsiglevel = max(rms[inoise]) * nfac
# minimum adjusted number of samples over minimum signal level # minimum adjusted number of samples over minimum signal level
minnum = len(isignal) * minpercent/100 minnum = len(isignal) * minpercent/100
# get number of samples above minimum adjusted signal level # 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: if numoverthr >= minnum:
print 'checksignallength: Signal reached required length.' print ("checksignallength: Signal reached required length.")
returnflag = 1 returnflag = 1
else: else:
print 'checksignallength: Signal shorter than required minimum signal length!' print ("checksignallength: Signal shorter than required minimum signal length!")
print 'Presumably picked noise peak, pick is rejected!' 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 returnflag = 0
if iplot == 2: if iplot == 2:
plt.figure(iplot) plt.figure(iplot)
p1, = plt.plot(t,x, 'k') p1, = plt.plot(t,rms, 'k')
p2, = plt.plot(t[inoise], e[inoise], 'c') p2, = plt.plot(t[inoise], rms[inoise], 'c')
p3, = plt.plot(t[isignal],e[isignal], 'r') p3, = plt.plot(t[isignal],rms[isignal], 'r')
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]]], \ 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) p5, = plt.plot([pick, pick], [min(rms), max(rms)], 'b', linewidth=2)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \ plt.legend([p1, p2, p3, p4, p5], ['RMS Data', 'RMS Noise Window', \
'Envelope Signal Window', 'Minimum Signal Level', \ 'RMS Signal Window', 'Minimum Signal Level', \
'Onset'], loc='best') 'Onset'], loc='best')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime) plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.ylabel('Counts') plt.ylabel('Counts')
@ -675,8 +681,8 @@ def checkPonsets(pickdic, dttolerance, iplot):
stations.append(key) stations.append(key)
# apply jackknife bootstrapping on variance of P onsets # apply jackknife bootstrapping on variance of P onsets
print '###############################################' print ("###############################################")
print 'checkPonsets: Apply jackknife bootstrapping on P-onset times ...' print ("checkPonsets: Apply jackknife bootstrapping on P-onset times ...")
[xjack,PHI_pseudo,PHI_sub] = jackknife(Ppicks, 'VAR', 1) [xjack,PHI_pseudo,PHI_sub] = jackknife(Ppicks, 'VAR', 1)
# get pseudo variances smaller than average variances # get pseudo variances smaller than average variances
# (times safety factor), these picks passed jackknife test # (times safety factor), these picks passed jackknife test
@ -684,7 +690,7 @@ def checkPonsets(pickdic, dttolerance, iplot):
# these picks did not pass jackknife test # these picks did not pass jackknife test
badjk = np.where(PHI_pseudo > 2 * xjack) badjk = np.where(PHI_pseudo > 2 * xjack)
badjkstations = np.array(stations)[badjk] 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 # calculate median from these picks
pmedian = np.median(np.array(Ppicks)[ij]) pmedian = np.median(np.array(Ppicks)[ij])
@ -696,9 +702,9 @@ def checkPonsets(pickdic, dttolerance, iplot):
goodstations = np.array(stations)[igood] goodstations = np.array(stations)[igood]
badstations = np.array(stations)[ibad] badstations = np.array(stations)[ibad]
print 'checkPonsets: %d pick(s) deviate too much from median!' % len(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) \ print ("checkPonsets: Skipped %d P pick(s) out of %d" % (len(badstations) \
+ len(badjkstations), len(stations)) + len(badjkstations), len(stations)))
goodmarker = 'goodPonsetcheck' goodmarker = 'goodPonsetcheck'
badmarker = 'badPonsetcheck' badmarker = 'badPonsetcheck'
@ -765,8 +771,8 @@ def jackknife(X, phi, h):
g = len(X) / h g = len(X) / h
if type(g) is not int: if type(g) is not int:
print 'jackknife: Cannot divide quantity X in equal sized subgroups!' print ("jackknife: Cannot divide quantity X in equal sized subgroups!")
print 'Choose another size for subgroups!' print ("Choose another size for subgroups!")
return PHI_jack, PHI_pseudo, PHI_sub return PHI_jack, PHI_pseudo, PHI_sub
else: else:
# estimator of undisturbed spot check # estimator of undisturbed spot check
@ -834,7 +840,7 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
assert isinstance(X, Stream), "%s is not a stream object" % str(X) 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 returnflag = 0
@ -875,9 +881,9 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
# vertical P-coda level must exceed horizontal P-coda level # vertical P-coda level must exceed horizontal P-coda level
# zfac times encodalevel # zfac times encodalevel
if zcodalevel < minsiglevel: if zcodalevel < minsiglevel:
print 'checkZ4S: Maybe S onset? Skip this P pick!' print ("checkZ4S: Maybe S onset? Skip this P pick!")
else: else:
print 'checkZ4S: P onset passes checkZ4S test!' print ("checkZ4S: P onset passes checkZ4S test!")
returnflag = 1 returnflag = 1
if iplot > 1: if iplot > 1:

10
pylot/core/read/data.py
View File

@ -228,6 +228,9 @@ class Data(object):
:param streams: :param streams:
:return: :return:
""" """
restflag = 0
if streams is None: if streams is None:
st_raw = self.getWFData() st_raw = self.getWFData()
st = st_raw.copy() st = st_raw.copy()
@ -236,7 +239,7 @@ class Data(object):
for tr in st: for tr in st:
# remove underscores # 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] tr.stats.station = tr.stats.station[0:3]
dlp = '%s/*.dless' % invdlpath dlp = '%s/*.dless' % invdlpath
invp = '%s/*.xml' % invdlpath invp = '%s/*.xml' % invdlpath
@ -273,6 +276,7 @@ class Data(object):
'date': st[ 'date': st[
i].stats.starttime, i].stats.starttime,
'units': "VEL"}) 'units': "VEL"})
restflag = 1
except ValueError as e: except ValueError as e:
vmsg = '{0}'.format(e) vmsg = '{0}'.format(e)
print(vmsg) print(vmsg)
@ -303,6 +307,7 @@ class Data(object):
st[i].attach_response(inv) st[i].attach_response(inv)
st[i].remove_response(output='VEL', st[i].remove_response(output='VEL',
pre_filt=prefilt) pre_filt=prefilt)
restflag = 1
except ValueError as e: except ValueError as e:
vmsg = '{0}'.format(e) vmsg = '{0}'.format(e)
print(vmsg) print(vmsg)
@ -334,6 +339,7 @@ class Data(object):
'units': "VEL"} 'units': "VEL"}
st[i].simulate(paz_remove=None, pre_filt=prefilt, st[i].simulate(paz_remove=None, pre_filt=prefilt,
seedresp=seedresp) seedresp=seedresp)
restflag = 1
except ValueError as e: except ValueError as e:
vmsg = '{0}'.format(e) vmsg = '{0}'.format(e)
print(vmsg) print(vmsg)
@ -346,7 +352,7 @@ class Data(object):
print("Go on processing data without source parameter " print("Go on processing data without source parameter "
"determination!") "determination!")
return st return st, restflag
def getEvtData(self): def getEvtData(self):
""" """