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

autoPyLoT.py

@@ -13,6 +13,7 @@ from pylot.core.read import Data, AutoPickParameter
 from pylot.core.pick.run_autopicking import run_autopicking
 from pylot.core.util.structure import DATASTRUCTURE
 from pylot.core.pick.utils import wadaticheck
+import pdb
 __version__ = _getVersionString()
 
 
@@ -30,6 +31,18 @@ def autoPyLoT(inputfile):
     .. rubric:: Example
 
     '''
+    print '************************************'
+    print '*********autoPyLoT starting*********'
+    print 'The Python picking and Location Tool'
+    print '      Version ', _getVersionString(), '2015'
+    print '**Authors:'
+    print '**S. Wehling-Benatelli'
+    print '** Ruhr-University Bochum'
+    print '**L. Kueperkoch'
+    print '** BESTEC GmbH'
+    print '**K. Olbert'
+    print '** Christian-Albrechts University Kiel'
+    print '************************************'
 
     # reading parameter file
 
@@ -115,7 +128,6 @@ def autoPyLoT(inputfile):
             station = wfdat[0].stats.station
             allonsets = {station: picks}
             for i in range(len(wfdat)):
-            #for i in range(0,5):
                 stationID = wfdat[i].stats.station
                 #check if station has already been processed
                 if stationID not in procstats:
@@ -139,6 +151,10 @@ def autoPyLoT(inputfile):
             print '-------Finished event %s!-------' % parameter.getParam('eventID')
             print '------------------------------------------'
             
+    print '************************************'
+    print '*********autoPyLoT terminates*******'
+    print 'The Python picking and Location Tool'
+    print '************************************'
 
 if __name__ == "__main__":
     # parse arguments

autoPyLoT_local.in

@@ -1,6 +1,7 @@
 %This is a parameter input file for autoPyLoT.
 %All main and special settings regarding data handling
 %and picking are to be set here!
+%Parameters are optimized for local data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 #main settings#

autoPyLoT_regional.in

@@ -1,13 +1,14 @@
 %This is a parameter input file for autoPyLoT.
 %All main and special settings regarding data handling
 %and picking are to be set here!
+%Parameters are optimized for regional data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 #main settings#
 /DATA/Egelados                      #rootpath#     %project path
 EVENT_DATA/LOCAL		    #datapath#     %data path
-2006.02_Nisyros                     #database#     %name of data base
+2006.01_Nisyros                     #database#     %name of data base
-e0032.033.06                        #eventID#      %event ID for single event processing
+e1412.008.06                        #eventID#      %event ID for single event processing
 PILOT				    #datastructure#  %choose data structure
 2                                   #iplot#        %flag for plotting: 0 none, 1, partly, >1 everything
 AUTOPHASES_AIC_HOS4_ARH             #phasefile#    %name of autoPILOT output phase file 
@@ -30,8 +31,8 @@ HYPOSAT                             #locrt#        %location routine used ("HYPO
 #common settings picker#
 20                                  #pstart#       %start time [s] for calculating CF for P-picking
 160                                 #pstop#        %end time [s] for calculating CF for P-picking
-1.0                                 #sstart#       %start time [s] after or before(-) P-onset for calculating CF for S-picking
+3.0                                 #sstart#       %start time [s] after or before(-) P-onset for calculating CF for S-picking
-50                                  #sstop#        %end time [s] after P-onset for calculating CF for S-picking
+100                                 #sstop#        %end time [s] after P-onset for calculating CF for S-picking
 3 10                                #bpz1#         %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
 3 12                                #bpz2#         %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
 3 8                                 #bph1#         %lower/upper corner freq. of first band pass filter H-comp. [Hz]
@@ -48,13 +49,12 @@ HOS                                 #algoP#        %choose algorithm for P-onset
 0.6                                 #tdet2z#       %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick
 0.2                                 #tpred2z#      %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick
 0.001                               #addnoise#     %add noise to seismogram for stable AR prediction
-4 0.1 1.0 0.5                       #tsnrz#        %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
+4 0.2 2.0 1.5                       #tsnrz#        %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
-4                                   #pickwinP#     %for initial AIC pick, length of P-pick window [s]
+4                                   #pickwinP#     %for initial AIC and refined pick, length of P-pick window [s]
 8                                   #Precalcwin#   %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
-0                                   #peps4aic#     %for HOS/AR, artificial uplift of samples of AIC-function (P)
+3.0                                 #aictsmooth#   %for HOS/AR, take average of samples for smoothing of AIC-function [s]
-0.2                                 #aictsmooth#   %for HOS/AR, take average of samples for smoothing of AIC-function [s]
+0.3                                 #tsmoothP#     %for HOS/AR, take average of samples for smoothing CF [s]
-0.1                                 #tsmoothP#     %for HOS/AR, take average of samples for smoothing CF [s]
+0.3                                 #ausP#         %for HOS/AR, artificial uplift of samples (aus) of CF (P)
-0.001                               #ausP#         %for HOS/AR, artificial uplift of samples (aus) of CF (P)
 1.3                                 #nfacP#        %for HOS/AR, noise factor for noise level determination (P)
 #H-components#
 ARH                                 #algoS#        %choose algorithm for S-onset determination (ARH or AR3)
@@ -63,18 +63,17 @@ ARH                                 #algoS#        %choose algorithm for S-onset
 0.6                                 #tdet2h#       %for HOS/AR, length of AR-determinaton window [s], H-components, 2nd pick
 0.3                                 #tpred2h#      %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick
 4                                   #Sarorder#     %for AR-picker, order of AR process of H-components
-20                                  #Srecalcwin#   %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
+10                                  #Srecalcwin#   %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
-5                                   #pickwinS#     %for initial AIC pick, length of S-pick window [s]
+6                                   #pickwinS#     %for initial AIC and refined pick, length of S-pick window [s]
-6 0.2 2.0 1.5                       #tsnrh#        %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
+5 0.2 3.0 3.0                       #tsnrh#        %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
-0.05                                #aictsmoothS#  %for AIC-picker, take average of samples for smoothing of AIC-function [s]
+3.0                                 #aictsmoothS#  %for AIC-picker, take average of samples for smoothing of AIC-function [s]
-0.02                                #tsmoothS#     %for AR-picker, take average of samples for smoothing CF [s] (S)
+1.0                                 #tsmoothS#     %for AR-picker, take average of samples for smoothing CF [s] (S)
-0.2                                 #pepsS#        %for AR-picker, artificial uplift of samples of CF (S) 
 0.4                                 #ausS#         %for HOS/AR, artificial uplift of samples (aus) of CF (S)
 1.5                                 #nfacS#        %for AR-picker, noise factor for noise level determination (S)
 %first-motion picker%
 1                                   #minfmweight#  %minimum required p weight for first-motion determination
 2                                   #minFMSNR#     %miniumum required SNR for first-motion determination
-5.0                                 #fmpickwin#    %pick window around P onset for calculating zero crossings
+6.0                                 #fmpickwin#    %pick window around P onset for calculating zero crossings
 %quality assessment%
 #inital AIC onset#
 0.04 0.08 0.16 0.32                 #timeerrorsP#   %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P

pylot/core/pick/run_autopicking.py

@@ -74,6 +74,10 @@ def run_autopicking(wfstream, pickparam):
     minFMSNR = pickparam.getParam('minFMSNR')
     fmpickwin = pickparam.getParam('fmpickwin')
     minfmweight = pickparam.getParam('minfmweight')
+    # parameters for checking signal length
+    minsiglength = pickparam.getParam('minsiglength')
+    minpercent = pickparam.getParam('minpercent')
+    nfacsl = pickparam.getParam('noisefactor')
 
     # initialize output 
     Pweight = 4   # weight for P onset 
@@ -94,6 +98,7 @@ def run_autopicking(wfstream, pickparam):
 
     aicSflag = 0
     aicPflag = 0
+    Pflag = 0
     Sflag = 0
 
     # split components
@@ -152,9 +157,15 @@ def run_autopicking(wfstream, pickparam):
         # of class AutoPicking
         aicpick = AICPicker(aiccf, tsnrz, pickwinP, iplot, None, tsmoothP)
         ##############################################################
+        # check signal length to detect spuriously picked noise peaks
+        z_copy[0].data = tr_filt.data
+        Pflag = checksignallength(z_copy, aicpick.getpick(), tsnrz, minsiglength, \
+                                          nfacsl, minpercent, iplot)
+        ##############################################################
         # go on with processing if AIC onset passes quality control
         if (aicpick.getSlope() >= minAICPslope and
-                aicpick.getSNR() >= minAICPSNR):
+                aicpick.getSNR() >= minAICPSNR and
+                Pflag == 1):
             aicPflag = 1
             print 'AIC P-pick passes quality control: Slope: %f, SNR: %f' % \
                   (aicpick.getSlope(), aicpick.getSNR())
@@ -190,7 +201,7 @@ def run_autopicking(wfstream, pickparam):
             mpickP = refPpick.getpick()
             #############################################################
             if mpickP is not None:
-            	# quality assessment
+           	# quality assessment
             	# get earliest and latest possible pick and symmetrized uncertainty
             	[lpickP, epickP, Perror] = earllatepicker(z_copy, nfacP, tsnrz, mpickP, iplot)
 
@@ -227,8 +238,8 @@ def run_autopicking(wfstream, pickparam):
             Sflag = 0
 
     else:
-        print 'run_autopicking: No vertical component data available, ' \
+        print 'run_autopicking: No vertical component data availabler!, ' \
-              'skipping station!'
+              'Skip station!'
 
     if edat is not None and ndat is not None and len(edat) > 0 and len(
             ndat) > 0 and Pweight < 4:

pylot/core/pick/utils.py

@@ -9,6 +9,7 @@
 """
 
 import numpy as np
+import scipy as sc
 import matplotlib.pyplot as plt
 from obspy.core import Stream, UTCDateTime
 import warnings
@@ -47,14 +48,11 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None):
     x = X[0].data
     t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
                   X[0].stats.delta)
-    # get latest possible pick
-    # get noise window
     inoise = getnoisewin(t, Pick1, TSNR[0], TSNR[1])
     # get signal window
     isignal = getsignalwin(t, Pick1, TSNR[2])
     # remove mean
-    meanwin = np.hstack((inoise, isignal))
+    x = x - np.mean(x[inoise])
-    x = x - np.mean(x[meanwin])
     # calculate noise level
     nlevel = np.sqrt(np.mean(np.square(x[inoise]))) * nfac
     # get time where signal exceeds nlevel
@@ -337,7 +335,7 @@ def getSNR(X, TSNR, t1):
         return
 
     # demean over entire snr window
-    x -= x[inoise[0]:isignal[-1]].mean()
+    x = x - np.mean(x[np.hstack([inoise, isignal])])
 
     # calculate ratios
     noiselevel = np.sqrt(np.mean(np.square(x[inoise])))
@@ -514,3 +512,88 @@ def wadaticheck(pickdic, dttolerance, iplot):
         plt.close(iplot)
 
     return checkedonsets
+
+
+def checksignallength(X, pick, TSNR, minsiglength, nfac, minpercent, iplot):
+    '''
+    Function to detect spuriously picked noise peaks.
+    Uses envelope to determine, how many samples [per cent] after
+    P onset are below certain threshold, calculated from noise 
+    level times noise factor.
+
+    : param: X, time series (seismogram)
+    : type:  `~obspy.core.stream.Stream`
+
+    : param: pick, initial (AIC) P onset time
+    : type:  float 
+
+    : param: TSNR, length of time windows around initial pick [s]
+    : type:  tuple (T_noise, T_gap, T_signal)
+
+    : param: minsiglength, minium required signal length [s] to
+             declare pick as P onset
+    : type:  float
+
+    : param: nfac, noise factor (nfac * noise level = threshold)
+    : type:  float
+
+    : param: minpercent, minimum required percentage of samples
+             above calculated threshold
+    : type:  float
+
+    : param: iplot, if iplot > 1, results are shown in figure
+    : type:  int
+    '''
+    
+    assert isinstance(X, Stream), "%s is not a stream object" % str(X)
+
+    print 'Checking signal length ...'
+
+    x = X[0].data
+    t = np.arange(0, X[0].stats.npts / X[0].stats.sampling_rate,
+                  X[0].stats.delta)
+
+    # generate envelope function from Hilbert transform
+    y = np.imag(sc.signal.hilbert(x))
+    e = np.sqrt(np.power(x, 2) + np.power(y, 2))
+    # get noise window
+    inoise = getnoisewin(t, pick, TSNR[0], TSNR[1])
+    # get signal window
+    isignal = getsignalwin(t, pick, TSNR[2])
+    # calculate minimum adjusted signal level
+    minsiglevel = max(e[inoise]) * nfac
+    # minimum adjusted number of samples over minimum signal level
+    minnum = len(isignal) * minpercent/100
+    # get number of samples above minimum adjusted signal level
+    numoverthr = len(np.where(e[isignal] >= minsiglevel)[0]) 
+
+    if numoverthr >= minnum:
+    	print 'checksignallength: Signal reached required length.'
+        returnflag = 1
+    else:
+        print 'checksignallength: Signal shorter than required minimum signal length!'
+        print 'Presumably picked picked noise peak, pick is rejected!'
+        returnflag = 0
+ 
+    if iplot == 2:
+        plt.figure(iplot)
+    	p1, = plt.plot(t,x, 'k')
+        p2, = plt.plot(t[inoise], e[inoise])
+        p3, = plt.plot(t[isignal],e[isignal], 'r') 
+        p4, = plt.plot([t[isignal[0]], t[isignal[len(isignal)-1]]], \
+                        [minsiglevel, minsiglevel], 'g')
+        p5, = plt.plot([pick, pick], [min(x), max(x)], 'c')
+        plt.legend([p1, p2, p3, p4, p5], ['Data', 'Envelope Noise Window', \
+                    'Envelope Signal Window', 'Minimum Signal Level', \
+                    'Onset'], loc='best')
+        plt.xlabel('Time [s] since %s' % X[0].stats.starttime)
+        plt.ylabel('Counts')
+        plt.title('Check for Signal Length, Station %s' % X[0].stats.station)
+        plt.yticks([])
+        plt.show()
+        raw_input()
+        plt.close(iplot)
+
+    return returnflag
+
+