diff --git a/autoPyLoT.py b/autoPyLoT.py index 658739fb..00068c73 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -107,7 +107,7 @@ def autoPyLoT(inputfile): ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) - + finalpicks = picks ########################################################## # locating if locflag == 1: @@ -199,8 +199,11 @@ def autoPyLoT(inputfile): # write phase files for various location routines # HYPO71 hypo71file = '%s/autoPyLoT_HYPO71.pha' % event - if finalpicks.getpicdic() is not None: - writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) + if hasattr(finalpicks, 'getpicdic'): + if finalpicks.getpicdic() is not None: + writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) + else: + writephases(picks, 'HYPO71', hypo71file) else: writephases(picks, 'HYPO71', hypo71file) @@ -222,7 +225,7 @@ def autoPyLoT(inputfile): ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) - + finalpicks = picks ########################################################## # locating if locflag == 1: @@ -312,8 +315,11 @@ def autoPyLoT(inputfile): # write phase files for various location routines # HYPO71 hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID')) - if finalpicks.getpicdic() is not None: - writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) + if hasattr(finalpicks, 'getpicdic'): + if finalpicks.getpicdic() is not None: + writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) + else: + writephases(picks, 'HYPO71', hypo71file) else: writephases(picks, 'HYPO71', hypo71file) diff --git a/autoPyLoT_local.in b/autoPyLoT_local.in index 1b845cbc..329a0fe2 100644 --- a/autoPyLoT_local.in +++ b/autoPyLoT_local.in @@ -6,8 +6,8 @@ #main settings# /DATA/Insheim #rootpath# %project path EVENT_DATA/LOCAL #datapath# %data path -2013.02_Insheim #database# %name of data base -e0019.048.13 #eventID# %event ID for single event processing +2015.11_Insheim #database# %name of data base + #eventID# %event ID for single event processing /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file PILOT #datastructure#%choose data structure 0 #iplot# %flag for plotting: 0 none, 1 partly, >1 everything @@ -38,7 +38,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file co 49 49.4 #blat# %lattitude bounding for location map %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% #common settings picker# -5.0 #pstart# %start time [s] for calculating CF for P-picking +15.0 #pstart# %start time [s] for calculating CF for P-picking 60.0 #pstop# %end time [s] for calculating CF for P-picking -1.0 #sstart# %start time [s] relative to P-onset for calculating CF for S-picking 10.0 #sstop# %end time [s] after P-onset for calculating CF for S-picking @@ -47,7 +47,7 @@ AUTOFOCMEC_AIC_HOS4_ARH.in #focmecin# %name of focmec input file co 2 15 #bph1# %lower/upper corner freq. of first band pass filter H-comp. [Hz] 2 20 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz] #special settings for calculating CF# -%!!Be careful when editing the following!! +%!!Edit the following only if you know what you are doing!!% #Z-component# HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3) 7.0 #tlta# %for HOS-/AR-AIC-picker, length of LTA window [s] @@ -60,7 +60,7 @@ HOS #algoP# %choose algorithm for P-onset 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.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) +6.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] @@ -94,8 +94,8 @@ ARH #algoS# %choose algorithm for S-onset 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected #check duration of signal using envelope function# 3 #minsiglength# %minimum required length of signal [s] -1.2 #noisefactor# %noiselevel*noisefactor=threshold -50 #minpercent# %required percentage of samples higher than threshold +1.0 #noisefactor# %noiselevel*noisefactor=threshold +40 #minpercent# %required percentage of samples higher than threshold #check for spuriously picked S-onsets# 2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude #check statistics of P onsets# diff --git a/pylot/core/active/activeSeismoPick.py b/pylot/core/active/activeSeismoPick.py index 5ae1cfa9..5632664f 100644 --- a/pylot/core/active/activeSeismoPick.py +++ b/pylot/core/active/activeSeismoPick.py @@ -147,6 +147,7 @@ class Survey(object): labelm = 'manual picks' labela = 'automatic picks' + fig = plt.figure() ax = fig.add_subplot(111) @@ -362,10 +363,19 @@ class Survey(object): count += 1 ttfile.close() srcfile.close() - print 'Wrote output for %s traces' %count - print 'WARNING: output generated for FMTOMO-obsdata. Obsdata seems to take Lat, Lon, Depth and creates output for FMTOMO as Depth, Lat, Lon' - print 'Dimensions of the seismic Array, transformed for FMTOMO, are Depth(%s, %s), Lat(%s, %s), Lon(%s, %s)'%( - min(DepthAll), max(DepthAll), min(LatAll), max(LatAll), min(LonAll), max(LonAll)) + msg = 'Wrote output for {0} traces\n' \ + 'WARNING: output generated for FMTOMO-obsdata. Obsdata seems ' \ + 'to take Lat, Lon, Depth and creates output for FMTOMO as ' \ + 'Depth, Lat, Lon\nDimensions of the seismic Array, ' \ + 'transformed for FMTOMO, are Depth({1}, {2}), Lat({3}, {4}), ' \ + 'Lon({5}, {6})'.format(count, + min(DepthAll), + max(DepthAll), + min(LatAll), + max(LatAll), + min(LonAll), + max(LonAll)) + print(msg) def countPickedTraces(self, shot): count = 0 diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index e6233b25..827d6055 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -1,7 +1,7 @@ #!/usr/bin/env python # -*- coding: utf-8 -*- """ -Created August/September 2015. +Created autumn/winter 2015. :author: Ludger Küperkoch / MAGS2 EP3 working group """ @@ -18,7 +18,7 @@ from pylot.core.read.data import Data class Magnitude(object): ''' Superclass for calculating Wood-Anderson peak-to-peak - amplitudes, local magnitudes, seismic moments + amplitudes, local magnitudes, source spectra, seismic moments and moment magnitudes. ''' @@ -52,7 +52,7 @@ class Magnitude(object): :param: vp [m/s], P-velocity :param: integer - :param: invdir, path to inventory or dataless-SEED file + :param: invdir, name and path to inventory or dataless-SEED file :type: string ''' @@ -240,8 +240,8 @@ class M0Mw(Magnitude): # call subfunction to estimate source spectrum # and to derive w0 and fc [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \ - self.getinvdir(), az, inc, self.getQp(), \ - self.getiplot()) + self.getinvdir(), self.getvp(), delta, az, \ + inc, self.getQp(), self.getiplot()) if w0 is not None: # call subfunction to calculate Mo and Mw @@ -265,9 +265,25 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv): ''' Subfunction of run_calcMoMw to calculate individual seismic moments and corresponding moment magnitudes. + + :param: wfstream + :type: `~obspy.core.stream.Stream` + + :param: w0, height of plateau of source spectrum + :type: float + + :param: rho, rock density [kg/m³] + :type: integer + + :param: delta, hypocentral distance [km] + :type: integer + + :param: inv, name/path of inventory or dataless-SEED file + :type: string ''' tr = wfstream[0] + delta = delta * 1000 # hypocentral distance in [m] print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ % tr.stats.station) @@ -278,7 +294,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv): Mo = w0 * 4 * np.pi * rho * np.power(vp, 3) * delta / (rP * freesurf) - Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 #after Hanks & Kanamori (1979), defined for [dyn*cm]! + #Mw = np.log10(Mo * 1e07) * 2 / 3 - 10.7 # after Hanks & Kanamori (1979), defined for [dyn*cm]! + Mw = np.log10(Mo) * 2 / 3 - 6.7 # for metric units print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw)) @@ -286,13 +303,40 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv): -def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): +def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, Qp, iplot): ''' Subfunction to calculate the source spectrum and to derive from that the plateau (usually called omega0) and the corner frequency assuming Aki's omega-square source model. Has to be derived from instrument corrected displacement traces, - thus restitution and integration necessary! Integrated traces have to be rotated - into ray-coordinate system ZNE => LQT! + thus restitution and integration necessary! Integrated traces are rotated + into ray-coordinate system ZNE => LQT using Obspy's rotate modul! + + :param: wfstream + :type: `~obspy.core.stream.Stream` + + :param: onset, P-phase onset time + :type: float + + :param: inventory, path/name of inventory or dataless-SEED file + :type: string + + :param: vp, Vp-wave velocity + :type: float + + :param: delta, hypocentral distance [km] + :type: integer + + :param: azimuth + :type: integer + + :param: incidence + :type: integer + + :param: Qp, quality factor for P-waves + :type: integer + + :param: iplot, show results (iplot>1) or not (iplot(<2) + :type: integer ''' print ("Calculating source spectrum ....") @@ -302,6 +346,7 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): Q = int(qu[0]) # A, i.e. power of frequency A = float(qu[1]) + delta = delta * 1000 # hypocentral distance in [m] fc = None w0 = None @@ -343,7 +388,11 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): LQT=cordat_copy.rotate('ZNE->LQT',azimuth, incidence) ldat = LQT.select(component="L") if len(ldat) == 0: - # yet Obspy's rotate can not handle channels 3/2/1 + # if horizontal channels are 2 and 3 + # no azimuth information is available and thus no + # rotation is possible! + print("calcsourcespec: Azimuth information is missing, " + "no rotation of components possible!") ldat = LQT.select(component="Z") # integrate to displacement @@ -407,8 +456,12 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): fi = np.where((f >= 1) & (f < 100)) F = f[fi] YY = Y[fi] + # correction for attenuation - YYcor = YY.real*Q**A + wa = 2 * np.pi * F #angular frequency + D = np.exp((wa * delta) / (2 * vp * Q*F**A)) + YYcor = YY.real*D + # get plateau (DC value) and corner frequency # initial guess of plateau w0in = np.mean(YYcor[0:100]) @@ -428,7 +481,8 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot): # use of conventional fitting [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot) - # get w0 and fc as median + # get w0 and fc as median of both + # source spectrum fits w0 = np.median([w01, w02]) fc = np.median([fc1, fc2]) print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index a54a94c2..fd86f600 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -638,6 +638,11 @@ def autopickstation(wfstream, pickparam, verbose=False): plt.legend([p1, p2], ['Data', 'CF1']) plt.title('%s, P Weight=%d, SNR=None, ' 'SNRdB=None' % (tr_filt.stats.channel, Pweight)) + else: + plt.title('%s, %s, P Weight=%d' % (tr_filt.stats.station, + tr_filt.stats.channel, + Pweight)) + plt.yticks([]) plt.ylim([-1.5, 1.5]) plt.ylabel('Normalized Counts') @@ -826,6 +831,9 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter): badpicks[i][1] = picks[badpicks[i][0]]['P']['mpp'] - float(res) # get corresponding waveform stream + msg = '#######################################################\n' \ + 'iteratepicker: Re-picking station {0}'.format(badpicks[i][0]) + print(msg) wf2pick = wf.select(station=badpicks[i][0]) # modify some picking parameters @@ -836,8 +844,8 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter): Precalcwin_old = pickparameter.getParam('Precalcwin') noisefactor_old = pickparameter.getParam('noisefactor') zfac_old = pickparameter.getParam('zfac') - pickparameter.setParam(pstart=badpicks[i][1] - wf2pick[0].stats.starttime \ - - pickparameter.getParam('tlta')) + pickparameter.setParam(pstart=max([0, badpicks[i][1] - wf2pick[0].stats.starttime \ + - pickparameter.getParam('tlta')])) pickparameter.setParam(pstop=pickparameter.getParam('pstart') + \ (3 * pickparameter.getParam('tlta'))) pickparameter.setParam(sstop=pickparameter.getParam('sstop') / 2)