diff --git a/pylot/core/analysis/coinctimes.py b/pylot/core/analysis/coinctimes.py new file mode 100644 index 00000000..b1486898 --- /dev/null +++ b/pylot/core/analysis/coinctimes.py @@ -0,0 +1,57 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from obspy.core import read +from obspy.signal.trigger import coincidenceTrigger + + + +class CoincidenceTimes(): + + def __init__(self, st, comp='Z', coinum=4, sta=1., lta=10., on=5., off=1.): + _type = 'recstalta' + self.coinclist = self.createCoincTriggerlist(data=st, trigcomp=comp, + coinum=coinum, sta=sta, + lta=lta, trigon=on, + trigoff=off, type=_type) + + def __str__(self): + n = 1 + out = '' + for time in self.getCoincTimes(): + out += 'event no. {0}: starttime is {1}\n'.format(n, time) + n += 1 + return out + + def getCoincTimes(self): + timelist = [] + for info in self.getCoincList(): + timelist.append(info['time']) + + return timelist + + def getCoincList(self): + return self.coinclist + + def createCoincTriggerlist(self, data, trigcomp, coinum, sta, lta, + trigon, trigoff, type): + ''' + uses a coincidence trigger to detect all events in the given + dataset + ''' + + triggerlist = coincidenceTrigger(type, trigon, trigoff, + data.select(component=trigcomp), + coinum, sta=sta, lta=lta) + return triggerlist + + +def main(): + data = read('/data/SDS/2014/1A/ZV??/?H?.D/*.365') + data.filter(type='bandpass', freqmin=5., freqmax=30.) + coincs = CoincidenceTimes(data) + print(coincs) + + +if __name__ == '__main__': + main() \ No newline at end of file diff --git a/pylot/core/analysis/trigger.py b/pylot/core/analysis/trigger.py new file mode 100644 index 00000000..add9a687 --- /dev/null +++ b/pylot/core/analysis/trigger.py @@ -0,0 +1,43 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +from obspy.signal.trigger import recSTALTA, triggerOnset + + +def createSingleTriggerlist(st, station='ZV01', trigcomp='Z', stalta=[1, 10], trigonoff=[6, 1]): + ''' + uses a single-station trigger to create a triggerlist for this station + :param st: + :param station: + :param trigcomp: + :param stalta: + :param trigonoff: + :return: + ''' + tr = st.copy().select(component=trigcomp, station=station)[0] + df = tr.stats.sampling_rate + + cft = recSTALTA(tr.data, int(stalta[0] * df), int(stalta[1] * df)) + triggers = triggerOnset(cft, trigonoff[0], trigonoff[1]) + trigg = [] + for time in triggers: + trigg.append(tr.stats.starttime + time[0] / df) + return trigg + + +def createSubCoincTriggerlist(trig, station='ZV01'): + ''' + makes a triggerlist with the events, that are triggered by the + coincidence trigger and are seen at the demanded station + :param trig: list containing triggers from coincidence trigger + :type trig: list + :param station: station name to get triggers for + :type station: str + :return: list of triggertimes + :rtype: list + ''' + trigg = [] + for tri in trig: + if station in tri['stations']: + trigg.append(tri['time']) + return trigg diff --git a/pylot/core/pick/Picker.py b/pylot/core/pick/Picker.py index 1a3bd2be..11e1b6a6 100644 --- a/pylot/core/pick/Picker.py +++ b/pylot/core/pick/Picker.py @@ -33,7 +33,7 @@ class AutoPicking(object): :type: `~pylot.core.pick.CharFuns.CharacteristicFunction` object :param: nfac (noise factor), nfac times noise level to calculate latest possible pick - in EarlLatePick + in EarlLatePicker :type: int :param: TSNR, length of time windows around pick used to determine SNR [s] @@ -121,6 +121,12 @@ class AutoPicking(object): def getpick(self): return self.Pick + def getSNR(self): + return self.SNR + + def getSlope(self): + return self.slope + def getLpick(self): return self.LPick @@ -151,7 +157,7 @@ class AICPicker(AutoPicking): def calcPick(self): - print 'Get onset time (pick) from AIC-CF ...' + print 'AICPicker: Get onset time (pick) from AIC-CF ...' self.Pick = None #find NaN's @@ -187,6 +193,45 @@ class AICPicker(AutoPicking): self.Pick = self.Tcf[i] break + #quality assessment using SNR and slope from CF + if self.Pick is not None: + #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 + tslope = self.TSNR[3] #slope determination window + #get noise window + inoise = np.where((self.Tcf <= max([self.Pick - tsafety, 0])) \ + & (self.Tcf >= max([self.Pick - tnoise - tsafety, 0]))) + #get signal window + isignal = np.where((self.Tcf <= min([self.Pick + tsignal + tsafety, len(self.Data[0].data)])) \ + & (self.Tcf >= self.Pick)) + #calculate SNR from CF + self.SNR = max(abs(self.cf[isignal] - np.mean(self.cf[isignal]))) / max(abs(self.cf[inoise] \ + - np.mean(self.cf[inoise]))) + #calculate slope from CF after initial pick + #get slope window + islope = np.where((self.Tcf <= min([self.Pick + tslope, len(self.Data[0].data)])) \ + & (self.Tcf >= self.Pick)) + #find maximum within slope determination window + #'cause slope should be calculated up to first local minimum only! + imax = np.argmax(self.Data[0].data[islope]) + islope = islope[0][0 :imax] + dataslope = self.Data[0].data[islope] + #calculate slope as polynomal fit of order 1 + xslope = np.arange(0, len(dataslope), 1) + P = np.polyfit(xslope, dataslope, 1) + datafit = np.polyval(P, xslope) + if datafit[0] >= datafit[len(datafit) - 1]: + print 'AICPicker: Negative slope, bad onset skipped!' + return + + self.slope = 1 / tslope * datafit[len(dataslope) - 1] - datafit[0] + + else: + self.SNR = None + self.slope = None + if self.iplot is not None: plt.figure(self.iplot) x = self.Data[0].data @@ -198,14 +243,30 @@ class AICPicker(AutoPicking): plt.yticks([]) plt.title(self.Data[0].stats.station) plt.show() + + if self.Pick is not None: + plt.figure(self.iplot + 1) + p11, = plt.plot(self.Tcf, x, 'k') + p12, = plt.plot(self.Tcf[inoise], self.Data[0].data[inoise]) + p13, = plt.plot(self.Tcf[isignal], self.Data[0].data[isignal], 'r') + p14, = plt.plot(self.Tcf[islope], dataslope, 'g--') + 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, \ + self.SNR, self.slope)) + plt.xlabel('Time [s] since %s' % self.Data[0].stats.starttime) + plt.ylabel('Counts') + ax = plt.gca() + ax.set_ylim([-10, max(self.Data[0].data)]) + ax.set_xlim([self.Tcf[inoise[0][0]] - 5, self.Tcf[isignal[0][len(isignal) - 1]] + 5]) + raw_input() plt.close(self.iplot) if self.Pick == None: print 'AICPicker: Could not find minimum, picking window too short?' - return self.Pick - class PragPicker(AutoPicking): ''' @@ -215,9 +276,11 @@ class PragPicker(AutoPicking): def calcPick(self): if self.getpick1() is not None: - print 'Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...' + print 'PragPicker: Get most likely pick from HOS- or AR-CF using pragmatic picking algorithm ...' self.Pick = None + self.SNR = None + self.slope = None #smooth CF ismooth = int(round(self.Tsmooth / self.dt)) cfsmooth = np.zeros(len(self.cf)) @@ -314,22 +377,26 @@ class EarlLatePicker(AutoPicking): self.LPick = None self.EPick = None + self.SNR = None + self.slope = None if self.getpick1() is not None: - print 'Get earliest and latest possible pick relative to most likely pick ...' + print 'EarlLatePicker: Get earliest and latest possible pick relative to most likely pick ...' ti = self.getpick1() x = self.Data t = self.Tcf - #some parmaters needed: + #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 <= ti - tsafety) & (self.Tcf >= ti - tnoise - tsafety)) + inoise = np.where((self.Tcf <= max([ti - tsafety, 0])) \ + & (self.Tcf >= max([ti - tnoise - tsafety, 0]))) #get signal window - isignal = np.where((self.Tcf <= ti + tsignal) & (self.Tcf >= ti)) + 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 @@ -337,7 +404,7 @@ class EarlLatePicker(AutoPicking): 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 'EarlLatePick: Signal lower than noise level, misspick?' + print 'EarlLatePicker: Signal lower than noise level, misspick?' return il = min([ilup[0][0], ildown[0][0]]) self.LPick = t[isignal][il] @@ -351,7 +418,7 @@ class EarlLatePicker(AutoPicking): ilup = min([ilup1[0][0], ilup2[0][0]]) ildown = min([ildown1[0][0], ildown2[0][0]]) if np.size(ilup) < 1 and np.size(ildown) < 1: - print 'EarlLatePick: Signal lower than noise level, misspick?' + print 'EarlLatePicker: Signal lower than noise level, misspick?' return il = min([ilup, ildown]) self.LPick = t[isignal][il] @@ -368,7 +435,7 @@ class EarlLatePicker(AutoPicking): ilup = min([ilup1[0][0], ilup2[0][0], ilup3[0][0]]) ildown = min([ildown1[0][0], ildown2[0][0], ildown3[0][0]]) if np.size(ilup) < 1 and np.size(ildown) < 1: - print 'EarlLatePick: Signal lower than noise level, misspick?' + print 'EarlLatePicker: Signal lower than noise level, misspick?' return il = min([ilup, ildown]) self.LPick = t[isignal][il] @@ -527,6 +594,6 @@ class EarlLatePicker(AutoPicking): plt.close(self.iplot) elif self.getpick1() == None: - print 'EarlLatePick: No initial onset time given! Check input!' + print 'EarlLatePicker: No initial onset time given! Check input!' return diff --git a/pylot/core/pick/run_makeCF.py b/pylot/core/pick/run_makeCF.py index 41c499b3..e57f80d2 100755 --- a/pylot/core/pick/run_makeCF.py +++ b/pylot/core/pick/run_makeCF.py @@ -2,7 +2,7 @@ # -*- coding: utf-8 -*- """ - Script to run PyLoT modules CharFuns.py and Picker.py. + Script to run autoPyLoT-script "makeCF.py". Only for test purposes! """ @@ -16,18 +16,22 @@ import argparse def run_makeCF(project, database, event, iplot, station=None): #parameters for CF calculation - t2 = 7 #length of moving window for HOS calculation [sec] - p = 4 #order of statistics - cuttimes = [10, 50] #start and end time for CF calculation - bpz = [2, 30] #corner frequencies of bandpass filter, vertical component - bph = [2, 15] #corner frequencies of bandpass filter, horizontal components - tdetz= 1.2 #length of AR-determination window [sec], vertical component - tdeth= 0.8 #length of AR-determination window [sec], horizontal components - tpredz = 0.4 #length of AR-prediction window [sec], vertical component - tpredh = 0.4 #length of AR-prediction window [sec], horizontal components - addnoise = 0.001 #add noise to seismogram for stable AR prediction - arzorder = 2 #chosen order of AR process, vertical component - arhorder = 4 #chosen order of AR process, horizontal components + t2 = 7 #length of moving window for HOS calculation [sec] + p = 4 #order of HOS + cuttimes = [10, 50] #start and end time for CF calculation + bpz = [2, 30] #corner frequencies of bandpass filter, vertical component + bph = [2, 15] #corner frequencies of bandpass filter, horizontal components + tdetz= 1.2 #length of AR-determination window [sec], vertical component + tdeth= 0.8 #length of AR-determination window [sec], horizontal components + tpredz = 0.4 #length of AR-prediction window [sec], vertical component + tpredh = 0.4 #length of AR-prediction window [sec], horizontal components + addnoise = 0.001 #add noise to seismogram for stable AR prediction + arzorder = 2 #chosen order of AR process, vertical component + arhorder = 4 #chosen order of AR process, horizontal components + TSNRhos = [5, 0.5, 1, 0.1] #window lengths [s] for calculating SNR for earliest/latest pick and quality assessment + #from HOS-CF [noise window, safety gap, signal window, slope determination window] + TSNRarz = [5, 0.5, 1, 0.5] #window lengths [s] for calculating SNR for earliest/lates pick and quality assessment + #from ARZ-CF #get waveform data if station: dpz = '/DATA/%s/EVENT_DATA/LOCAL/%s/%s/%s*HZ.msd' % (project, database, event, station) @@ -66,12 +70,12 @@ def run_makeCF(project, database, event, iplot, station=None): aiccf = AICcf(st_copy, cuttimes) #instance of AICcf ############################################################## #get prelimenary onset time from AIC-HOS-CF using subclass AICPicker of class AutoPicking - aicpick = AICPicker(aiccf, None, [5, 0.5, 1], 3, 10, None, 0.1) + aicpick = AICPicker(aiccf, None, TSNRhos, 3, 10, None, 0.1) ############################################################## #get refined onset time from HOS-CF using class Picker - hospick = PragPicker(hoscf, None, [5, 0.5, 1], 2, 10, 0.001, 0.2, aicpick.getpick()) + hospick = PragPicker(hoscf, None, TSNRhos, 2, 10, 0.001, 0.2, aicpick.getpick()) #get earliest and latest possible picks - hosELpick = EarlLatePicker(hoscf, 1.5, [5, 0.5, 1], None, 10, None, None, hospick.getpick()) + hosELpick = EarlLatePicker(hoscf, 1.5, TSNRhos, None, 10, None, None, hospick.getpick()) ############################################################## #calculate ARZ-CF using subclass ARZcf of class CharcteristicFunction #get stream object of filtered data @@ -86,12 +90,12 @@ def run_makeCF(project, database, event, iplot, station=None): araiccf = AICcf(st_copy, cuttimes, tpredz, 0, tdetz) #instance of AICcf ############################################################## #get onset time from AIC-ARZ-CF using subclass AICPicker of class AutoPicking - aicarzpick = AICPicker(araiccf, 1.5, [5, 0.5, 2], 2, 10, None, 0.1) + aicarzpick = AICPicker(araiccf, 1.5, TSNRarz, 2, 10, None, 0.1) ############################################################## #get refined onset time from ARZ-CF using class Picker - arzpick = PragPicker(arzcf, 1.5, [5, 0.5, 2], 2.0, 10, 0.1, 0.05, aicarzpick.getpick()) + arzpick = PragPicker(arzcf, 1.5, TSNRarz, 2.0, 10, 0.1, 0.05, aicarzpick.getpick()) #get earliest and latest possible picks - arzELpick = EarlLatePicker(arzcf, 1.5, [5, 0.5, 1], None, 10, None, None, arzpick.getpick()) + arzELpick = EarlLatePicker(arzcf, 1.5, TSNRarz, None, 10, None, None, arzpick.getpick()) elif not wfzfiles: print 'No vertical component data found!' @@ -127,12 +131,12 @@ def run_makeCF(project, database, event, iplot, station=None): arhaiccf = AICcf(H_copy, cuttimes, tpredh, 0, tdeth) #instance of AICcf ############################################################## #get onset time from AIC-ARH-CF using subclass AICPicker of class AutoPicking - aicarhpick = AICPicker(arhaiccf, 1.5, [5, 0.5, 2], 4, 10, None, 0.1) + aicarhpick = AICPicker(arhaiccf, 1.5, TSNRarz, 4, 10, None, 0.1) ############################################################### #get refined onset time from ARH-CF using class Picker - arhpick = PragPicker(arhcf, 1.5, [5, 0.5, 2], 2.5, 10, 0.1, 0.05, aicarhpick.getpick()) + arhpick = PragPicker(arhcf, 1.5, TSNRarz, 2.5, 10, 0.1, 0.05, aicarhpick.getpick()) #get earliest and latest possible picks - arhELpick = EarlLatePicker(arhcf, 1.5, [5, 0.5, 1], None, 10, None, None, arhpick.getpick()) + arhELpick = EarlLatePicker(arhcf, 1.5, TSNRarz, None, 10, None, None, arhpick.getpick()) #create stream with 3 traces #merge streams @@ -155,7 +159,7 @@ def run_makeCF(project, database, event, iplot, station=None): #calculate AR3C-CF using subclass AR3Ccf of class CharacteristicFunction ar3ccf = AR3Ccf(AllC, cuttimes, tpredz, arhorder, tdetz, addnoise) #instance of AR3Ccf #get earliest and latest possible pick from initial ARH-pick - ar3cELpick = EarlLatePicker(ar3ccf, 1.5, [5, 0.5, 1], None, 10, None, None, arhpick.getpick()) + ar3cELpick = EarlLatePicker(ar3ccf, 1.5, TSNRarz, None, 10, None, None, arhpick.getpick()) ############################################################## if iplot: #plot vertical trace @@ -187,7 +191,8 @@ def run_makeCF(project, database, event, iplot, station=None): plt.ylim([-1.5, 1.5]) plt.xlabel('Time [s]') plt.ylabel('Normalized Counts') - plt.title([tr.stats.station, tr.stats.channel]) + plt.title('%s, %s, AIC-SNR=%7.2f, AIC-Slope=%12.2f' % (tr.stats.station, \ + tr.stats.channel, aicpick.getSNR(), aicpick.getSlope())) plt.suptitle(tr.stats.starttime) plt.legend([p1, p2, p3, p4, p5], ['Data', 'HOS-CF', 'HOSAIC-CF', 'ARZ-CF', 'ARZAIC-CF']) #plot horizontal traces