diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 1d620f47..5dc320ec 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -10,7 +10,7 @@ import numpy as np import matplotlib.pyplot as plt -from obspy.core import Stream +from obspy.core import Stream, UTCDateTime def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): @@ -66,9 +66,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): #get earliest possible pick #determine all zero crossings in signal window - # remove mean from signal window - signal = x[isignal] - x[isignal].mean() - zc = crossings_nonzero_all(signal) + zc = crossings_nonzero_all(x[isignal]) #calculate mean half period T0 of signal as the average of the T0 = np.mean(np.diff(zc)) * X[0].stats.delta #this is half wave length! #T0/4 is assumed as time difference between most likely and earliest possible pick! @@ -161,7 +159,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): index1 = [] i = 0 for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]): - i += 1 + i = i + 1 if xraw[j - 1] <= 0 and xraw[j] >= 0: zc1.append(t[ipick][i]) index1.append(i) @@ -196,7 +194,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=None): index2 = [] i = 0 for j in range(ipick[0][1], ipick[0][len(t[ipick]) - 1]): - i += 1 + i = i + 1 if xfilt[j - 1] <= 0 and xfilt[j] >= 0: zc2.append(t[ipick][i]) index2.append(i) @@ -382,4 +380,117 @@ def getsignalwin(t, t1, tsignal): return isignal +def wadaticheck(pickdic, dttolerance, iplot): + ''' + Function to calculate Wadati-diagram from given P and S onsets in order + to detect S pick outliers. If a certain S-P time deviates from regression + of S-P time the S pick is marked and down graded. + : param: pickdic, dictionary containing picks and quality parameters + : type: dictionary + + : param: dttolerance, maximum adjusted deviation of S-P time from + S-P time regression + : type: float + + : param: iplot, if iplot > 1, Wadati diagram is shown + : type: int + ''' + + checkedonsets = pickdic + + # search for good quality picks and calculate S-P time + Ppicks = [] + Spicks = [] + SPtimes = [] + vpvs = [] + for key in pickdic: + if pickdic[key]['P']['weight'] < 4 and pickdic[key]['S']['weight'] < 4: + # calculate S-P time + spt = UTCDateTime(pickdic[key]['S']['mpp']) - UTCDateTime(pickdic[key]['P']['mpp']) + # add S-P time to dictionary + pickdic[key]['SPt'] = spt + # add P onsets and corresponding S-P times to list + UTCPpick = UTCDateTime(pickdic[key]['P']['mpp']) - UTCDateTime(1970,1,1,0,0,0) + UTCSpick = UTCDateTime(pickdic[key]['S']['mpp']) - UTCDateTime(1970,1,1,0,0,0) + Ppicks.append(UTCPpick) + Spicks.append(UTCSpick) + SPtimes.append(spt) + vpvs.append(UTCPpick/UTCSpick) + + # calculate average vp/vs ratio before check + vpvsr = np.mean(vpvs) + print 'wadaticheck: Average Vp/Vs ratio before check:', vpvsr + + if len(SPtimes) >= 3: + # calculate slope + p1 = np.polyfit(Ppicks, SPtimes, 1) + wdfit = np.polyval(p1, Ppicks) + wfitflag = 0 + + checkedPpicks = [] + checkedSpicks = [] + checkedSPtimes = [] + checkedvpvs = [] + # calculate deviations from Wadati regression + for key in pickdic: + if pickdic[key].has_key('SPt'): + ii = 0 + wddiff = abs(pickdic[key]['SPt'] - wdfit[ii]) + ii += 1 + # check, if deviation is larger than adjusted + if wddiff >= dttolerance: + # mark onset and downgrade S-weight to 4 + # (not used anymore) + marker = 'badWadatiCheck' + pickdic[key]['S']['weight'] = 4 + else: + marker = 'goodWadatiCheck' + checkedPpick = UTCDateTime(pickdic[key]['P']['mpp']) - \ + UTCDateTime(1970,1,1,0,0,0) + checkedPpicks.append(checkedPpick) + checkedSpick = UTCDateTime(pickdic[key]['S']['mpp']) - \ + UTCDateTime(1970,1,1,0,0,0) + checkedSpicks.append(checkedSpick) + checkedSPtime = UTCDateTime(pickdic[key]['S']['mpp']) - \ + UTCDateTime(pickdic[key]['P']['mpp']) + checkedSPtimes.append(checkedSPtime) + checkedvpvs.append(checkedPpick/checkedSpick) + + pickdic[key]['S']['marked'] = marker + + # calculate average vp/vs ratio after check + cvpvsr = np.mean(checkedvpvs) + print 'wadaticheck: Average Vp/Vs ratio after check:', cvpvsr + + # calculate new slope + p2 = np.polyfit(checkedPpicks, checkedSPtimes, 1) + wdfit2 = np.polyval(p2, checkedPpicks) + + checkedonsets = pickdic + + else: + print 'wadaticheck: Not enough S-P times available for reliable regression!' + print 'Skip wadati check!' + wfitflag = 1 + + # plot results + iplot = 2 + if iplot > 1: + f = plt.figure(iplot) + f1, = plt.plot(Ppicks, SPtimes, 'ro') + if wfitflag == 0: + f2, = plt.plot(Ppicks, wdfit, 'k') + f3, = plt.plot(checkedPpicks, checkedSPtimes, 'ko') + f4, = plt.plot(checkedPpicks, wdfit2, 'g') + plt.ylabel('S-P Times [s]') + plt.xlabel('P Times [s]') + plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(old)=%5.2f, Vp/Vs(checked)=%5.2f' \ + % (len(SPtimes), vpvsr, cvpvsr)) + plt.legend([f1, f2, f3, f4], ['Skipped S-Picks', 'Wadati 1', 'Reliable S-Picks', \ + 'Wadati 2'], loc='best') + plt.show() + raw_input() + plt.close(f) + + return checkedonsets