marcel/pylot
3
0
Fork 0
Code Issues 7 Pull Requests 1 Releases Wiki Activity

[refactor] major refactoring of autoPyLoT and making use of the newly introduced Magnitude objects

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2016-09-28 11:03:53 +02:00
parent 231e7dafa9
commit 5f92d1f0db
1 changed files with 56 additions and 157 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
autoPyLoT.py
View File

@ -5,17 +5,18 @@ from __future__ import print_function
import argparse import argparse
import glob import glob
import string
import os import os
import numpy as np from obspy import read_events
from pylot.core.analysis.magnitude import MomentMagnitude import pylot.core.loc.hsat as hsat
import pylot.core.loc.nll as nll
from pylot.core.analysis.magnitude import MomentMagnitude, RichterMagnitude
from pylot.core.io.data import Data from pylot.core.io.data import Data
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import AutoPickParameter
import pylot.core.loc.nll as nll from pylot.core.io.phases import add_amplitudes
import pylot.core.loc.hsat as hsat
from pylot.core.pick.autopick import autopickevent, iteratepicker from pylot.core.pick.autopick import autopickevent, iteratepicker
from pylot.core.util.dataprocessing import restitute_data, read_metadata
from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.structure import DATASTRUCTURE
from pylot.core.util.version import get_git_version as _getVersionString from pylot.core.util.version import get_git_version as _getVersionString
@ -54,6 +55,7 @@ def autoPyLoT(inputfile):
data = Data() data = Data()
evt = None
# getting information on data structure # getting information on data structure
if parameter.hasParam('datastructure'): if parameter.hasParam('datastructure'):
@ -97,19 +99,25 @@ def autoPyLoT(inputfile):
print("!!No source parameter estimation possible!!") print("!!No source parameter estimation possible!!")
print(" !!! ") print(" !!! ")
# multiple event processing
# read each event in database
datapath = datastructure.expandDataPath() datapath = datastructure.expandDataPath()
if not parameter.hasParam('eventID'): if not parameter.hasParam('eventID'):
for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]: # multiple event processing
# read each event in database
events = [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]
else:
# single event processing
events = glob.glob(os.path.join(datapath, parameter.get('eventID')))
for event in events:
data.setWFData(glob.glob(os.path.join(datapath, event, '*'))) data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
evID = os.path.split(event)[-1]
print('Working on event %s' % event) print('Working on event %s' % event)
print(data) print(data)
wfdat = data.getWFData() # all available streams wfdat = data.getWFData() # all available streams
metadata = read_metadata(parameter.get('invdir'))
corr_dat, rest_flag = restitute_data(wfdat.copy(), *metadata)
########################################################## ##########################################################
# !automated picking starts here! # !automated picking starts here!
picks = autopickevent(wfdat, parameter) picks = autopickevent(wfdat, parameter)
finalpicks = picks
########################################################## ##########################################################
# locating # locating
if locflag == 1: if locflag == 1:
@ -117,7 +125,6 @@ def autoPyLoT(inputfile):
nll.export(picks, phasefile) nll.export(picks, phasefile)
# For locating the event the NLLoc-control file has to be modified! # 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) nllocout = '%s_%s' % (evID, nllocoutpatter)
# create comment line for NLLoc-control file # create comment line for NLLoc-control file
nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, nll.modify_inputs(ctrf, nllocroot, nllocout, phasef,
@ -135,6 +142,9 @@ def autoPyLoT(inputfile):
if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4:
badpicks.append([key, picks[key]['P']['mpp']]) badpicks.append([key, picks[key]['P']['mpp']])
# TODO keep code DRY (Don't Repeat Yourself) the following part is written twice
# suggestion: delete block and modify the later similar block to work properly
if len(badpicks) == 0: if len(badpicks) == 0:
print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!")
# get NLLoc-location file # get NLLoc-location file
@ -142,11 +152,18 @@ def autoPyLoT(inputfile):
if len(glob.glob(locsearch)) > 0: if len(glob.glob(locsearch)) > 0:
# get latest NLLoc-location file if several are available # get latest NLLoc-location file if several are available
nllocfile = max(glob.glob(locsearch), key=os.path.getctime) nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
evt = read_events(nllocfile)
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
nllocfile, picks, parameter.get('rho'), \ parameter.get('Qp'), parameter.get('rho'), True, 2)
parameter.get('vp'), parameter.get('Qp'), \ local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2)
parameter.get('invdir')) for station, amplitude in local_mag.amplitudes:
picks[station]['S']['Ao'] = amplitude
evt = add_amplitudes(evt, local_mag.amplitudes)
netML = local_mag.net_magnitude()
netMw = moment_mag.net_magnitude()
evt.origins[0].magnitudes.append(netMw)
evt.origins[0].magnitudes.append(netML)
else: else:
print("autoPyLoT: No NLLoc-location file available!") print("autoPyLoT: No NLLoc-location file available!")
print("No source parameter estimation possible!") print("No source parameter estimation possible!")
@ -183,31 +200,30 @@ def autoPyLoT(inputfile):
if len(badpicks) == 0: if len(badpicks) == 0:
print("autoPyLoT: No more bad onsets found, stop iterative picking!") print("autoPyLoT: No more bad onsets found, stop iterative picking!")
nlloccounter = maxnumit nlloccounter = maxnumit
evt = read_events(nllocfile)
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'),
nllocfile, picks, parameter.get('rho'), \ parameter.get('Qp'), parameter.get('rho'), True, 2)
parameter.get('vp'), parameter.get('Qp'), \ local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2)
parameter.get('invdir')) for station, amplitude in local_mag.amplitudes:
picks[station]['S']['Ao'] = amplitude
evt = add_amplitudes(evt, local_mag.amplitudes)
netML = local_mag.net_magnitude()
# get network moment magntiude # get network moment magntiude
netMw = [] netMw = moment_mag.net_magnitude()
for key in finalpicks.getpicdic(): evt.origins[0].magnitudes.append(netMw)
if finalpicks.getpicdic()[key]['P']['Mw'] is not None: evt.origins[0].magnitudes.append(netML)
netMw.append(finalpicks.getpicdic()[key]['P']['Mw']) print("Network moment magnitude: %4.1f" % netMw.mag)
netMw = np.median(netMw)
print("Network moment magnitude: %4.1f" % netMw)
else: else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!") print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
########################################################## ##########################################################
# write phase files for various location routines # write phase files for various location routines
# HYPO71 # HYPO71
hypo71file = '%s/autoPyLoT_HYPO71.pha' % event 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) hsat.export(picks, hypo71file)
data.applyEVTData(picks) data.applyEVTData(picks)
if evt is not None:
data.applyEVTData(evt, 'event')
fnqml = '%s/autoPyLoT' % event fnqml = '%s/autoPyLoT' % event
data.exportEvent(fnqml) data.exportEvent(fnqml)
@ -219,123 +235,6 @@ def autoPyLoT(inputfile):
if locflag == 0: if locflag == 0:
print("autoPyLoT was running in non-location mode!") print("autoPyLoT was running in non-location mode!")
# single event processing
else:
data.setWFData(glob.glob(os.path.join(datapath, parameter.get('eventID'), '*')))
print("Working on event {0}".format(parameter.get('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
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:
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 = MomentMagnitude(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
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
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
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 = MomentMagnitude(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 = []
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.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)
endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n'
------------------------------------------'''.format \
(version=_getVersionString()) % parameter.get('eventID')
print(endsplash)
if locflag == 0:
print("autoPyLoT was running in non-location mode!")
endsp = '''####################################\n endsp = '''####################################\n
************************************\n ************************************\n
*********autoPyLoT terminates*******\n *********autoPyLoT terminates*******\n