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

Merge branch 'develop' of ariadne.geophysik.ruhr-uni-bochum.de:/data/git/pylot into develop

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2015-11-23 11:41:14 +01:00
commit 1f1d0aa118
6 changed files with 161 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

51
autoPyLoT.py
View File

@ -6,11 +6,11 @@ import argparse
import glob
import subprocess
import string
from obspy.core import read
from obspy.core import read, UTCDateTime
from pylot.core.read.data import Data
from pylot.core.read.inputs import AutoPickParameter
from pylot.core.util.structure import DATASTRUCTURE
from pylot.core.pick.autopick import autopickevent
from pylot.core.pick.autopick import autopickevent, iteratepicker
from pylot.core.loc.nll import *
from pylot.core.util.version import get_git_version as _getVersionString
@ -121,7 +121,28 @@ def autoPyLoT(inputfile):
# !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
# 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!")
else:
# get theoretical P-onset times from NLLoc-location file
locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
# get latest file if several are available
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
if os.path.isfile(nllocfile):
picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
# write phases to NLLoc-phase file
picksExport(picks, 'NLLoc', phasefile)
# locate the event
locate(nlloccall, ctrfile)
else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
##########################################################
# write phase files for various location routines
# HYPO71
@ -160,10 +181,30 @@ def autoPyLoT(inputfile):
# 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
# 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!")
else:
# get theoretical P-onset times from NLLoc-location file
locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout)
# get latest file if several are available
nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
if os.path.isfile(nllocfile):
picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
# write phases to NLLoc-phase file
picksExport(picks, 'NLLoc', phasefile)
# locate the event
locate(nlloccall, ctrfile)
else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
##########################################################
# write phase files for various location routines
# HYPO71

18
autoPyLoT_local.in
View File

@ -17,7 +17,7 @@ PILOT #datastructure#%choose data structure
/home/ludger/NLLOC/Insheim #nllocroot# %root of NLLoc-processing directory
AUTOPHASES.obs #phasefile# %name of autoPyLoT-output phase file for NLLoc
%(in nllocroot/obs)
Insheim_min1d2015_auto.in #locfile# %name of autoPyLoT-output control file for NLLoc
Insheim_min1d2015_auto.in #ctrfile# %name of autoPyLoT-output control file for NLLoc
%(in nllocroot/run)
ttime #ttpatter# %pattern of NLLoc ttimes from grid
%(in nllocroot/times)
@ -41,7 +41,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file co
#common settings picker#
15 #pstart# %start time [s] for calculating CF for P-picking
60 #pstop# %end time [s] for calculating CF for P-picking
-1.0 #sstart# %start time [s] after or before(-) P-onset for calculating CF for S-picking
-1.0 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking
7 #sstop# %end time [s] after P-onset for calculating CF for S-picking
2 20 #bpz1# %lower/upper corner freq. of first band pass filter Z-comp. [Hz]
2 30 #bpz2# %lower/upper corner freq. of second band pass filter Z-comp. [Hz]
@ -51,7 +51,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file co
%!!Be careful when editing the following!!
#Z-component#
HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3)
7 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s]
7.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s]
4 #hosorder# %for HOS-picker, order of Higher Order Statistics
2 #Parorder# %for AR-picker, order of AR process of Z-component
1.2 #tdet1z# %for AR-picker, length of AR determination window [s] for Z-component, 1st pick
@ -60,8 +60,8 @@ HOS #algoP# %choose algorithm for P-onset
0.2 #tpred2z# %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick
0.001 #addnoise# %add noise to seismogram for stable AR prediction
3 0.1 0.5 0.5 #tsnrz# %for HOS/AR, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
3 #pickwinP# %for initial AIC pick, length of P-pick window [s]
8 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
3.0 #pickwinP# %for initial AIC pick, length of P-pick window [s]
8.0 #Precalcwin# %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)
0 #peps4aic# %for HOS/AR, artificial uplift of samples of AIC-function (P)
0.2 #aictsmooth# %for HOS/AR, take average of samples for smoothing of AIC-function [s]
0.1 #tsmoothP# %for HOS/AR, take average of samples for smoothing CF [s]
@ -74,8 +74,8 @@ ARH #algoS# %choose algorithm for S-onset
0.6 #tdet2h# %for HOS/AR, length of AR-determinaton window [s], H-components, 2nd pick
0.3 #tpred2h# %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick
4 #Sarorder# %for AR-picker, order of AR process of H-components
6 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
3 #pickwinS# %for initial AIC pick, length of S-pick window [s]
6.0 #Srecalcwin# %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)
3.0 #pickwinS# %for initial AIC pick, length of S-pick window [s]
2 0.2 1.5 0.5 #tsnrh# %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise,tsafetey,tsignal,tslope] [s]
0.05 #aictsmoothS# %for AIC-picker, take average of samples for smoothing of AIC-function [s]
0.02 #tsmoothS# %for AR-picker, take average of samples for smoothing CF [s] (S)
@ -83,7 +83,7 @@ ARH #algoS# %choose algorithm for S-onset
0.4 #ausS# %for HOS/AR, artificial uplift of samples (aus) of CF (S)
1.5 #nfacS# %for AR-picker, noise factor for noise level determination (S)
%first-motion picker%
1 #minfmweight# %minimum required p weight for first-motion determination
1 #minfmweight# %minimum required P weight for first-motion determination
2 #minFMSNR# %miniumum required SNR for first-motion determination
0.2 #fmpickwin# %pick window around P onset for calculating zero crossings
%quality assessment%
@ -92,7 +92,7 @@ ARH #algoS# %choose algorithm for S-onset
0.04 0.08 0.16 0.32 #timeerrorsS# %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S
10 #minAICPslope# %below this slope [counts/s] the initial P pick is rejected
1.2 #minAICPSNR# %below this SNR the initial P pick is rejected
6 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected
3 #minAICSslope# %below this slope [counts/s] the initial S pick is rejected
1.5 #minAICSSNR# %below this SNR the initial S pick is rejected
#check duration of signal using envelope function#
5 #minsiglength# %minimum required length of signal [s]

2
pylot/core/active/fmtomoUtils.py
View File

@ -341,7 +341,7 @@ def addCheckerboard(spacing = 10., pertubation = 0.1, inputfile = 'vgrids.in',
_update_progress(progress)
print('Added checkerboard to the grid in file %s with a spacing of %s and a pertubation of %s %%. '
'Outputfile: %s.'%(inputfile, spacing, pertubation, outputfile))
'Outputfile: %s.'%(inputfile, spacing, pertubation*100, outputfile))
outfile.close()
def addBox(x = (None, None), y = (None, None), z = (None, None),

32
pylot/core/active/seismicArrayPreparation.py
View File

@ -381,14 +381,34 @@ class SeisArray(object):
return surface
def generateFMTOMOinputFromArray(self, nRP, nThetaP, nPhiP, nRI, nThetaI, nPhiI,
Rbt, cushionfactor, interpolationMethod = 'linear',
def generateFMTOMOinputFromArray(self, nPointsPropgrid, nPointsInvgrid,
zBotTop, cushionfactor, interpolationMethod = 'linear',
customgrid = 'mygrid.in', writeVTK = True):
'''
Generate FMTOMO input files from the SeisArray dimensions.
Generates: vgrids.in, interfaces.in, propgrid.in
:param: nPointsPropgrid, number of points in each direction of the propagation grid (z, y, x)
:type: tuple
:param: nPointsInvgrid, number of points in each direction of the inversion grid (z, y, x)
:type: tuple
:param: zBotTop, (bottom, top) dimensions of the model
:type: tuple
:param: cushionfactor, adds cushioning around the model (0.1 = 10%)
:type: float
'''
nRP, nThetaP, nPhiP = nPointsPropgrid
nRI, nThetaI, nPhiI = nPointsInvgrid
print('\n------------------------------------------------------------')
print('Automatically generating input for FMTOMO from array size.')
print('Propgrid: nR = %s, nTheta = %s, nPhi = %s'%(nRP, nThetaP, nPhiP))
print('Interpolation Grid and Interfaces Grid: nR = %s, nTheta = %s, nPhi = %s'%(nRI, nThetaI, nPhiI))
print('Bottom and Top of model: (%s, %s)'%(Rbt[0], Rbt[1]))
print('Bottom and Top of model: (%s, %s)'%(zBotTop[0], zBotTop[1]))
print('Method: %s, customgrid = %s'%(interpolationMethod, customgrid))
print('------------------------------------------------------------')
@ -398,13 +418,13 @@ class SeisArray(object):
z.append(point[2])
return min(z)
self.generatePropgrid(nThetaP, nPhiP, nRP, Rbt, cushionfactor = cushionfactor,
self.generatePropgrid(nThetaP, nPhiP, nRP, zBotTop, cushionfactor = cushionfactor,
cushionpropgrid = 0.05)
surface = self.generateVgrid(nThetaI, nPhiI, nRI, Rbt, method = interpolationMethod,
surface = self.generateVgrid(nThetaI, nPhiI, nRI, zBotTop, method = interpolationMethod,
cushionfactor = cushionfactor, infilename = customgrid,
returnTopo = True)
depthmax = abs(Rbt[0] - getZmin(surface)) - 1.0 # cushioning for the bottom interface
depthmax = abs(zBotTop[0] - getZmin(surface)) - 1.0 # cushioning for the bottom interface
interf1, interf2 = self.generateInterfaces(nThetaI, nPhiI, depthmax, cushionfactor = cushionfactor,
returnInterfaces = True, method = interpolationMethod)

82
pylot/core/pick/autopick.py
View File

@ -3,7 +3,7 @@
"""
Function to run automated picking algorithms using AIC,
HOS and AR prediction. Uses object CharFuns and Picker and
HOS and AR prediction. Uses objects CharFuns and Picker and
function conglomerate utils.
:author: MAGS2 EP3 working group / Ludger Kueperkoch
@ -16,6 +16,7 @@ from pylot.core.pick.Picker import AICPicker, PragPicker
from pylot.core.pick.CharFuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker,\
getSNR, fmpicker, checkPonsets, wadaticheck, crossings_nonzero_all
from pylot.core.util.utils import getPatternLine
from pylot.core.read.data import Data
from pylot.core.analysis.magnitude import WApp, DCfc
@ -316,12 +317,11 @@ def autopickstation(wfstream, pickparam):
# from P pulse
# initialize Data object
data = Data()
[corzdat, restflag] = data.restituteWFData(invdir, zdat)
z_copy = zdat.copy()
[corzdat, restflag] = data.restituteWFData(invdir, z_copy)
if restflag == 1:
# integrate to displacement
corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta)
# class needs stream object => build it
z_copy = zdat.copy()
z_copy[0].data = corintzdat
# largest detectable period == window length
# after P pulse for calculating source spectrum
@ -799,3 +799,77 @@ def autopickstation(wfstream, pickparam):
picks[phase]['Ao'] = Ao
return picks
def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
'''
Repicking of bad onsets. Uses theoretical onset times from NLLoc-location file.
:param wf: waveform, obspy stream object
:param NLLocfile: path/name of NLLoc-location file
:param picks: dictionary of available onset times
:param badpicks: picks to be repicked
:param pickparameter: picking parameters from autoPyLoT-input file
'''
print("#######################################################")
print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \
% badpicks
newpicks = {}
for i in range(0, len(badpicks)):
if len(badpicks[i][0]) > 4:
Ppattern = '%s ? ? ? P' % badpicks[i][0]
elif len(badpicks[i][0]) == 4:
Ppattern = '%s ? ? ? P' % badpicks[i][0]
elif len(badpicks[i][0]) < 4:
Ppattern = '%s ? ? ? P' % badpicks[i][0]
nllocline = getPatternLine(NLLocfile, Ppattern)
res = nllocline.split(None)[16]
# get theoretical P-onset time from residuum
badpicks[i][1] = picks[badpicks[i][0]]['P']['mpp'] - float(res)
# get corresponding waveform stream
wf2pick = wf.select(station=badpicks[i][0])
# modify some picking parameters
pstart_old = pickparameter.getParam('pstart')
pstop_old = pickparameter.getParam('pstop')
pickwinP_old = pickparameter.getParam('pickwinP')
Precalcwin_old = pickparameter.getParam('Precalcwin')
pickparameter.setParam(pstart=badpicks[i][1] - wf2pick[0].stats.starttime \
- pickparameter.getParam('tlta'))
pickparameter.setParam(pstop=pickparameter.getParam('pstart') + \
(3 * pickparameter.getParam('tlta')))
pickparameter.setParam(pickwinP=pickparameter.getParam('pickwinP') / 2)
pickparameter.setParam(Precalcwin=pickparameter.getParam('Precalcwin') / 2)
print("iteratepicker: The following picking parameters have been modified for iterative picking:")
print("pstart: %fs => %fs" % (pstart_old, pickparameter.getParam('pstart')))
print("pstop: %fs => %fs" % (pstop_old, pickparameter.getParam('pstop')))
print("pickwinP: %fs => %fs" % (pickwinP_old, pickparameter.getParam('pickwinP')))
print("Precalcwin: %fs => %fs" % (Precalcwin_old, pickparameter.getParam('Precalcwin')))
# repick station
newpicks = autopickstation(wf2pick, pickparameter)
# replace old dictionary with new one
picks[badpicks[i][0]] = newpicks
# reset temporary change of picking parameters
print("iteratepicker: Resetting picking parameters ...")
pickparameter.setParam(pstart=pstart_old)
pickparameter.setParam(pstop=pstop_old)
pickparameter.setParam(pickwinP=pickwinP_old)
pickparameter.setParam(Precalcwin=Precalcwin_old)
return picks

2
pylot/core/read/inputs.py
View File

@ -145,7 +145,7 @@ class AutoPickParameter(object):
def setParam(self, **kwargs):
for param, value in kwargs.items():
self.__setitem__(param, value)
print(self)
#print(self)
@staticmethod
def _printParameterError(errmsg):