Merge branch 'develop' of ariadne:/data/git/pylot into develop

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()

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

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
 

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]
+    t2 = 7                   #length of moving window for HOS calculation [sec]
-    p = 4               #order of statistics
+    p = 4                    #order of HOS
-    cuttimes = [10, 50] #start and end time for CF calculation
+    cuttimes = [10, 50]      #start and end time for CF calculation
-    bpz = [2, 30]       #corner frequencies of bandpass filter, vertical component
+    bpz = [2, 30]            #corner frequencies of bandpass filter, vertical component
-    bph = [2, 15]       #corner frequencies of bandpass filter, horizontal components
+    bph = [2, 15]            #corner frequencies of bandpass filter, horizontal components
-    tdetz= 1.2          #length of AR-determination window [sec], vertical component
+    tdetz= 1.2               #length of AR-determination window [sec], vertical component
-    tdeth= 0.8          #length of AR-determination window [sec], horizontal components
+    tdeth= 0.8               #length of AR-determination window [sec], horizontal components
-    tpredz = 0.4        #length of AR-prediction window [sec], vertical component
+    tpredz = 0.4             #length of AR-prediction window [sec], vertical component
-    tpredh = 0.4        #length of AR-prediction window [sec], horizontal components
+    tpredh = 0.4             #length of AR-prediction window [sec], horizontal components
-    addnoise = 0.001    #add noise to seismogram for stable AR prediction
+    addnoise = 0.001         #add noise to seismogram for stable AR prediction
-    arzorder = 2        #chosen order of AR process, vertical component
+    arzorder = 2             #chosen order of AR process, vertical component
-    arhorder = 4        #chosen order of AR process, horizontal components
+    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