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

Uses now UTCDateTime.timestamp as this is more efficient and shorter.

Browse Source
This commit is contained in:
Ludger Küperkoch 2015-06-22 11:07:22 +02:00
parent 6b14c452e2
commit fd6e4cb02a
1 changed files with 24 additions and 29 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -13,6 +13,7 @@ import matplotlib.pyplot as plt
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
import warnings import warnings
def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
''' '''
Function to derive earliest and latest possible pick after Diehl & Kissling (2009) Function to derive earliest and latest possible pick after Diehl & Kissling (2009)
@ -65,8 +66,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
#get earliest possible pick #get earliest possible pick
#determine all zero crossings in signal window (demeaned) #determine all zero crossings in signal window
zc = crossings_nonzero_all(x[isignal] - x[isignal].mean()) zc = crossings_nonzero_all(x[isignal])
#calculate mean half period T0 of signal as the average of the #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 = 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! #T0/4 is assumed as time difference between most likely and earliest possible pick!
@ -385,13 +386,13 @@ def getsignalwin(t, t1, tsignal):
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
to detect S pick outliers. If a certain S-P time deviates from regression to detect S pick outliers. If a certain S-P time deviates by dttolerance
of S-P time the S pick is marked and down graded. from regression of S-P time the S pick is marked and down graded.
: param: pickdic, dictionary containing picks and quality parameters : param: pickdic, dictionary containing picks and quality parameters
: type: dictionary : type: dictionary
: param: dttolerance, maximum adjusted deviation of S-P time from : param: dttolerance, maximum adjusted deviation of S-P time from
S-P time regression S-P time regression
: type: float : type: float
@ -405,7 +406,6 @@ def wadaticheck(pickdic, dttolerance, iplot):
Ppicks = [] Ppicks = []
Spicks = [] Spicks = []
SPtimes = [] SPtimes = []
vpvs = []
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
@ -413,65 +413,60 @@ def wadaticheck(pickdic, dttolerance, iplot):
# 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,1,1,0,0,0) UTCPpick = UTCDateTime(pickdic[key]['P']['mpp'])
UTCSpick = UTCDateTime(pickdic[key]['S']['mpp']) - UTCDateTime(1970,1,1,0,0,0) UTCSpick = UTCDateTime(pickdic[key]['S']['mpp'])
Ppicks.append(UTCPpick) Ppicks.append(UTCPpick.timestamp)
Spicks.append(UTCSpick) Spicks.append(UTCSpick.timestamp)
SPtimes.append(spt) 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 average 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 = []
checkedvpvs = []
# calculate deviations from Wadati regression # calculate deviations from Wadati regression
for key in pickdic: for key in pickdic:
if pickdic[key].has_key('SPt'): if pickdic[key].has_key('SPt'):
ii = 0 ii = 0
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
# (not used anymore) # (not used anymore)
marker = 'badWadatiCheck' marker = 'badWadatiCheck'
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)
checkedSpicks.append(checkedSpick)
checkedSPtime = pickdic[key]['S']['mpp'] - pickdic[key]['P']['mpp'] 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 average 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:', cvpvsr
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!'
@ -487,7 +482,7 @@ def wadaticheck(pickdic, dttolerance, iplot):
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('Wadati-Diagram, %d S-P Times, Vp/Vs(old)=%5.2f, Vp/Vs(checked)=%5.2f' \ plt.title('Wadati-Diagram, %d S-P Times, Vp/Vs(raw)=%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], ['Skipped S-Picks', 'Wadati 1', 'Reliable S-Picks', \
'Wadati 2'], loc='best') 'Wadati 2'], loc='best')