Removed class EarlLatePicker, replaced by new function earllatepicker.

This commit is contained in:
Ludger Küperkoch 2015-03-18 14:45:49 +01:00
parent 16ae4bdfe9
commit 787cac7d68

View File

@ -27,15 +27,11 @@ class AutoPicking(object):
Superclass of different, automated picking algorithms applied on a CF determined
using AIC, HOS, or AR prediction.
'''
def __init__(self, cf, nfac, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, Pick1=None):
def __init__(self, cf, TSNR, PickWindow, iplot=None, aus=None, Tsmooth=None, Pick1=None):
'''
:param: cf, characteristic function, on which the picking algorithm is applied
:type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object
:param: nfac (noise factor), nfac times noise level to calculate latest possible pick
in EarlLatePicker
:type: int
:param: TSNR, length of time windows around pick used to determine SNR [s]
:type: tuple (T_noise, T_gap, T_signal)
@ -63,7 +59,6 @@ class AutoPicking(object):
self.Tcf = cf.getTimeArray()
self.Data = cf.getXCF()
self.dt = cf.getIncrement()
self.setnfac(nfac)
self.setTSNR(TSNR)
self.setPickWindow(PickWindow)
self.setiplot(iplot)
@ -74,25 +69,18 @@ class AutoPicking(object):
def __str__(self):
return '''\n\t{name} object:\n
nfac:\t{nfac}\n
TSNR:\t\t\t{TSNR}\n
PickWindow:\t{PickWindow}\n
aus:\t{aus}\n
Tsmooth:\t{Tsmooth}\n
Pick1:\t{Pick1}\n
'''.format(name=type(self).__name__,
nfac=self.getnfac(),
TSNR=self.getTSNR(),
PickWindow=self.getPickWindow(),
aus=self.getaus(),
Tsmooth=self.getTsmooth(),
Pick1=self.getpick1())
def getnfac(self):
return self.nfac
def setnfac(self, nfac):
self.nfac = nfac
def getTSNR(self):
return self.TSNR
@ -127,15 +115,6 @@ class AutoPicking(object):
def getSlope(self):
return self.slope
def getLpick(self):
return self.LPick
def getEpick(self):
return self.EPick
def getPickError(self):
return self.PickError
def getiplot(self):
return self.iplot
@ -165,7 +144,6 @@ class AICPicker(AutoPicking):
print 'AICPicker: Get initial onset time (pick) from AIC-CF ...'
self.Pick = None
self.PickError = None
#find NaN's
nn = np.isnan(self.cf)
if len(nn) > 1:
@ -263,7 +241,7 @@ class AICPicker(AutoPicking):
p1, = plt.plot(self.Tcf, x / max(x), 'k')
p2, = plt.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r')
if self.Pick is not None:
p3, = plt.plot([self.Pick, self.Pick], [-1 , 1], 'b', linewidth=2)
p3, = plt.plot([self.Pick, self.Pick], [-0.1 , 0.5], 'b', linewidth=2)
plt.legend([p1, p2, p3], ['(HOS-/AR-) Data', 'Smoothed AIC-CF', 'AIC-Pick'])
else:
plt.legend([p1, p2], ['(HOS-/AR-) Data', 'Smoothed AIC-CF'])
@ -281,7 +259,7 @@ class AICPicker(AutoPicking):
p15, = plt.plot(self.Tcf[islope], datafit, 'g', linewidth=2)
plt.legend([p11, p12, p13, p14, p15], ['Data', 'Noise Window', 'Signal Window', 'Slope Window', 'Slope'], \
loc='best')
plt.title('SNR and Slope, Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \
plt.title('Station %s, SNR=%7.2f, Slope= %12.2f counts/s' % (self.Data[0].stats.station, \
self.SNR, self.slope))
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.ylabel('Counts')
@ -307,7 +285,6 @@ class PragPicker(AutoPicking):
print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...'
self.Pick = None
self.PickError = None
self.SNR = None
self.slope = None
#smooth CF
@ -392,313 +369,3 @@ class PragPicker(AutoPicking):
self.Pick = None
print 'PragPicker: No initial onset time given! Check input!'
return
class EarlLatePicker(AutoPicking):
'''
Method to derive earliest and latest possible pick after Diehl & Kissling (2009)
as reasonable uncertainties. Latest possible pick is based on noise level,
earliest possible pick is half a signal wavelength in front of most likely
pick given by PragPicker. Most likely pick (initial pick) must be given.
'''
def calcPick(self):
self.LPick = None
self.EPick = None
self.PickError = None
self.SNR = None
self.slope = None
if self.getpick1() is not None:
print 'EarlLatePicker: Get earliest and latest possible pick relative to most likely pick ...'
ti = self.getpick1()
x = self.Data
t = self.Tcf
#some parameters needed:
tnoise = self.TSNR[0] #noise window length for calculating noise level
tsignal = self.TSNR[2] #signal window length
tsafety = self.TSNR[1] #safety gap between signal onset and noise window
#get latest possible pick
#get noise window
inoise = np.where((self.Tcf <= max([ti - tsafety, 0])) \
& (self.Tcf >= max([ti - tnoise - tsafety, 0])))
#get signal window
isignal = np.where((self.Tcf <= min([ti + tsignal + tsafety, len(x[0].data)])) \
& (self.Tcf >= ti))
#calculate noise level
if len(x) == 1:
nlevel = max(abs(x[0].data[inoise])) * self.nfac
#get time where signal exceeds nlevel
ilup = np.where(x[0].data[isignal] > nlevel)
ildown = np.where(x[0].data[isignal] < -nlevel)
if len(ilup[0]) <= 1 and len(ildown[0]) <= 1:
print 'EarlLatePicker: Signal lower than noise level, misspick?'
return
il = min([ilup[0][0], ildown[0][0]])
self.LPick = t[isignal][il]
elif len(x) == 2:
nlevel = max(np.sqrt(np.power(x[0].data[inoise], 2) + np.power(x[1].data[inoise], 2)))
#get earliest time where signal exceeds nlevel
ilup1 = np.where(x[0].data[isignal] > nlevel)
ilup2 = np.where(x[1].data[isignal] > nlevel)
ildown1 = np.where(x[0].data[isignal] < -nlevel)
ildown2 = np.where(x[1].data[isignal] < -nlevel)
if np.size(ilup1) < 1 and np.size(ilup2) > 1:
ilup = ilup2
elif np.size(ilup1) > 1 and np.size(ilup2) < 1:
ilup = ilup1
elif np.size(ilup1) < 1 and np.size(ilup2) < 1:
ilup = None
else:
ilup = min([ilup1[0][0], ilup2[0][0]])
if np.size(ildown1) < 1 and np.size(ildown2) > 1:
ildown = ildown2
elif np.size(ildown1) > 1 and np.size(ildown2) < 1:
ildown = ildown1
elif np.size(ildown1) < 1 and np.size(ildown2) < 1:
ildown = None
else:
ildown = min([ildown1[0][0], ildown2[0][0]])
if ilup == None and ildown == None:
print 'EarlLatePicker: Signal lower than noise level, misspick?'
return
il = min([ilup, ildown])
self.LPick = t[isignal][il]
elif len(x) == 3:
nlevel = max(np.sqrt(np.power(x[0].data[inoise], 2) + np.power(x[1].data[inoise], 2) + \
np.power(x[2].data[inoise], 2)))
#get earliest time where signal exceeds nlevel
ilup1 = np.where(x[0].data[isignal] > nlevel)
ilup2 = np.where(x[1].data[isignal] > nlevel)
ilup3 = np.where(x[2].data[isignal] > nlevel)
ildown1 = np.where(x[0].data[isignal] < -nlevel)
ildown2 = np.where(x[1].data[isignal] < -nlevel)
ildown3 = np.where(x[2].data[isignal] < -nlevel)
if np.size(ilup1) > 1 and np.size(ilup2) < 1 and np.size(ilup3) < 1:
ilup = ilup1
elif np.size(ilup1) > 1 and np.size(ilup2) > 1 and np.size(ilup3) < 1:
ilup = min([ilup1[0][0], ilup2[0][0]])
elif np.size(ilup1) > 1 and np.size(ilup2) > 1 and np.size(ilup3) > 1:
ilup = min([ilup1[0][0], ilup2[0][0], ilup3[0][0]])
elif np.size(ilup1) < 1 and np.size(ilup2) > 1 and np.size(ilup3) > 1:
ilup = min([ilup2[0][0], ilup3[0][0]])
elif np.size(ilup1) > 1 and np.size(ilup2) < 1 and np.size(ilup3) > 1:
ilup = min([ilup1[0][0], ilup3[0][0]])
elif np.size(ilup1) < 1 and np.size(ilup2) < 1 and np.size(ilup3) < 1:
ilup = None
else:
ilup = min([ilup1[0][0], ilup2[0][0], ilup3[0][0]])
if np.size(ildown1) > 1 and np.size(ildown2) < 1 and np.size(ildown3) < 1:
ildown = ildown1
elif np.size(ildown1) > 1 and np.size(ildown2) > 1 and np.size(ildown3) < 1:
ildown = min([ildown1[0][0], ildown2[0][0]])
elif np.size(ildown1) > 1 and np.size(ildown2) > 1 and np.size(ildown3) > 1:
ildown = min([ildown1[0][0], ildown2[0][0], ildown3[0][0]])
elif np.size(ildown1) < 1 and np.size(ildown2) > 1 and np.size(ildown3) > 1:
ildown = min([ildown2[0][0], ildown3[0][0]])
elif np.size(ildown1) > 1 and np.size(ildown2) < 1 and np.size(ildown3) > 1:
ildown = min([ildown1[0][0], ildown3[0][0]])
elif np.size(ildown1) < 1 and np.size(ildown2) < 1 and np.size(ildown3) < 1:
ildown = None
else:
ildown = min([ildown1[0][0], ildown2[0][0], ildown3[0][0]])
if ilup == None and ildown == None:
print 'EarlLatePicker: Signal lower than noise level, misspick?'
return
il = min([ilup, ildown])
self.LPick = t[isignal][il]
#get earliest possible pick
#get next 2 zero crossings after most likely pick
#if there is one trace in stream
if len(x) == 1:
zc = []
zc.append(ti)
i = 0
for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
i = i+ 1
if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
zc.append(t[isignal][i])
elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
zc.append(t[isignal][i])
if len(zc) == 3:
break
#calculate maximum period of signal out of zero crossings
Ts = max(np.diff(zc))
#if there are two traces in stream
#get maximum of two signal periods
if len(x) == 2:
zc1 = []
zc2 = []
zc1.append(ti)
zc2.append(ti)
i = 0
for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
i = i+ 1
if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
zc1.append(t[isignal][i])
elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
zc1.append(t[isignal][i])
if x[1].data[j-1] <= 0 and x[1].data[j] >= 0:
zc2.append(t[isignal][i])
elif x[1].data[j-1] > 0 and x[1].data[j] <= 0:
zc2.append(t[isignal][i])
if len(zc1) >= 3 and len(zc2) >= 3:
break
Ts = max([max(np.diff(zc1)), max(np.diff(zc2))])
#if there are three traces in stream
#get maximum of three signal periods
if len(x) == 3:
zc1 = []
zc2 = []
zc3 = []
zc1.append(ti)
zc2.append(ti)
zc3.append(ti)
i = 0
for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]):
i = i+ 1
if x[0].data[j-1] <= 0 and x[0].data[j] >= 0:
zc1.append(t[isignal][i])
elif x[0].data[j-1] > 0 and x[0].data[j] <= 0:
zc1.append(t[isignal][i])
if x[1].data[j-1] <= 0 and x[1].data[j] >= 0:
zc2.append(t[isignal][i])
elif x[1].data[j-1] > 0 and x[1].data[j] <= 0:
zc2.append(t[isignal][i])
if x[2].data[j-1] <= 0 and x[2].data[j] >= 0:
zc3.append(t[isignal][i])
elif x[2].data[j-1] > 0 and x[2].data[j] <= 0:
zc3.append(t[isignal][i])
if len(zc1) >= 3 and len(zc2) >= 3 and len(zc3) >= 3:
break
Ts = max([max(np.diff(zc1)), max(np.diff(zc2)), max(np.diff(zc3))])
#Ts/4 is assumed as time difference between most likely and earliest possible pick!
self.EPick = ti - Ts/4
#get symmetric pick error as mean from earliest and latest possible pick
#by weighting latest possible pick tow times earliest possible pick
diffti_tl = self.LPick - ti
diffti_te = ti - self.EPick
self.PickError = (diffti_te + 2 * diffti_tl) / 3
if self.iplot is not None:
plt.figure(self.iplot)
if len(x) == 1:
p1, = plt.plot(t, x[0].data, 'k')
p2, = plt.plot(t[inoise], x[0].data[inoise])
p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
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'], \
loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
elif len(x) == 2:
plt.subplot(2,1,1)
p1, = plt.plot(t, x[0].data, 'k')
p2, = plt.plot(t[inoise], x[0].data[inoise])
p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(zc1[0:3], [0, 0, 0], '*g', markersize=14)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot(zc1[0:3], [0, 0, 0], '*g')
plt.yticks([])
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
plt.subplot(2,1,2)
plt.plot(t, x[1].data, 'k')
plt.plot(t[inoise], x[1].data[inoise])
plt.plot(t[isignal], x[1].data[isignal], 'r')
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([ti, ti], [max(x[1].data), -max(x[1].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot(zc2[0:3], [0, 0, 0], '*g', markersize=14)
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.yticks([])
elif len(x) == 3:
plt.subplot(3,1,1)
p1, = plt.plot(t, x[0].data, 'k')
p2, = plt.plot(t[inoise], x[0].data[inoise])
p3, = plt.plot(t[isignal], x[0].data[isignal], 'r')
p4, = plt.plot([t[0], t[int(len(t)) - 1]], [nlevel, nlevel], '--k')
p5, = plt.plot(zc1[0:3], [0, 0, 0], '*g', markersize=14)
plt.legend([p1, p2, p3, p4, p5], ['Data', 'Noise Window', 'Signal Window', 'Noise Level', 'Zero Crossings'], \
loc='best')
plt.plot([t[0], t[int(len(t)) - 1]], [-nlevel, -nlevel], '--k')
plt.plot([ti, ti], [max(x[0].data), -max(x[0].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[0].data)/2, -max(x[0].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.yticks([])
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.title('Earliest-/Latest Possible/Most Likely Pick & Symmetric Pick Error, %s' % self.Data[0].stats.station)
plt.subplot(3,1,2)
plt.plot(t, x[1].data, 'k')
plt.plot(t[inoise], x[1].data[inoise])
plt.plot(t[isignal], x[1].data[isignal], 'r')
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([ti, ti], [max(x[1].data), -max(x[1].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[1].data)/2, -max(x[1].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot(zc2[0:3], [0, 0, 0], '*g', markersize=14)
plt.yticks([])
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.subplot(3,1,3)
plt.plot(t, x[2].data, 'k')
plt.plot(t[inoise], x[2].data[inoise])
plt.plot(t[isignal], x[2].data[isignal], 'r')
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([ti, ti], [max(x[2].data), -max(x[2].data)], 'b', linewidth=2)
plt.plot([self.LPick, self.LPick], [max(x[2].data)/2, -max(x[2].data)/2], '--k')
plt.plot([self.EPick, self.EPick], [max(x[2].data)/2, -max(x[2].data)/2], '--k')
plt.plot([ti + self.PickError, ti + self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot([ti - self.PickError, ti - self.PickError], [max(x[0].data)/2, -max(x[0].data)/2], 'r--')
plt.plot(zc3[0:3], [0, 0, 0], '*g', markersize=14)
plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
plt.yticks([])
ax = plt.gca()
ax.set_xlim([self.Tcf[inoise[0][0]] - 2, self.Tcf[isignal[0][len(isignal) - 1]] + 3])
plt.show()
raw_input()
plt.close(self.iplot)
elif self.getpick1() == None:
print 'EarlLatePicker: No initial onset time given! Check input!'
return