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

Modified earllatepicker: Mean is removed from trace calculated from noise + signal window.

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-06-22 12:39:29 +02:00
parent 635ac1686b
commit aba3997b20
1 changed files with 36 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -53,6 +53,9 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1]) inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1])
# get signal window # get signal window
isignal = getsignalwin(t, Pick1, TSNR[2]) isignal = getsignalwin(t, Pick1, TSNR[2])
# remove mean
meanwin = np.hstack((inoise, isignal))
x = x - np.mean(x[meanwin])
# calculate noise level # calculate noise level
nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac
# get time where signal exceeds nlevel # get time where signal exceeds nlevel
@ -412,34 +415,28 @@ def wadaticheck(pickdic, dttolerance, iplot):
SPtimes = [] SPtimes = []
for key in pickdic: for key in pickdic:
if pickdic[key]['P']['weight'] < 4 and pickdic[key]['S']['weight'] < 4: if pickdic[key]['P']['weight'] < 4 and pickdic[key]['S']['weight'] < 4:
# calculate S-P time # calculate S-P time
spt = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp'] spt = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp']
# add S-P time to dictionary # add S-P time to dictionary
pickdic[key]['SPt'] = spt pickdic[key]['SPt'] = spt
# add P onsets and corresponding S-P times to list # add P onsets and corresponding S-P times to list
UTCPpick = UTCDateTime(pickdic[key]['P']['mpp']) - UTCDateTime(1970, UTCPpick = UTCDateTime(pickdic[key]['P']['mpp'])
1, 1, UTCSpick = UTCDateTime(pickdic[key]['S']['mpp'])
0, 0, Ppicks.append(UTCPpick.timestamp)
0) Spicks.append(UTCSpick.timestamp)
UTCSpick = UTCDateTime(pickdic[key]['S']['mpp']) - UTCDateTime(1970, SPtimes.append(spt)
1, 1,
0, 0,
0)
Ppicks.append(UTCPpick)
Spicks.append(UTCSpick)
SPtimes.append(spt)
vpvs.append(UTCPpick / UTCSpick)
if len(SPtimes) >= 3: if len(SPtimes) >= 3:
# calculate slope # calculate slope
p1 = np.polyfit(Ppicks, SPtimes, 1) p1 = np.polyfit(Ppicks, SPtimes, 1)
wdfit = np.polyval(p1, Ppicks) wdfit = np.polyval(p1, Ppicks)
wfitflag = 0 wfitflag = 0
# calculate vp/vs ratio before check # calculate vp/vs ratio before check
vpvsr = p1[0] + 1 vpvsr = p1[0] + 1
print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr
checkedPpicks = [] checkedPpicks = []
checkedSpicks = [] checkedSpicks = []
checkedSPtimes = [] checkedSPtimes = []
@ -457,23 +454,19 @@ def wadaticheck(pickdic, dttolerance, iplot):
pickdic[key]['S']['weight'] = 9 pickdic[key]['S']['weight'] = 9
else: else:
marker = 'goodWadatiCheck' marker = 'goodWadatiCheck'
checkedPpick = UTCDateTime(pickdic[key]['P']['mpp']) - \ checkedPpick = UTCDateTime(pickdic[key]['P']['mpp'])
UTCDateTime(1970, 1, 1, 0, 0, 0) checkedPpicks.append(checkedPpick.timestamp)
checkedPpicks.append(checkedPpick) checkedSpick = UTCDateTime(pickdic[key]['S']['mpp'])
checkedSpick = UTCDateTime(pickdic[key]['S']['mpp']) - \ checkedSpicks.append(checkedSpick.timestamp)
UTCDateTime(1970, 1, 1, 0, 0, 0) checkedSPtime = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp']
checkedSpicks.append(checkedSpick)
checkedSPtime = pickdic[key]['S']['mpp'] - \
pickdic[key]['P']['mpp']
checkedSPtimes.append(checkedSPtime) checkedSPtimes.append(checkedSPtime)
checkedvpvs.append(checkedPpick / checkedSpick)
pickdic[key]['S']['marked'] = marker pickdic[key]['S']['marked'] = marker
# calculate new slope # calculate new slope
p2 = np.polyfit(checkedPpicks, checkedSPtimes, 1) p2 = np.polyfit(checkedPpicks, checkedSPtimes, 1)
wdfit2 = np.polyval(p2, checkedPpicks) wdfit2 = np.polyval(p2, checkedPpicks)
# calculate vp/vs ratio after check # calculate vp/vs ratio after check
cvpvsr = p2[0] + 1 cvpvsr = p2[0] + 1
@ -482,26 +475,24 @@ def wadaticheck(pickdic, dttolerance, iplot):
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
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)
f1, = plt.plot(Ppicks, SPtimes, 'ro') f1, = plt.plot(Ppicks, SPtimes, 'ro')
if wfitflag == 0: if wfitflag == 0:
f2, = plt.plot(Ppicks, wdfit, 'k') f2, = plt.plot(Ppicks, wdfit, 'k')
f3, = plt.plot(checkedPpicks, checkedSPtimes, 'ko') f3, = plt.plot(checkedPpicks, checkedSPtimes, 'ko')
f4, = plt.plot(checkedPpicks, wdfit2, 'g') f4, = plt.plot(checkedPpicks, wdfit2, 'g')
plt.ylabel('S-P Times [s]') plt.ylabel('S-P Times [s]')
plt.xlabel('P Times [s]') plt.xlabel('P Times [s]')
plt.title( plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(raw)=%5.2f, Vp/Vs(checked)=%5.2f' \
'Wadati-Diagram, %d S-P Times, Vp/Vs(old)=%5.2f, Vp/Vs(checked)=%5.2f' \ % (len(SPtimes), vpvsr, cvpvsr))
% (len(SPtimes), vpvsr, cvpvsr)) plt.legend([f1, f2, f3, f4], ['Skipped S-Picks', 'Wadati 1', 'Reliable S-Picks', \
plt.legend([f1, f2, f3, f4], 'Wadati 2'], loc='best')
['Skipped S-Picks', 'Wadati 1', 'Reliable S-Picks', \
'Wadati 2'], loc='best')
plt.show() plt.show()
raw_input() raw_input()
plt.close(iplot) plt.close(iplot)