From 7f85ed99c0d9685b3e087f4ab808524e4756aa5f Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Thu, 19 Nov 2015 11:38:54 +0100 Subject: [PATCH 1/8] bugfix textoutput --- pylot/core/active/fmtomoUtils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/active/fmtomoUtils.py b/pylot/core/active/fmtomoUtils.py index 3ca0f5d3..558f2fd4 100644 --- a/pylot/core/active/fmtomoUtils.py +++ b/pylot/core/active/fmtomoUtils.py @@ -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), From 734eca30dbde13f57ef9997d6517ee039b693e89 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 15:49:34 +0100 Subject: [PATCH 2/8] Implementation of new function iteratepicker of module autopick. --- autoPyLoT.py | 55 +++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 50 insertions(+), 5 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 85d65a48..fcfca838 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -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,30 @@ 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: + print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \ + % badpicks + # 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 +183,32 @@ 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: + print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \ + % badpicks + # 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 From 07bbc2926e80d4b7855312fab143c5c1e14e7d89 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 15:51:22 +0100 Subject: [PATCH 3/8] New function iteratepicker for iterative picking. --- pylot/core/pick/autopick.py | 79 +++++++++++++++++++++++++++++++++++-- 1 file changed, 75 insertions(+), 4 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 988c8e5c..7685596f 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -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,74 @@ 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 + ''' + + 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) + pickparameter.setParam(iplot=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 + + + + + + + From 8197d8f3d57cdff4fafb839605434c9a5f9f8510 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 15:52:14 +0100 Subject: [PATCH 4/8] Some input parameters changed from integers to floats. --- autoPyLoT_local.in | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index baabf0ce..93e56f16 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -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] From 67ac5807788d244e52be2efe9fa71c5b0d9c1396 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 15:53:25 +0100 Subject: [PATCH 5/8] Suppressed print output in setParam. --- pylot/core/read/inputs.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/read/inputs.py b/pylot/core/read/inputs.py index e6d609e3..cd66e845 100644 --- a/pylot/core/read/inputs.py +++ b/pylot/core/read/inputs.py @@ -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): From 8a16643bd81926bae6649389141403fec6fd0874 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 16:02:25 +0100 Subject: [PATCH 6/8] Marginal changes. --- pylot/core/pick/autopick.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 7685596f..be31c183 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -815,6 +815,10 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter): :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: @@ -842,7 +846,6 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter): (3 * pickparameter.getParam('tlta'))) pickparameter.setParam(pickwinP=pickparameter.getParam('pickwinP') / 2) pickparameter.setParam(Precalcwin=pickparameter.getParam('Precalcwin') / 2) - pickparameter.setParam(iplot=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'))) From 8173e44b0497e0c96e22642394027f59ee9f5e5e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 20 Nov 2015 16:02:54 +0100 Subject: [PATCH 7/8] Matginal changes. --- autoPyLoT.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index fcfca838..6a89a660 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -131,8 +131,6 @@ def autoPyLoT(inputfile): if len(badpicks) == 0: print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") else: - print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \ - % badpicks # 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 @@ -195,8 +193,6 @@ def autoPyLoT(inputfile): if len(badpicks) == 0: print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") else: - print("autoPyLoT: Found bad onsets at station(s) %s, starting re-picking them ...") \ - % badpicks # 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 From bd0d96c2ffd25b7ab8cfee05c84d9aea07208ba5 Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Mon, 23 Nov 2015 11:35:15 +0100 Subject: [PATCH 8/8] changed input for generateFMTOMOinpu --- pylot/core/active/seismicArrayPreparation.py | 32 ++++++++++++++++---- 1 file changed, 26 insertions(+), 6 deletions(-) diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 47f50835..607d904f 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -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)