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

new function getResolutionWindow and doc testing added

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2015-06-25 10:21:52 +02:00
parent c9f07b6540
commit 7ec28664b4
1 changed files with 49 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -398,6 +398,43 @@ def getsignalwin(t, t1, tsignal):
return isignal return isignal
def getResolutionWindow(snr):
"""
Number -> Float
produce the half of the time resolution window width from given SNR
value
SNR >= 3 -> 2 sec HRW
3 > SNR >= 2 -> 5 sec MRW
2 > SNR >= 1.5 -> 10 sec LRW
1.5 > SNR -> 15 sec VLRW
see also Diehl et al. 2009
>>> getResolutionWindow(0.5)
7.5
>>> getResolutionWindow(1.8)
5.0
>>> getResolutionWindow(2.3)
2.5
>>> getResolutionWindow(4)
1.0
>>> getResolutionWindow(2)
2.5
"""
res_wins = {'HRW': 2., 'MRW': 5., 'LRW': 10., 'VLRW': 15.}
if snr < 1.5:
time_resolution = res_wins['VLRW']
elif snr < 2.:
time_resolution = res_wins['LRW']
elif snr < 3.:
time_resolution = res_wins['MRW']
else:
time_resolution = res_wins['HRW']
return time_resolution/2
def wadaticheck(pickdic, dttolerance, iplot): def wadaticheck(pickdic, dttolerance, iplot):
''' '''
Function to calculate Wadati-diagram from given P and S onsets in order Function to calculate Wadati-diagram from given P and S onsets in order
@ -453,7 +490,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
for key in pickdic: for key in pickdic:
if pickdic[key].has_key('SPt'): if pickdic[key].has_key('SPt'):
wddiff = abs(pickdic[key]['SPt'] - wdfit[ii]) wddiff = abs(pickdic[key]['SPt'] - wdfit[ii])
ii += 1 ii += 1
# check, if deviation is larger than adjusted # check, if deviation is larger than adjusted
if wddiff >= dttolerance: if wddiff >= dttolerance:
# mark onset and downgrade S-weight to 9 # mark onset and downgrade S-weight to 9
@ -481,7 +518,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr
else: else:
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
@ -489,7 +526,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
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
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)
@ -518,14 +555,14 @@ 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 envelope to determine, how many samples [per cent] after
P onset are below certain threshold, calculated from noise P onset are below certain threshold, calculated from noise
level times noise factor. level times noise factor.
: param: X, time series (seismogram) : param: X, time series (seismogram)
: type: `~obspy.core.stream.Stream` : type: `~obspy.core.stream.Stream`
: param: pick, initial (AIC) P onset time : param: pick, initial (AIC) P onset time
: type: float : type: float
: param: TSNR, length of time windows around initial pick [s] : param: TSNR, length of time windows around initial pick [s]
: type: tuple (T_noise, T_gap, T_signal) : type: tuple (T_noise, T_gap, T_signal)
@ -544,7 +581,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
: param: iplot, if iplot > 1, results are shown in figure : param: iplot, if iplot > 1, results are shown in figure
: type: int : type: int
''' '''
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 ...'
@ -565,7 +602,7 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
# 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(e[isignal] >= minsiglevel)[0])
if numoverthr >= minnum: if numoverthr >= minnum:
print 'checksignallength: Signal reached required length.' print 'checksignallength: Signal reached required length.'
@ -574,12 +611,12 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
print 'checksignallength: Signal shorter than required minimum signal length!' print 'checksignallength: Signal shorter than required minimum signal length!'
print 'Presumably picked picked noise peak, pick is rejected!' print 'Presumably picked picked noise peak, pick is rejected!'
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,x, 'k')
p2, = plt.plot(t[inoise], e[inoise]) p2, = plt.plot(t[inoise], e[inoise])
p3, = plt.plot(t[isignal],e[isignal], 'r') 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')
p5, = plt.plot([pick, pick], [min(x), max(x)], 'c') p5, = plt.plot([pick, pick], [min(x), max(x)], 'c')
@ -596,4 +633,6 @@ def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
return returnflag return returnflag
if __name__ == '__main__':
import doctest
doctest.testmod()