diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index dc928a85..b3579663 100644 --- a/pylot/core/pick/Picker.py +++ b/pylot/core/pick/Picker.py @@ -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 -