Cleaned up source code, debuged: calculates now T/4 instead of T/8 out of zero crossings.

This commit is contained in:
Ludger Küperkoch 2015-03-30 14:18:19 +02:00
parent 0ad7e629ce
commit bebe3a3c45

View File

@ -43,86 +43,81 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
LPick = None LPick = None
EPick = None EPick = None
PickError = None PickError = None
if Pick1 is not 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, X[0].stats.delta) t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta)
#some parameters needed: #some parameters needed:
tnoise = TSNR[0] #noise window length for calculating noise level tnoise = TSNR[0] #noise window length for calculating noise level
tsignal = TSNR[2] #signal window length tsignal = TSNR[2] #signal window length
tsafety = TSNR[1] #safety gap between signal onset and noise window tsafety = TSNR[1] #safety gap between signal onset and noise window
#get latest possible pick #get latest possible pick
#get noise window #get noise window
inoise = np.where((t <= max([Pick1 - tsafety, 0])) \ inoise = np.where((t <= max([Pick1 - tsafety, 0])) \
& (t >= max([Pick1 - tnoise - tsafety, 0]))) & (t >= max([Pick1 - tnoise - tsafety, 0])))
#get signal window #get signal window
isignal = np.where((t <= min([Pick1 + tsignal + tsafety, len(x)])) \ isignal = np.where((t <= min([Pick1 + tsignal + tsafety, len(x)])) \
& (t >= Pick1)) & (t >= Pick1))
#calculate noise level #calculate noise level
nlevel = max(abs(x[inoise])) * nfac nlevel = max(abs(x[inoise])) * nfac
#get time where signal exceeds nlevel #get time where signal exceeds nlevel
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 len(ilup[0]) <= 1 and len(ildown[0]) <= 1: if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
print 'earllatepicker: Signal lower than noise level, misspick?' print 'earllatepicker: Signal lower than noise level, misspick?'
return return
il = min([ilup[0][0], ildown[0][0]]) il = min([ilup[0][0], ildown[0][0]])
LPick = t[isignal][il] LPick = t[isignal][il]
#get earliest possible pick #get earliest possible pick
#get next 2 zero crossings after most likely pick #get next 2 zero crossings after most likely pick
#initial onset is assumed to be the first zero crossing #initial onset is assumed to be the first zero crossing
zc = [] zc = []
zc.append(Pick1) zc.append(Pick1)
i = 0 i = 0
for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]): for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
i = i+ 1 i = i+ 1
if x[j-1] <= 0 and x[j] >= 0: if x[j-1] <= 0 and x[j] >= 0:
zc.append(t[isignal][i]) zc.append(t[isignal][i])
elif x[j-1] > 0 and x[j] <= 0: elif x[j-1] > 0 and x[j] <= 0:
zc.append(t[isignal][i]) zc.append(t[isignal][i])
if len(zc) == 3: if len(zc) == 3:
break break
#calculate maximum period of signal out of zero crossings #calculate maximum period T0 of signal out of zero crossings
Ts = max(np.diff(zc)) T0 = max(np.diff(zc)) #this is half wave length!
#Ts/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!
EPick = Pick1 - Ts/4 EPick = Pick1 - T0/2
#get symmetric pick error as mean from earliest and latest possible pick #get symmetric pick error as mean from earliest and latest possible pick
#by weighting latest possible pick two times earliest possible pick #by weighting latest possible pick two times earliest possible pick
diffti_tl = LPick -Pick1 diffti_tl = LPick - Pick1
diffti_te = Pick1 - EPick diffti_te = Pick1 - EPick
PickError = (diffti_te + 2 * diffti_tl) / 3 PickError = (diffti_te + 2 * diffti_tl) / 3
if iplot is not None: if iplot is not None:
plt.figure(iplot) plt.figure(iplot)
p1, = plt.plot(t, x, 'k') p1, = plt.plot(t, x, 'k')
p2, = plt.plot(t[inoise], x[inoise]) p2, = plt.plot(t[inoise], x[inoise])
p3, = plt.plot(t[isignal], x[isignal], 'r') p3, = plt.plot(t[isignal], x[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k') p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14) p5, = plt.plot(zc, [0, 0, 0], '*g', markersize=14)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \ plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
loc='best') loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k') plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2) plt.plot([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2)
plt.plot([LPick, LPick], [max(x)/2, -max(x)/2], '--k') plt.plot([LPick, LPick], [max(x)/2, -max(x)/2], '--k')
plt.plot([EPick, EPick], [max(x)/2, -max(x)/2], '--k') plt.plot([EPick, EPick], [max(x)/2, -max(x)/2], '--k')
plt.plot([Pick1 + PickError, Pick1 + PickError], [max(x)/2, -max(x)/2], 'r--') plt.plot([Pick1 + PickError, Pick1 + PickError], [max(x)/2, -max(x)/2], 'r--')
plt.plot([Pick1 - PickError, Pick1 - PickError], [max(x)/2, -max(x)/2], 'r--') plt.plot([Pick1 - PickError, Pick1 - PickError], [max(x)/2, -max(x)/2], 'r--')
plt.xlabel('Time [s] since %s' % X[0].stats.starttime) plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
plt.yticks([]) plt.yticks([])
ax = plt.gca() ax = plt.gca()
ax.set_xlim([t[inoise[0][0]] - 2, t[isignal[0][len(isignal) - 1]] + 3]) ax.set_xlim([t[inoise[0][0]] - 2, t[isignal[0][len(isignal) - 1]] + 3])
plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % X[0].stats.station) plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % X[0].stats.station)
plt.show() plt.show()
raw_input() raw_input()
plt.close(iplot) plt.close(iplot)
elif Pick1 == None:
print 'earllatepicker: No initial onset time given! Check input!'
return
return EPick, LPick, PickError return EPick, LPick, PickError