#!/usr/bin/python # -*- coding: utf-8 -*- from __future__ import print_function import argparse import glob import string import numpy as np from pylot.core.analysis.magnitude import M0Mw from pylot.core.io.data import Data from pylot.core.io.inputs import AutoPickParameter from pylot.core.loc.nll import * from pylot.core.pick.autopick import autopickevent, iteratepicker from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.version import get_git_version as _getVersionString __version__ = _getVersionString() def autoPyLoT(inputfile): """ Determine phase onsets automatically utilizing the automatic picking algorithms by Kueperkoch et al. 2010/2012. :param inputfile: path to the input file containing all parameter information for automatic picking (for formatting details, see. `~pylot.core.io.inputs.AutoPickParameter` :type inputfile: str :return: .. rubric:: Example """ splash = '''************************************\n *********autoPyLoT starting*********\n The Python picking and Location Tool\n Version {version} 2015\n \n Authors:\n S. Wehling-Benatelli (Ruhr-Universität Bochum)\n L. Küperkoch (BESTEC GmbH, Landau i. d. Pfalz)\n K. Olbert (Christian-Albrechts Universität zu Kiel)\n ***********************************'''.format(version=_getVersionString()) print(splash) # reading parameter file parameter = AutoPickParameter(inputfile) data = Data() # getting information on data structure if parameter.hasParam('datastructure'): datastructure = DATASTRUCTURE[parameter.getParam('datastructure')]() dsfields = {'root': parameter.getParam('rootpath'), 'dpath': parameter.getParam('datapath'), 'dbase': parameter.getParam('database')} exf = ['root', 'dpath', 'dbase'] if parameter.hasParam('eventID'): dsfields['eventID'] = parameter.getParam('eventID') exf.append('eventID') datastructure.modifyFields(**dsfields) datastructure.setExpandFields(exf) # check if default location routine NLLoc is available if parameter.hasParam('nllocbin'): locflag = 1 # get NLLoc-root path nllocroot = parameter.getParam('nllocroot') # get path to NLLoc executable nllocbin = parameter.getParam('nllocbin') nlloccall = '%s/NLLoc' % nllocbin # get name of phase file phasef = parameter.getParam('phasefile') phasefile = '%s/obs/%s' % (nllocroot, phasef) # get name of NLLoc-control file ctrf = parameter.getParam('ctrfile') ctrfile = '%s/run/%s' % (nllocroot, ctrf) # pattern of NLLoc ttimes from location grid ttpat = parameter.getParam('ttpatter') # pattern of NLLoc-output file nllocoutpatter = parameter.getParam('outpatter') maxnumit = 3 # maximum number of iterations for re-picking else: locflag = 0 print(" !!! ") print("!!No location routine available, autoPyLoT is running in non-location mode!!") print("!!No source parameter estimation possible!!") print(" !!! ") # multiple event processing # read each event in database datapath = datastructure.expandDataPath() if not parameter.hasParam('eventID'): for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]: 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! picks = autopickevent(wfdat, parameter) finalpicks = picks ########################################################## # locating if locflag == 1: # write phases to NLLoc-phase file picksExport(picks, 'NLLoc', phasefile) # For locating the event the NLLoc-control file has to be modified! evID = event[string.rfind(event, "/") + 1: len(events) - 1] nllocout = '%s_%s' % (evID, nllocoutpatter) # create comment line for NLLoc-control file modifyInputFile(ctrf, nllocroot, nllocout, phasef, ttpat) # 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 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!") # get NLLoc-location file locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) if len(glob.glob(locsearch)) > 0: # get latest NLLoc-location file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # calculating seismic moment Mo and moment magnitude Mw finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ nllocfile, picks, parameter.getParam('rho'), \ parameter.getParam('vp'), parameter.getParam('Qp'), \ parameter.getParam('invdir')) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") else: # get theoretical P-onset times from NLLoc-location file locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) if len(glob.glob(locsearch)) > 0: # get latest file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) nlloccounter = 0 while len(badpicks) > 0 and nlloccounter <= maxnumit: nlloccounter += 1 if nlloccounter > maxnumit: print("autoPyLoT: Number of maximum iterations reached, stop iterative picking!") break print("autoPyLoT: Starting with iteration No. %d ..." % nlloccounter) picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter) # write phases to NLLoc-phase file picksExport(picks, 'NLLoc', phasefile) # remove actual NLLoc-location file to keep only the last os.remove(nllocfile) # locate the event locate(nlloccall, ctrfile) print("autoPyLoT: Iteration No. %d finished." % nlloccounter) # get updated NLLoc-location file nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # check for bad picks badpicks = [] for key in picks: if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: badpicks.append([key, picks[key]['P']['mpp']]) print("autoPyLoT: After iteration No. %d: %d bad onsets found ..." % (nlloccounter, \ len(badpicks))) if len(badpicks) == 0: print("autoPyLoT: No more bad onsets found, stop iterative picking!") nlloccounter = maxnumit # calculating seismic moment Mo and moment magnitude Mw finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ nllocfile, picks, parameter.getParam('rho'), \ parameter.getParam('vp'), parameter.getParam('Qp'), \ parameter.getParam('invdir')) # get network moment magntiude netMw = [] for key in finalpicks.getpicdic(): if finalpicks.getpicdic()[key]['P']['Mw'] is not None: netMw.append(finalpicks.getpicdic()[key]['P']['Mw']) netMw = np.median(netMw) print("Network moment magnitude: %4.1f" % netMw) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") ########################################################## # write phase files for various location routines # HYPO71 hypo71file = '%s/autoPyLoT_HYPO71.pha' % event if hasattr(finalpicks, 'getpicdic'): if finalpicks.getpicdic() is not None: writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) data.applyEVTData(finalpicks.getpicdic()) else: writephases(picks, 'HYPO71', hypo71file) data.applyEVTData(picks) else: writephases(picks, 'HYPO71', hypo71file) data.applyEVTData(picks) fnqml = '%s/autoPyLoT' % event data.exportEvent(fnqml) 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 {0}".format(parameter.getParam('eventID'))) print(data) wfdat = data.getWFData() # all available streams ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) finalpicks = picks ########################################################## # locating if locflag == 1: # write phases to NLLoc-phase file picksExport(picks, 'NLLoc', phasefile) # For locating the event the NLLoc-control file has to be modified! nllocout = '%s_%s' % (parameter.getParam('eventID'), nllocoutpatter) # create comment line for NLLoc-control file modifyInputFile(ctrf, nllocroot, nllocout, phasef, ttpat) # 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 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!") # get NLLoc-location file locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) if len(glob.glob(locsearch)) > 0: # get latest NLLOc-location file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # calculating seismic moment Mo and moment magnitude Mw finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ nllocfile, picks, parameter.getParam('rho'), \ parameter.getParam('vp'), parameter.getParam('Qp'), \ parameter.getParam('invdir')) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") else: # get theoretical P-onset times from NLLoc-location file locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) if len(glob.glob(locsearch)) > 0: # get latest file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) nlloccounter = 0 while len(badpicks) > 0 and nlloccounter <= maxnumit: nlloccounter += 1 if nlloccounter > maxnumit: print("autoPyLoT: Number of maximum iterations reached, stop iterative picking!") break print("autoPyLoT: Starting with iteration No. %d ..." % nlloccounter) picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter) # write phases to NLLoc-phase file picksExport(picks, 'NLLoc', phasefile) # remove actual NLLoc-location file to keep only the last os.remove(nllocfile) # locate the event locate(nlloccall, ctrfile) print("autoPyLoT: Iteration No. %d finished." % nlloccounter) # get updated NLLoc-location file nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # check for bad picks badpicks = [] for key in picks: if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: badpicks.append([key, picks[key]['P']['mpp']]) print("autoPyLoT: After iteration No. %d: %d bad onsets found ..." % (nlloccounter, \ len(badpicks))) if len(badpicks) == 0: print("autoPyLoT: No more bad onsets found, stop iterative picking!") nlloccounter = maxnumit # calculating seismic moment Mo and moment magnitude Mw finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ nllocfile, picks, parameter.getParam('rho'), \ parameter.getParam('vp'), parameter.getParam('Qp'), \ parameter.getParam('invdir')) # get network moment magntiude netMw = [] for key in finalpicks.getpicdic(): if finalpicks.getpicdic()[key]['P']['Mw'] is not None: netMw.append(finalpicks.getpicdic()[key]['P']['Mw']) netMw = np.median(netMw) print("Network moment magnitude: %4.1f" % netMw) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") ########################################################## # write phase files for various location routines # HYPO71 hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID')) if hasattr(finalpicks, 'getpicdic'): if finalpicks.getpicdic() is not None: writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) data.applyEVTData(finalpicks.getpicdic()) else: writephases(picks, 'HYPO71', hypo71file) data.applyEVTData(picks) else: writephases(picks, 'HYPO71', hypo71file) data.applyEVTData(picks) fnqml = '%s/%s/autoPyLoT' % (datapath, parameter.getParam('eventID')) data.exportEvent(fnqml) 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!") endsp = '''####################################\n ************************************\n *********autoPyLoT terminates*******\n The Python picking and Location Tool\n ************************************'''.format(version=_getVersionString()) print(endsp) if __name__ == "__main__": from pylot.core.util.defaults import AUTOMATIC_DEFAULTS # parse arguments parser = argparse.ArgumentParser( description='''autoPyLoT automatically picks phase onset times using higher order statistics, autoregressive prediction and AIC''') parser.add_argument('-i', '-I', '--inputfile', type=str, action='store', help='''full path to the file containing the input parameters for autoPyLoT''', default=AUTOMATIC_DEFAULTS ) parser.add_argument('-v', '-V', '--version', action='version', version='autoPyLoT ' + __version__, help='show version information and exit') cla = parser.parse_args() autoPyLoT(str(cla.inputfile))