Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

autoPyLoT.py

@@ -6,11 +6,11 @@ import argparse
 import glob
 import subprocess
 import string
-from obspy.core import read
+from obspy.core import read, UTCDateTime
 from pylot.core.read.data import Data
 from pylot.core.read.inputs import AutoPickParameter
 from pylot.core.util.structure import DATASTRUCTURE
-from pylot.core.pick.autopick import autopickevent
+from pylot.core.pick.autopick import autopickevent, iteratepicker
 from pylot.core.loc.nll import *
 from pylot.core.util.version import get_git_version as _getVersionString
 
@@ -121,7 +121,28 @@ def autoPyLoT(inputfile):
 
                     # !iterative picking if traces remained unpicked or occupied with bad picks!
                     # get theoretical onset times for picks with weights >= 4
-                    # in order to reprocess them using smaller time windows
+                    # in order to reprocess them using smaller time windows around theoretical onset
+                    # get stations with bad onsets
+                    badpicks = []
+                    for key in picks:
+                         if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
+                             badpicks.append([key, picks[key]['P']['mpp']])
+
+                    if len(badpicks) == 0:
+                        print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") 
+                    else:
+                        # get theoretical P-onset times from NLLoc-location file
+                        locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
+                        # get latest file if several are available
+                        nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                        if os.path.isfile(nllocfile):
+                            picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
+                            # write phases to NLLoc-phase file
+                            picksExport(picks, 'NLLoc', phasefile)
+                            # locate the event
+                            locate(nlloccall, ctrfile)
+                        else:
+                            print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
                 ##########################################################
                 # write phase files for various location routines
                 # HYPO71
@@ -160,10 +181,30 @@ def autoPyLoT(inputfile):
 
                 # locate the event
                 locate(nlloccall, ctrfile)
-
                 # !iterative picking if traces remained unpicked or occupied with bad picks!
                 # get theoretical onset times for picks with weights >= 4
-                # in order to reprocess them using smaller time windows
+                # in order to reprocess them using smaller time windows around theoretical onset
+                # get stations with bad onsets
+                badpicks = []
+                for key in picks:
+                     if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
+                         badpicks.append([key, picks[key]['P']['mpp']])
+
+                if len(badpicks) == 0:
+                    print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") 
+                else:
+                    # get theoretical P-onset times from NLLoc-location file
+                    locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
+                    # get latest file if several are available
+                    nllocfile = max(glob.glob(locsearch), key=os.path.getctime) 
+                    if os.path.isfile(nllocfile):
+                        picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
+                        # write phases to NLLoc-phase file
+                        picksExport(picks, 'NLLoc', phasefile)
+                        # locate the event
+                        locate(nlloccall, ctrfile)
+                    else:
+                        print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
             ##########################################################
             # write phase files for various location routines
             # HYPO71

autoPyLoT_local.in

@@ -17,7 +17,7 @@ PILOT				    #datastructure#%choose data structure
 /home/ludger/NLLOC/Insheim          #nllocroot#    %root of NLLoc-processing directory
 AUTOPHASES.obs                      #phasefile#    %name of autoPyLoT-output phase file for NLLoc
                                                    %(in nllocroot/obs)
-Insheim_min1d2015_auto.in           #locfile#      %name of autoPyLoT-output control file for NLLoc
+Insheim_min1d2015_auto.in           #ctrfile#      %name of autoPyLoT-output control file for NLLoc
                                                    %(in nllocroot/run)
 ttime                               #ttpatter#     %pattern of NLLoc ttimes from grid
                                                    %(in nllocroot/times)
@@ -41,7 +41,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in          #focmecin#     %name of focmec input file co
 #common settings picker#
 15                                  #pstart#       %start time [s] for calculating CF for P-picking
 60                                  #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
+-1.0                                #sstart#       %start time [s] relative to P-onset for calculating CF for S-picking
 7                                   #sstop#        %end time [s] after P-onset for calculating CF for S-picking
 2 20                                #bpz1#         %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
 2 30                                #bpz2#         %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
@@ -51,7 +51,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in          #focmecin#     %name of focmec input file co
 %!!Be careful when editing the following!!
 #Z-component#
 HOS                                 #algoP#        %choose algorithm for P-onset determination (HOS, ARZ, or AR3)
-7                                   #tlta#         %for HOS-/AR-AIC-picker, length of LTA window [s] 
+7.0                                 #tlta#         %for HOS-/AR-AIC-picker, length of LTA window [s] 
 4                                   #hosorder#     %for HOS-picker, order of Higher Order Statistics
 2                                   #Parorder#     %for AR-picker, order of AR process of Z-component
 1.2                                 #tdet1z#       %for AR-picker, length of AR determination window [s] for Z-component, 1st pick
@@ -60,8 +60,8 @@ HOS                                 #algoP#        %choose algorithm for P-onset
 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
 3 0.1 0.5 0.5                       #tsnrz#        %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
-3                                   #pickwinP#     %for initial AIC pick, length of P-pick window [s]
+3.0                                 #pickwinP#     %for initial AIC pick, length of P-pick window [s]
-8                                   #Precalcwin#   %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
+8.0                                 #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)
 0.2                                 #aictsmooth#   %for HOS/AR, take average of samples for smoothing of AIC-function [s]
 0.1                                 #tsmoothP#     %for HOS/AR, take average of samples for smoothing CF [s]
@@ -74,8 +74,8 @@ 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
-6                                   #Srecalcwin#   %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
+6.0                                 #Srecalcwin#   %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
-3                                   #pickwinS#     %for initial AIC pick, length of S-pick window [s]
+3.0                                 #pickwinS#     %for initial AIC pick, length of S-pick window [s]
 2 0.2 1.5 0.5                       #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]
 0.02                                #tsmoothS#     %for AR-picker, take average of samples for smoothing CF [s] (S)
@@ -83,7 +83,7 @@ ARH                                 #algoS#        %choose algorithm for S-onset
 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
+1                                   #minfmweight#  %minimum required P weight for first-motion determination
 2                                   #minFMSNR#     %miniumum required SNR for first-motion determination
 0.2                                 #fmpickwin#    %pick window around P onset for calculating zero crossings
 %quality assessment%
@@ -92,7 +92,7 @@ ARH                                 #algoS#        %choose algorithm for S-onset
 0.04 0.08 0.16 0.32                 #timeerrorsS#   %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S
 10                                  #minAICPslope#  %below this slope [counts/s] the initial P pick is rejected
 1.2                                 #minAICPSNR#    %below this SNR the initial P pick is rejected
-6                                   #minAICSslope#  %below this slope [counts/s] the initial S pick is rejected
+3                                   #minAICSslope#  %below this slope [counts/s] the initial S pick is rejected
 1.5                                 #minAICSSNR#    %below this SNR the initial S pick is rejected
 #check duration of signal using envelope function#
 5                                   #minsiglength# %minimum required length of signal [s]

pylot/core/pick/autopick.py

@@ -3,7 +3,7 @@
 
 """
 Function to run automated picking algorithms using AIC,
-HOS and AR prediction. Uses object CharFuns and Picker and
+HOS and AR prediction. Uses objects CharFuns and Picker and
 function conglomerate utils.
 
 :author: MAGS2 EP3 working group / Ludger Kueperkoch
@@ -16,6 +16,7 @@ from pylot.core.pick.Picker import AICPicker, PragPicker
 from pylot.core.pick.CharFuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
 from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker,\
     getSNR, fmpicker, checkPonsets, wadaticheck, crossings_nonzero_all
+from pylot.core.util.utils import getPatternLine
 from pylot.core.read.data import Data
 from pylot.core.analysis.magnitude import WApp, DCfc
 
@@ -316,12 +317,11 @@ def autopickstation(wfstream, pickparam):
                 # from P pulse
                 # initialize Data object
                 data = Data()
-                [corzdat, restflag] = data.restituteWFData(invdir, zdat)
+                z_copy = zdat.copy()
+                [corzdat, restflag] = data.restituteWFData(invdir, z_copy)
                 if restflag == 1:
                     # integrate to displacement
                     corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta)
-                    # class needs stream object => build it
-                    z_copy = zdat.copy()
                     z_copy[0].data = corintzdat
                     # largest detectable period == window length
                     # after P pulse for calculating source spectrum
@@ -799,3 +799,77 @@ def autopickstation(wfstream, pickparam):
     picks[phase]['Ao'] = Ao
 
     return picks
+
+def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
+    '''
+    Repicking of bad onsets. Uses theoretical onset times from NLLoc-location file.
+
+    :param wf: waveform, obspy stream object
+
+    :param NLLocfile: path/name of NLLoc-location file
+
+    :param picks: dictionary of available onset times
+
+    :param badpicks: picks to be repicked
+
+    :param pickparameter: picking parameters from autoPyLoT-input file 
+    '''
+
+    print("#######################################################")
+    print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \
+                               % badpicks
+
+    newpicks = {}
+    for i in range(0, len(badpicks)):
+         if len(badpicks[i][0]) > 4:
+             Ppattern = '%s  ?    ?    ? P' % badpicks[i][0]
+         elif len(badpicks[i][0]) == 4:
+             Ppattern = '%s   ?    ?    ? P' % badpicks[i][0]
+         elif len(badpicks[i][0]) < 4:
+             Ppattern = '%s    ?    ?    ? P' % badpicks[i][0]
+         nllocline = getPatternLine(NLLocfile, Ppattern)
+         res = nllocline.split(None)[16]
+         # get theoretical P-onset time from residuum
+         badpicks[i][1] = picks[badpicks[i][0]]['P']['mpp'] - float(res)
+
+         # get corresponding waveform stream
+         wf2pick = wf.select(station=badpicks[i][0])
+
+         # modify some picking parameters
+         pstart_old = pickparameter.getParam('pstart')
+         pstop_old = pickparameter.getParam('pstop')
+         pickwinP_old = pickparameter.getParam('pickwinP')
+         Precalcwin_old = pickparameter.getParam('Precalcwin')
+         pickparameter.setParam(pstart=badpicks[i][1] - wf2pick[0].stats.starttime \
+          - pickparameter.getParam('tlta'))
+         pickparameter.setParam(pstop=pickparameter.getParam('pstart') + \
+          (3 * pickparameter.getParam('tlta')))
+         pickparameter.setParam(pickwinP=pickparameter.getParam('pickwinP') / 2)
+         pickparameter.setParam(Precalcwin=pickparameter.getParam('Precalcwin') / 2)
+         print("iteratepicker: The following picking parameters have been modified for iterative picking:")
+         print("pstart: %fs => %fs" % (pstart_old, pickparameter.getParam('pstart')))
+         print("pstop: %fs => %fs" % (pstop_old, pickparameter.getParam('pstop')))
+         print("pickwinP: %fs => %fs" % (pickwinP_old, pickparameter.getParam('pickwinP')))
+         print("Precalcwin: %fs => %fs" % (Precalcwin_old, pickparameter.getParam('Precalcwin')))
+
+         # repick station
+         newpicks = autopickstation(wf2pick, pickparameter)
+         
+         # replace old dictionary with new one
+         picks[badpicks[i][0]] = newpicks
+ 
+    # reset temporary change of picking parameters
+    print("iteratepicker: Resetting picking parameters ...")
+    pickparameter.setParam(pstart=pstart_old)
+    pickparameter.setParam(pstop=pstop_old)
+    pickparameter.setParam(pickwinP=pickwinP_old)
+    pickparameter.setParam(Precalcwin=Precalcwin_old)
+
+    return picks
+
+
+
+
+
+
+

pylot/core/read/inputs.py

@@ -145,7 +145,7 @@ class AutoPickParameter(object):
     def setParam(self, **kwargs):
         for param, value in kwargs.items():
             self.__setitem__(param, value)
-        print(self)
+        #print(self)
 
     @staticmethod
     def _printParameterError(errmsg):