From 8ba34db05c6fbdf4e4d1573464b5040020cf55b1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Wed, 18 Mar 2015 14:44:08 +0100 Subject: [PATCH] New function to calculate earliest and latest possible pick from a given initial (most likely) pick. Replaces class EarlLatePicker of object Picker. --- pylot/core/pick/earllatepicker.py | 138 ++++++++++++++++++++++++++++++ 1 file changed, 138 insertions(+) create mode 100755 pylot/core/pick/earllatepicker.py diff --git a/pylot/core/pick/earllatepicker.py b/pylot/core/pick/earllatepicker.py new file mode 100755 index 00000000..1a2663f7 --- /dev/null +++ b/pylot/core/pick/earllatepicker.py @@ -0,0 +1,138 @@ +#!/usr/bin/python +# -*- coding: utf-8 -*- +""" + Created Mar 2015 + Transcription of the rezipe of Diehl et al. (2009) for consistent phase + picking. For a given inital (the most likely) pick, the corresponding earliest + and latest possible picks are calculated based on noise measurements in front of + the most likely pick and signal wavelength derived from zero crossings. + + :author: MAGS2 EP3 working group / Ludger Kueperkoch +""" +import numpy as np +import matplotlib.pyplot as plt +from obspy.core import Stream +import argparse + +def earllatepicker(X, nfac, TSNR, Pick1, iplot=None): + ''' + Function 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. + + :param: x, time series (seismogram) + :type: `~obspy.core.stream.Stream` + + :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) + + :param: Pick1, initial (prelimenary) onset time, starting point for EarlLatePicker + :type: float + + ''' + + assert isinstance(X, Stream), "%s is not a stream object" % str(X) + + LPick = None + EPick = None + PickError = None + if Pick1 is not None: + print 'earllatepicker: Get earliest and latest possible pick relative to most likely pick ...' + + x =X[0].data + t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate, X[0].stats.delta) + #some parameters needed: + tnoise = TSNR[0] #noise window length for calculating noise level + tsignal = TSNR[2] #signal window length + tsafety = TSNR[1] #safety gap between signal onset and noise window + + #get latest possible pick + #get noise window + inoise = np.where((t <= max([Pick1 - tsafety, 0])) \ + & (t >= max([Pick1 - tnoise - tsafety, 0]))) + #get signal window + isignal = np.where((t <= min([Pick1 + tsignal + tsafety, len(x)])) \ + & (t >= Pick1)) + #calculate noise level + nlevel = max(abs(x[inoise])) * nfac + #get time where signal exceeds nlevel + ilup = np.where(x[isignal] > nlevel) + ildown = np.where(x[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]]) + LPick = t[isignal][il] + + #get earliest possible pick + #get next 2 zero crossings after most likely pick + #if there is one trace in stream + zc = [] + zc.append(Pick1) + i = 0 + for j in range(isignal[0][1],isignal[0][len(t[isignal]) - 1]): + i = i+ 1 + if x[j-1] <= 0 and x[j] >= 0: + zc.append(t[isignal][i]) + elif x[j-1] > 0 and x[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)) + #Ts/4 is assumed as time difference between most likely and earliest possible pick! + EPick = Pick1 - 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 = LPick -Pick1 + diffti_te = Pick1 - EPick + PickError = (diffti_te + 2 * diffti_tl) / 3 + + if iplot is not None: + plt.figure(iplot) + p1, = plt.plot(t, x, 'k') + p2, = plt.plot(t[inoise], x[inoise]) + p3, = plt.plot(t[isignal], x[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([Pick1, Pick1], [max(x), -max(x)], 'b', linewidth=2) + plt.plot([LPick, LPick], [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.xlabel('Time [s] since %s' % X[0].stats.starttime) + plt.yticks([]) + ax = plt.gca() + 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.show() + raw_input() + plt.close(iplot) + + elif Pick1 == None: + print 'earllatepicker: No initial onset time given! Check input!' + return + + return EPick, LPick, PickError + +if __name__ == "__main__": + parser = argparse.ArgumentParser() + parser.add_argument('--X', type=~obspy.core.stream.Stream, help='time series (seismogram) read with obspy module read') + parser.add_argument('--nfac', type=int, help='(noise factor), nfac times noise level to calculate latest possible pick') + parser.add_argument('--TSNR', type=tuple, help='length of time windows around pick used to determine SNR \ + [s] (Tnoise, Tgap, Tsignal)') + parser.add_argument('--Pick1', type=float, help='Onset time of most likely pick') + parser.add_argument('--iplot', type=int, help='if set, figure no. iplot occurs') + args = parser.parse_args() + earllatepicker(args.X, args.nfac, args.TSNR, args.Pick1, args.iplot) + +#earllatepicker(X, nfac, TSNR, Pick1, iplot)