diff --git a/autoPyLoT.py b/autoPyLoT.py index 71a8031f..e3f696f5 100755 --- a/autoPyLoT.py +++ b/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 '------------------------------------------' - print '-----Finished event %s!-----' % event - print '------------------------------------------' + endsplash = '''------------------------------------------\n' + -----Finished event %s!-----\n' + ------------------------------------------'''.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 '------------------------------------------' - print '-------Finished event %s!-------' % parameter.getParam('eventID') - print '------------------------------------------' - - print '####################################' - print '************************************' - print '*********autoPyLoT terminates*******' - print 'The Python picking and Location Tool' - print '************************************' + endsp = '''####################################\n + ************************************\n + *********autoPyLoT terminates*******\n + The Python picking and Location Tool\n + ************************************'''.format(version=_getVersionString()) + print(endsp) if __name__ == "__main__": # parse arguments diff --git a/pylot/core/loc/nll.py b/pylot/core/loc/nll.py index f411906c..3ecd1211 100644 --- a/pylot/core/loc/nll.py +++ b/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) - locfiles = 'LOCFILES %s NLLOC_OBS %s %s 0' % (phasefile, tttable, nllocout) + nllocout = os.path.join(root,'loc', outpath) + 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): ''' diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 5a5cc4e5..0e1c3030 100755 --- a/pylot/core/pick/autopick.py +++ b/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 - specpara = DCfc(z_copy, mpickP, calcwin, iplot) - w0 = specpara.getw0() - fc = specpara.getfc() + # calculate source spectrum and get w0 and fc + specpara = DCfc(z_copy, mpickP, calcwin, iplot) + w0 = specpara.getw0() + 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 diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 8319d90a..9887c78d 100644 --- a/pylot/core/pick/utils.py +++ b/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/4 is assumed as time difference between most likely and earliest possible pick! - EPick = Pick1 - T0 / 2 + T0 = np.mean(np.diff(zc)) * X[0].stats.delta # this is half wave length! + EPick = Pick1 - T0 # half wavelength as suggested by Diehl et al. # 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 diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index bbfe4d4d..483aabcb 100644 --- a/pylot/core/util/utils.py +++ b/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()