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

autoPyLoT.py

@@ -6,7 +6,7 @@ import os
 import argparse
 import glob
 import subprocess
-
+import string
 from obspy.core import read
 from pylot.core.read.data import Data
 from pylot.core.read.inputs import AutoPickParameter
@@ -88,7 +88,9 @@ def autoPyLoT(inputfile):
             nllocoutpatter = parameter.getParam('outpatter')
         else:
             locflag = 0
-            print ("!!No location routine available, autoPyLoT just picks the events without locating them!!")
+            print ("                 !!!              ")
+            print ("!!No location routine available, autoPyLoT is running in non-location mode!!")
+            print ("                 !!!              ")
 
 
         # multiple event processing
@@ -99,7 +101,6 @@ def autoPyLoT(inputfile):
                 data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
                 print('Working on event %s' % event)
                 print(data)
-
                 wfdat = data.getWFData()  # all available streams
                 ##########################################################
                 # !automated picking starts here!
@@ -114,7 +115,8 @@ def autoPyLoT(inputfile):
                     # For locating the event the NLLoc-control file has to be modified!
                     # create comment line for NLLoc-control file
                     # NLLoc-output file
-                    nllocout = '%s/loc/%s_%s' % (nllocroot, event, nllocoutpatter)
+                    evID = event[string.rfind(event, "/") + 1 : len(events) - 1]
+                    nllocout = '%s/loc/%s_%s' % (nllocroot, evID, nllocoutpatter)
                     locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0' % (phasefile, ttpatter, nllocout)
                     print ("Modifying  NLLoc-control file %s ..." % locfile)
                     # modification of NLLoc-control file
@@ -130,20 +132,28 @@ def autoPyLoT(inputfile):
 
                     # locate the event
                     subprocess.call([nlloccall, locfile])
+
+                    # !iterative picking if traces remained unpicked or with bad picks!
+                    # get theoretical onset times for picks with weights >= 4
+                    # in order to reprocess them using smaller time windows
                 ##########################################################
                 # write phase files for various location routines
                 # HYPO71
-                hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, eventID)
+                hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, evID)
                 writephases(picks, 'HYPO71', hypo71file)
 
-                print '------------------------------------------'
+                endsplash = '''------------------------------------------\n'
-                print '-----Finished event %s!-----' % event
+                               -----Finished event %s!-----\n' 
-                print '------------------------------------------'
+                               ------------------------------------------'''.format \
+                                (version=_getVersionString()) % evID
+                print(endsplash)
+                if locflag == 0:
+                    print("autoPyLoT was running in non-location mode!")
 
         # single event processing
         else:
             data.setWFData(glob.glob(os.path.join(datapath, parameter.getParam('eventID'), '*')))
-            print 'Working on event ', parameter.getParam('eventID')
+            print("Working on event "), parameter.getParam('eventID')
             print data
 
             wfdat = data.getWFData()  # all available streams
@@ -175,22 +185,30 @@ def autoPyLoT(inputfile):
 
                 # locate the event
                 subprocess.call([nlloccall, locfile])
+
+                # !iterative picking if traces remained unpicked or with bad picks!
+                # get theoretical onset times for picks with weights >= 4
+                # in order to reprocess them using smaller time windows
             ##########################################################
             # write phase files for various location routines
             # HYPO71
             hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID'))
             writephases(picks, 'HYPO71', hypo71file)
            
+            endsplash = '''------------------------------------------\n'
+                           -----Finished event %s!-----\n' 
+                           ------------------------------------------'''.format \
+                            (version=_getVersionString()) % parameter.getParam('eventID')
+            print(endsplash)
+            if locflag == 0:
+                print("autoPyLoT was running in non-location mode!")
                    
-            print '------------------------------------------'
+    endsp = '''####################################\n
-            print '-------Finished event %s!-------' % parameter.getParam('eventID')
+               ************************************\n
-            print '------------------------------------------'
+               *********autoPyLoT terminates*******\n
-
+               The Python picking and Location Tool\n
-    print '####################################'
+               ************************************'''.format(version=_getVersionString())
-    print '************************************'
+    print(endsp)
-    print '*********autoPyLoT terminates*******'
-    print 'The Python picking and Location Tool'
-    print '************************************'
 
 if __name__ == "__main__":
     # parse arguments

pylot/core/loc/nll.py

@@ -2,8 +2,13 @@
 # -*- coding: utf-8 -*-
 
 import subprocess
+import os
 from obspy.core.event import readEvents
 from pylot.core.pick.utils import writephases
+from pylot.core.util.utils import getPatternLine
+from pylot.core.util.version import get_git_version as _getVersionString
+
+__version__ = _getVersionString()
 
 def picksExport(picks, phasefile):
     '''
@@ -16,7 +21,7 @@ def picksExport(picks, phasefile):
     # write phases to NLLoc-phase file
     writephases(picks, 'NLLoc', phasefile)
 
-def modfiyInputFile(fn, root, outpath, phasefile, tttable):
+def modfiyInputFile(fn, root, outpath, phasefn, tttn):
     '''
 
     :param fn:
@@ -29,18 +34,21 @@ def modfiyInputFile(fn, root, outpath, phasefile, tttable):
     # For locating the event we have to modify the NLLoc-control file!
     # create comment line for NLLoc-control file NLLoc-output file
     print ("Modifying  NLLoc-control file %s ..." % fn)
-    nllocout = '%s/loc/%s_%s' % (root, outpath)
+    nllocout = os.path.join(root,'loc', outpath)
-    locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0' % (phasefile, tttable, nllocout)
+    phasefile = os.path.join(root, 'obs', phasefn)
+    tttable = os.path.join(root, 'time', tttn)
+    locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0\n' % (phasefile, tttable, nllocout)
 
     # modification of NLLoc-control file
+    curlocfiles = getPatternLine(fn, 'LOCFILES')
     nllfile = open(fn, 'r')
     filedata = nllfile.read()
     if filedata.find(locfiles) < 0:
         # replace old command
-        filedata = filedata.replace('LOCFILES', locfiles)
+        filedata = filedata.replace(curlocfiles, locfiles)
         nllfile = open(fn, 'w')
         nllfile.write(filedata)
-        nllfile.close()
+    nllfile.close()
 
 def locate(call, fnin):
     '''

pylot/core/pick/autopick.py

@@ -28,7 +28,6 @@ def autopickevent(data, param):
     wdttolerance = param.getParam('wdttolerance')
     mdttolerance = param.getParam('mdttolerance')
     iplot = param.getParam('iplot')
-
     for n in range(len(data)):
         station = data[n].stats.station
         if station not in stations:
@@ -141,6 +140,8 @@ def autopickstation(wfstream, pickparam):
 
     # split components
     zdat = wfstream.select(component="Z")
+    if len(zdat) == 0:  # check for other components
+        zdat = wfstream.select(component="3")
     edat = wfstream.select(component="E")
     if len(edat) == 0:  # check for other components
         edat = wfstream.select(component="2")
@@ -330,17 +331,17 @@ def autopickstation(wfstream, pickparam):
                     # calculate spectrum using only first cycles of
                     # waveform after P onset!
                     zc = crossings_nonzero_all(wfzc)
-                    if np.size(zc) == 0:
+                    if np.size(zc) == 0 or len(zc) <= 3:
                         print ("Something is wrong with the waveform, "
                                "no zero crossings derived!")
                         print ("Cannot calculate source spectrum!")
                     else:
                         index = min([3, len(zc) - 1])
                         calcwin = (zc[index] - zc[0]) * z_copy[0].stats.delta
-                    # calculate source spectrum and get w0 and fc
+                        # calculate source spectrum and get w0 and fc
-                    specpara = DCfc(z_copy, mpickP, calcwin, iplot)
+                        specpara = DCfc(z_copy, mpickP, calcwin, iplot)
-                    w0 = specpara.getw0()
+                        w0 = specpara.getw0()
-                    fc = specpara.getfc()
+                        fc = specpara.getfc()
 
                 print ("autopickstation: P-weight: %d, SNR: %f, SNR[dB]: %f, "
                        "Polarity: %s" % (Pweight, SNRP, SNRPdB, FM))
@@ -759,12 +760,12 @@ def autopickstation(wfstream, pickparam):
             plt.close()
     ##########################################################################
     # calculate "real" onset times
-    if mpickP is not None:
+    if mpickP is not None and epickP is not None and mpickP is not None:
         lpickP = zdat[0].stats.starttime + lpickP
         epickP = zdat[0].stats.starttime + epickP
         mpickP = zdat[0].stats.starttime + mpickP
 
-    if mpickS is not None:
+    if mpickS is not None and epickS is not None and mpickS is not None:
         lpickS = edat[0].stats.starttime + lpickS
         epickS = edat[0].stats.starttime + epickS
         mpickS = edat[0].stats.starttime + mpickS

pylot/core/pick/utils.py

@@ -8,7 +8,7 @@
 
    :author: Ludger Kueperkoch / MAGS2 EP3 working group
 """
-import pdb
+
 import numpy as np
 import matplotlib.pyplot as plt
 from obspy.core import Stream, UTCDateTime
@@ -91,9 +91,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealthMode = False):
         # determine all zero crossings in signal window (demeaned)
         zc = crossings_nonzero_all(x[pis] - x[pis].mean())
         # calculate mean half period T0 of signal as the average of the
-        T0 = np.mean(np.diff(zc)) * X[0].stats.delta  # this is half wave length
+        T0 = np.mean(np.diff(zc)) * X[0].stats.delta  # this is half wave length!
-        # T0/4 is assumed as time difference between most likely and earliest possible pick!
+        EPick = Pick1 - T0 # half wavelength as suggested by Diehl et al.
-        EPick = Pick1 - T0 / 2
 
 
     # get symmetric pick error as mean from earliest and latest possible pick
@@ -867,6 +866,8 @@ def checkZ4S(X, pick, zfac, checkwin, iplot):
 
     # split components
     zdat = X.select(component="Z")
+    if len(zdat) == 0:  # check for other components
+        zdat = X.select(component="3")
     edat = X.select(component="E")
     if len(edat) == 0:  # check for other components
         edat = X.select(component="2")
@@ -962,6 +963,8 @@ def writephases(arrivals, fformat, filename):
         for key in arrivals:
             if arrivals[key]['P']['weight'] < 4:
                 fm = arrivals[key]['P']['fm']
+                if fm == None:
+                   fm = '?'
                 onset = arrivals[key]['P']['mpp']
                 year = onset.year
                 month = onset.month

pylot/core/util/utils.py

@@ -1,7 +1,5 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-#
-# -*- coding: utf-8 -*-
 
 import os
 import subprocess
@@ -69,6 +67,30 @@ def getHash(time):
 def getOwner(fn):
     return pwd.getpwuid(os.stat(fn).st_uid).pw_name
 
+def getPatternLine(fn, pattern):
+    """
+    Takes a file name and a pattern string to search for in the file and
+    returns the first line which contains the pattern string otherwise None.
+
+    :param fn: file name
+    :type fn: str
+    :param pattern: pattern string to search for
+    :type pattern: str
+    :return: the complete line containing pattern or None
+
+    >>> getPatternLine('utils.py', 'python')
+    '#!/usr/bin/env python\\n'
+    >>> print(getPatternLine('version.py', 'palindrome'))
+    None
+    """
+    fobj = open(fn, 'r')
+    for line in fobj.readlines():
+        if pattern in line:
+            fobj.close()
+            return line
+
+    return None
+
 
 def prepTimeAxis(stime, trace):
     nsamp = trace.stats.npts
@@ -342,3 +364,7 @@ def createAmplitude(pickID, amp, unit, category, cinfo):
     amplitude.type = ope.AmplitudeCategory(category)
     amplitude.pick_id = pickID
     return amplitude
+
+if __name__ == "__main__":
+    import doctest
+    doctest.testmod()