2015-06-19 09:08:28 +02:00
|
|
|
#!/usr/bin/python
|
2015-04-22 12:46:24 +02:00
|
|
|
# -*- coding: utf-8 -*-
|
|
|
|
|
|
|
|
|
|
|
|
import os
|
|
|
|
import argparse
|
2015-05-20 09:38:25 +02:00
|
|
|
import glob
|
2015-04-22 12:46:24 +02:00
|
|
|
|
2015-05-29 16:26:31 +02:00
|
|
|
import matplotlib.pyplot as plt
|
|
|
|
from obspy.core import read
|
2015-04-22 12:46:24 +02:00
|
|
|
from pylot.core.util import _getVersionString
|
|
|
|
from pylot.core.read import Data, AutoPickParameter
|
2015-05-29 16:26:31 +02:00
|
|
|
from pylot.core.pick.run_autopicking import run_autopicking
|
2015-04-29 06:29:08 +02:00
|
|
|
from pylot.core.util.structure import DATASTRUCTURE
|
2015-06-19 15:27:24 +02:00
|
|
|
from pylot.core.pick.utils import wadaticheck
|
2015-04-22 12:46:24 +02:00
|
|
|
__version__ = _getVersionString()
|
|
|
|
|
|
|
|
|
2015-05-20 09:43:32 +02:00
|
|
|
def autoPyLoT(inputfile):
|
2015-04-29 06:29:08 +02:00
|
|
|
'''
|
|
|
|
Determine phase onsets automatically utilizing the automatic picking
|
2015-06-19 09:08:28 +02:00
|
|
|
algorithms by Kueperkoch et al. 2010/2012.
|
2015-04-29 06:29:08 +02:00
|
|
|
|
|
|
|
:param inputfile: path to the input file containing all parameter
|
|
|
|
information for automatic picking (for formatting details, see.
|
|
|
|
`~pylot.core.read.input.AutoPickParameter`
|
|
|
|
:type inputfile: str
|
|
|
|
:return:
|
|
|
|
|
|
|
|
.. rubric:: Example
|
|
|
|
|
|
|
|
'''
|
2015-05-20 09:38:25 +02:00
|
|
|
|
|
|
|
# reading parameter file
|
|
|
|
|
2015-04-22 12:46:24 +02:00
|
|
|
parameter = AutoPickParameter(inputfile)
|
2015-04-29 06:29:08 +02:00
|
|
|
|
2015-06-19 15:27:24 +02:00
|
|
|
# get some parameters for quality control from
|
|
|
|
# parameter input file (usually autoPyLoT.in).
|
|
|
|
wdttolerance = parameter.getParam('wdttolerance')
|
|
|
|
iplot = parameter.getParam('iplot')
|
|
|
|
|
2015-04-22 12:46:24 +02:00
|
|
|
data = Data()
|
|
|
|
|
2015-05-20 09:38:25 +02:00
|
|
|
# getting information on data structure
|
|
|
|
|
2015-04-29 06:29:08 +02:00
|
|
|
if parameter.hasParam('datastructure'):
|
|
|
|
datastructure = DATASTRUCTURE[parameter.getParam('datastructure')]()
|
2015-05-20 09:38:25 +02:00
|
|
|
dsfields = {'root':parameter.getParam('rootpath'),
|
|
|
|
'dpath':parameter.getParam('datapath'),
|
|
|
|
'dbase':parameter.getParam('database')}
|
2015-04-29 07:57:52 +02:00
|
|
|
|
2015-05-20 09:59:06 +02:00
|
|
|
exf = ['root', 'dpath', 'dbase']
|
|
|
|
|
|
|
|
if parameter.hasParam('eventID'):
|
|
|
|
dsfields['eventID'] = parameter.getParam('eventID')
|
|
|
|
exf.append('eventID')
|
2015-04-22 12:46:24 +02:00
|
|
|
|
2015-06-01 16:28:27 +02:00
|
|
|
datastructure.modifyFields(**dsfields)
|
2015-05-20 09:59:06 +02:00
|
|
|
datastructure.setExpandFields(exf)
|
2015-04-22 12:46:24 +02:00
|
|
|
|
2015-06-19 09:08:28 +02:00
|
|
|
# multiple event processing
|
2015-05-29 16:26:31 +02:00
|
|
|
# read each event in database
|
2015-05-20 09:38:25 +02:00
|
|
|
datapath = datastructure.expandDataPath()
|
2015-05-20 09:59:06 +02:00
|
|
|
if not parameter.hasParam('eventID'):
|
2015-05-29 16:26:31 +02:00
|
|
|
for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]:
|
2015-05-20 11:29:55 +02:00
|
|
|
data.setWFData(glob.glob(os.path.join(datapath, event, '*')))
|
2015-06-01 16:28:27 +02:00
|
|
|
print 'Working on event %s' %event
|
|
|
|
print data
|
|
|
|
|
|
|
|
wfdat = data.getWFData() # all available streams
|
|
|
|
##########################################################
|
|
|
|
# !automated picking starts here!
|
|
|
|
procstats = []
|
2015-06-19 15:27:24 +02:00
|
|
|
# initialize dictionary for onsets
|
|
|
|
picks = None
|
|
|
|
station = wfdat[0].stats.station
|
|
|
|
allonsets = {station: picks}
|
2015-06-01 16:28:27 +02:00
|
|
|
for i in range(len(wfdat)):
|
|
|
|
stationID = wfdat[i].stats.station
|
2015-06-19 09:08:28 +02:00
|
|
|
# check if station has already been processed
|
2015-06-01 16:28:27 +02:00
|
|
|
if stationID not in procstats:
|
|
|
|
procstats.append(stationID)
|
2015-06-19 09:08:28 +02:00
|
|
|
# find corresponding streams
|
2015-06-01 16:28:27 +02:00
|
|
|
statdat = wfdat.select(station=stationID)
|
2015-06-19 15:27:24 +02:00
|
|
|
######################################################
|
2015-06-19 09:08:28 +02:00
|
|
|
# get onset times and corresponding picking errors
|
|
|
|
picks = run_autopicking(statdat, parameter)
|
2015-06-19 15:27:24 +02:00
|
|
|
######################################################
|
|
|
|
# add station and corresponding onsets to dictionary
|
|
|
|
station = stationID
|
|
|
|
allonsets[station] = picks
|
2015-06-19 09:08:28 +02:00
|
|
|
|
|
|
|
# quality control
|
|
|
|
# jackknife on P onset times
|
2015-06-19 15:27:24 +02:00
|
|
|
# check S-P times (Wadati)
|
|
|
|
checkedonsets = wadaticheck(allonsets, wdttolerance, iplot)
|
2015-06-19 09:08:28 +02:00
|
|
|
# jackknife on S onset times
|
2015-06-02 13:42:44 +02:00
|
|
|
print '------------------------------------------'
|
|
|
|
print '-----Finished event %s!-----' % event
|
|
|
|
print '------------------------------------------'
|
2015-06-01 16:28:27 +02:00
|
|
|
|
2015-06-19 09:08:28 +02:00
|
|
|
# single event processing
|
2015-05-20 09:59:06 +02:00
|
|
|
else:
|
2015-06-02 13:43:39 +02:00
|
|
|
data.setWFData(glob.glob(os.path.join(datapath, parameter.getParam('eventID'), '*')))
|
2015-06-01 16:28:27 +02:00
|
|
|
print 'Working on event ', parameter.getParam('eventID')
|
2015-05-20 10:56:53 +02:00
|
|
|
print data
|
2015-05-29 16:26:31 +02:00
|
|
|
|
2015-06-01 16:28:27 +02:00
|
|
|
wfdat = data.getWFData() # all available streams
|
|
|
|
##########################################################
|
|
|
|
# !automated picking starts here!
|
|
|
|
procstats = []
|
2015-06-19 09:08:28 +02:00
|
|
|
# initialize dictionary for onsets
|
|
|
|
picks = None
|
|
|
|
station = wfdat[0].stats.station
|
|
|
|
allonsets = {station: picks}
|
2015-06-01 16:28:27 +02:00
|
|
|
for i in range(len(wfdat)):
|
2015-06-19 15:27:24 +02:00
|
|
|
#for i in range(0,5):
|
2015-06-01 16:28:27 +02:00
|
|
|
stationID = wfdat[i].stats.station
|
|
|
|
#check if station has already been processed
|
|
|
|
if stationID not in procstats:
|
|
|
|
procstats.append(stationID)
|
|
|
|
#find corresponding streams
|
|
|
|
statdat = wfdat.select(station=stationID)
|
2015-06-19 15:27:24 +02:00
|
|
|
######################################################
|
|
|
|
# get onset times and corresponding picking parameters
|
2015-06-19 09:08:28 +02:00
|
|
|
picks = run_autopicking(statdat, parameter)
|
2015-06-19 15:27:24 +02:00
|
|
|
######################################################
|
2015-06-19 09:08:28 +02:00
|
|
|
# add station and corresponding onsets to dictionary
|
2015-06-19 15:27:24 +02:00
|
|
|
station = stationID
|
2015-06-19 09:08:28 +02:00
|
|
|
allonsets[station] = picks
|
2015-06-19 15:27:24 +02:00
|
|
|
|
2015-06-19 09:08:28 +02:00
|
|
|
#quality control
|
|
|
|
#jackknife on P onset times
|
2015-06-19 15:27:24 +02:00
|
|
|
#check S-P times (Wadati)
|
|
|
|
checkedonsets = wadaticheck(allonsets, wdttolerance, iplot)
|
2015-06-19 09:08:28 +02:00
|
|
|
#jackknife on S onset times
|
2015-06-02 13:42:44 +02:00
|
|
|
print '------------------------------------------'
|
|
|
|
print '-------Finished event %s!-------' % parameter.getParam('eventID')
|
|
|
|
print '------------------------------------------'
|
2015-05-29 16:26:31 +02:00
|
|
|
|
2015-04-29 07:57:52 +02:00
|
|
|
|
2015-04-22 12:46:24 +02:00
|
|
|
if __name__ == "__main__":
|
|
|
|
# parse arguments
|
|
|
|
parser = argparse.ArgumentParser(
|
2015-05-29 16:26:31 +02:00
|
|
|
description='''autoPyLoT automatically picks phase onset times using higher order statistics,
|
|
|
|
autoregressive prediction and AIC''')
|
2015-04-22 12:46:24 +02:00
|
|
|
|
|
|
|
parser.add_argument('-i', '-I', '--inputfile', type=str,
|
|
|
|
action='store',
|
|
|
|
help='''full path to the file containing the input
|
|
|
|
parameters for autoPyLoT''',
|
2015-05-04 05:31:10 +02:00
|
|
|
default=os.path.join(os.path.expanduser('~'), '.pylot',
|
|
|
|
'autoPyLoT.in')
|
2015-04-22 12:46:24 +02:00
|
|
|
)
|
|
|
|
parser.add_argument('-v', '-V', '--version', action='version',
|
|
|
|
version='autoPyLoT ' + __version__,
|
|
|
|
help='show version information and exit')
|
|
|
|
|
|
|
|
cla = parser.parse_args()
|
|
|
|
|
2015-05-20 10:01:19 +02:00
|
|
|
autoPyLoT(str(cla.inputfile))
|