pylot/autoPyLoT.py

367 lines
19 KiB
Python
Raw Normal View History

#!/usr/bin/python
# -*- coding: utf-8 -*-
2016-03-30 08:14:58 +02:00
from __future__ import print_function
import argparse
import glob
import string
import os
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
import pylot.core.loc.nll as nll
import pylot.core.loc.hsat as hsat
from pylot.core.pick.autopick import autopickevent, iteratepicker
from pylot.core.util.structure import DATASTRUCTURE
2015-10-19 05:32:10 +02:00
from pylot.core.util.version import get_git_version as _getVersionString
2015-07-02 09:27:11 +02:00
__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.get('datastructure')]()
dsfields = {'root': parameter.get('rootpath'),
'dpath': parameter.get('datapath'),
'dbase': parameter.get('database')}
2015-05-20 09:59:06 +02:00
exf = ['root', 'dpath', 'dbase']
if parameter.hasParam('eventID'):
dsfields['eventID'] = parameter.get('eventID')
2015-05-20 09:59:06 +02:00
exf.append('eventID')
datastructure.modifyFields(**dsfields)
2015-05-20 09:59:06 +02:00
datastructure.setExpandFields(exf)
# check if default location routine NLLoc is available
if parameter.hasParam('nllocbin'):
locflag = 1
# get NLLoc-root path
nllocroot = parameter.get('nllocroot')
# get path to NLLoc executable
nllocbin = parameter.get('nllocbin')
nlloccall = '%s/NLLoc' % nllocbin
# get name of phase file
phasef = parameter.get('phasefile')
phasefile = '%s/obs/%s' % (nllocroot, phasef)
# get name of NLLoc-control file
ctrf = parameter.get('ctrfile')
ctrfile = '%s/run/%s' % (nllocroot, ctrf)
# pattern of NLLoc ttimes from location grid
ttpat = parameter.get('ttpatter')
# pattern of NLLoc-output file
nllocoutpatter = parameter.get('outpatter')
2016-03-30 08:14:58 +02:00
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(" !!! ")
2015-10-27 10:00:16 +01:00
# multiple event processing
# read each event in database
datapath = datastructure.expandDataPath()
2015-05-20 09:59:06 +02:00
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
2015-10-27 15:26:25 +01:00
##########################################################
# locating
if locflag == 1:
# write phases to NLLoc-phase file
nll.export(picks, phasefile)
# For locating the event the NLLoc-control file has to be modified!
2016-03-30 08:14:58 +02:00
evID = event[string.rfind(event, "/") + 1: len(events) - 1]
nllocout = '%s_%s' % (evID, nllocoutpatter)
# create comment line for NLLoc-control file
nll.modify_inputs(ctrf, nllocroot, nllocout, phasef,
ttpat)
# locate the event
nll.locate(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:
2016-03-30 08:14:58 +02:00
if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
badpicks.append([key, picks[key]['P']['mpp']])
if len(badpicks) == 0:
2016-03-30 08:14:58 +02:00
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
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \
parameter.get('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
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
nlloccounter = 0
2016-03-30 08:14:58 +02:00
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
nll.export(picks, phasefile)
# remove actual NLLoc-location file to keep only the last
os.remove(nllocfile)
# locate the event
nll.locate(ctrfile)
print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
# get updated NLLoc-location file
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# check for bad picks
badpicks = []
for key in picks:
2016-03-30 08:14:58 +02:00
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, \
2016-03-30 08:14:58 +02:00
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.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir'))
# get network moment magntiude
netMw = []
2016-03-30 08:14:58 +02:00
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!")
2015-10-27 15:26:25 +01:00
##########################################################
# write phase files for various location routines
# HYPO71
hypo71file = '%s/autoPyLoT_HYPO71.pha' % event
if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None:
hsat.export(finalpicks.getpicdic(), hypo71file)
data.applyEVTData(finalpicks.getpicdic())
else:
hsat.export(picks, hypo71file)
data.applyEVTData(picks)
fnqml = '%s/autoPyLoT' % event
data.exportEvent(fnqml)
2015-10-27 15:26:25 +01:00
endsplash = '''------------------------------------------\n'
2016-03-30 08:14:58 +02:00
-----Finished event %s!-----\n'
------------------------------------------'''.format \
(version=_getVersionString()) % evID
print(endsplash)
if locflag == 0:
print("autoPyLoT was running in non-location mode!")
# single event processing
2015-05-20 09:59:06 +02:00
else:
data.setWFData(glob.glob(os.path.join(datapath, parameter.get('eventID'), '*')))
print("Working on event {0}".format(parameter.get('eventID')))
2016-03-30 08:14:58 +02:00
print(data)
wfdat = data.getWFData() # all available streams
##########################################################
# !automated picking starts here!
picks = autopickevent(wfdat, parameter)
finalpicks = picks
2015-10-27 15:26:25 +01:00
##########################################################
# locating
if locflag == 1:
# write phases to NLLoc-phase file
nll.export(picks, phasefile)
# For locating the event the NLLoc-control file has to be modified!
nllocout = '%s_%s' % (parameter.get('eventID'), nllocoutpatter)
# create comment line for NLLoc-control file
nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, ttpat)
# locate the event
nll.locate(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:
2016-03-30 08:14:58 +02:00
if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
badpicks.append([key, picks[key]['P']['mpp']])
if len(badpicks) == 0:
2016-03-30 08:14:58 +02:00
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
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \
parameter.get('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
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
nlloccounter = 0
2016-03-30 08:14:58 +02:00
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
nll.export(picks, phasefile)
# remove actual NLLoc-location file to keep only the last
os.remove(nllocfile)
# locate the event
nll.locate(ctrfile)
print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
# get updated NLLoc-location file
2016-03-30 08:14:58 +02:00
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# check for bad picks
badpicks = []
for key in picks:
2016-03-30 08:14:58 +02:00
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, \
2016-03-30 08:14:58 +02:00
len(badpicks)))
if len(badpicks) == 0:
print("autoPyLoT: No more bad onsets found, stop iterative picking!")
nlloccounter = maxnumit
2016-03-30 08:14:58 +02:00
# calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir'))
# get network moment magntiude
netMw = []
2016-03-30 08:14:58 +02:00
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!")
2015-10-27 15:26:25 +01:00
##########################################################
# write phase files for various location routines
# HYPO71
hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.get('eventID'))
if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None:
hsat.export(finalpicks.getpicdic(), hypo71file)
data.applyEVTData(finalpicks.getpicdic())
else:
hsat.export(picks, hypo71file)
data.applyEVTData(picks)
fnqml = '%s/%s/autoPyLoT' % (datapath, parameter.get('eventID'))
data.exportEvent(fnqml)
2016-03-30 08:14:58 +02:00
endsplash = '''------------------------------------------\n'
2016-03-30 08:14:58 +02:00
-----Finished event %s!-----\n'
------------------------------------------'''.format \
(version=_getVersionString()) % parameter.get('eventID')
print(endsplash)
if locflag == 0:
print("autoPyLoT was running in non-location mode!")
2016-03-30 08:14:58 +02:00
endsp = '''####################################\n
************************************\n
*********autoPyLoT terminates*******\n
The Python picking and Location Tool\n
************************************'''.format(version=_getVersionString())
print(endsp)
2016-03-30 08:14:58 +02:00
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))