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 
Conflicts:
	pylot/core/active/activeSeismoPick.py
This commit is contained in:
Marcel Paffrath 2016-01-19 10:29:45 +01:00
commit 7edee1266a
5 changed files with 109 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
autoPyLoT.py
View File

@ -107,7 +107,7 @@ def autoPyLoT(inputfile):
########################################################## ##########################################################
# !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:
@ -199,10 +199,13 @@ def autoPyLoT(inputfile):
# 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'):
if finalpicks.getpicdic() is not None: if finalpicks.getpicdic() is not None:
writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
else: else:
writephases(picks, 'HYPO71', hypo71file) writephases(picks, 'HYPO71', hypo71file)
else:
writephases(picks, 'HYPO71', hypo71file)
endsplash = '''------------------------------------------\n' endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n' -----Finished event %s!-----\n'
@ -222,7 +225,7 @@ def autoPyLoT(inputfile):
########################################################## ##########################################################
# !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:
@ -312,10 +315,13 @@ def autoPyLoT(inputfile):
# write phase files for various location routines # write phase files for various location routines
# HYPO71 # HYPO71
hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID')) hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.getParam('eventID'))
if hasattr(finalpicks, 'getpicdic'):
if finalpicks.getpicdic() is not None: if finalpicks.getpicdic() is not None:
writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
else: else:
writephases(picks, 'HYPO71', hypo71file) writephases(picks, 'HYPO71', hypo71file)
else:
writephases(picks, 'HYPO71', hypo71file)
endsplash = '''------------------------------------------\n' endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n' -----Finished event %s!-----\n'

14
autoPyLoT_local.in
View File

@ -6,8 +6,8 @@
#main settings# #main settings#
/DATA/Insheim #rootpath# %project path /DATA/Insheim #rootpath# %project path
EVENT_DATA/LOCAL #datapath# %data path EVENT_DATA/LOCAL #datapath# %data path
2013.02_Insheim #database# %name of data base 2015.11_Insheim #database# %name of data base
e0019.048.13 #eventID# %event ID for single event processing #eventID# %event ID for single event processing
/DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file /DATA/Insheim/STAT_INFO #invdir# %full path to inventory or dataless-seed file
PILOT #datastructure#%choose data structure PILOT #datastructure#%choose data structure
0 #iplot# %flag for plotting: 0 none, 1 partly, >1 everything 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 49 49.4 #blat# %lattitude bounding for location map
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
#common settings picker# #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 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 -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 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 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] 2 20 #bph2# %lower/upper corner freq. of second band pass filter z-comp. [Hz]
#special settings for calculating CF# #special settings for calculating CF#
%!!Be careful when editing the following!! %!!Edit the following only if you know what you are doing!!%
#Z-component# #Z-component#
HOS #algoP# %choose algorithm for P-onset determination (HOS, ARZ, or AR3) 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] 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 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.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] 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 #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.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] 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 1.5 #minAICSSNR# %below this SNR the initial S pick is rejected
#check duration of signal using envelope function# #check duration of signal using envelope function#
3 #minsiglength# %minimum required length of signal [s] 3 #minsiglength# %minimum required length of signal [s]
1.2 #noisefactor# %noiselevel*noisefactor=threshold 1.0 #noisefactor# %noiselevel*noisefactor=threshold
50 #minpercent# %required percentage of samples higher than threshold 40 #minpercent# %required percentage of samples higher than threshold
#check for spuriously picked S-onsets# #check for spuriously picked S-onsets#
2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude 2.0 #zfac# %P-amplitude must exceed at least zfac times RMS-S amplitude
#check statistics of P onsets# #check statistics of P onsets#

18
pylot/core/active/activeSeismoPick.py
View File

@ -147,6 +147,7 @@ class Survey(object):
labelm = 'manual picks' labelm = 'manual picks'
labela = 'automatic picks' labela = 'automatic picks'
fig = plt.figure() fig = plt.figure()
ax = fig.add_subplot(111) ax = fig.add_subplot(111)
@ -362,10 +363,19 @@ class Survey(object):
count += 1 count += 1
ttfile.close() ttfile.close()
srcfile.close() srcfile.close()
print 'Wrote output for %s traces' %count msg = 'Wrote output for {0} traces\n' \
print 'WARNING: output generated for FMTOMO-obsdata. Obsdata seems to take Lat, Lon, Depth and creates output for FMTOMO as Depth, Lat, Lon' 'WARNING: output generated for FMTOMO-obsdata. Obsdata seems ' \
print 'Dimensions of the seismic Array, transformed for FMTOMO, are Depth(%s, %s), Lat(%s, %s), Lon(%s, %s)'%( 'to take Lat, Lon, Depth and creates output for FMTOMO as ' \
min(DepthAll), max(DepthAll), min(LatAll), max(LatAll), min(LonAll), max(LonAll)) '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): def countPickedTraces(self, shot):
count = 0 count = 0

78
pylot/core/analysis/magnitude.py
View File

@ -1,7 +1,7 @@
#!/usr/bin/env python #!/usr/bin/env python
# -*- coding: utf-8 -*- # -*- coding: utf-8 -*-
""" """
Created August/September 2015. Created autumn/winter 2015.
:author: Ludger Küperkoch / MAGS2 EP3 working group :author: Ludger Küperkoch / MAGS2 EP3 working group
""" """
@ -18,7 +18,7 @@ from pylot.core.read.data import Data
class Magnitude(object): class Magnitude(object):
''' '''
Superclass for calculating Wood-Anderson peak-to-peak Superclass for calculating Wood-Anderson peak-to-peak
amplitudes, local magnitudes, seismic moments amplitudes, local magnitudes, source spectra, seismic moments
and moment magnitudes. and moment magnitudes.
''' '''
@ -52,7 +52,7 @@ class Magnitude(object):
:param: vp [m/s], P-velocity :param: vp [m/s], P-velocity
:param: integer :param: integer
:param: invdir, path to inventory or dataless-SEED file :param: invdir, name and path to inventory or dataless-SEED file
:type: string :type: string
''' '''
@ -240,8 +240,8 @@ class M0Mw(Magnitude):
# call subfunction to estimate source spectrum # call subfunction to estimate source spectrum
# and to derive w0 and fc # and to derive w0 and fc
[w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \ [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \
self.getinvdir(), az, inc, self.getQp(), \ self.getinvdir(), self.getvp(), delta, az, \
self.getiplot()) inc, self.getQp(), self.getiplot())
if w0 is not None: if w0 is not None:
# call subfunction to calculate Mo and Mw # 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 Subfunction of run_calcMoMw to calculate individual
seismic moments and corresponding moment magnitudes. 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] tr = wfstream[0]
delta = delta * 1000 # hypocentral distance in [m]
print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \
% tr.stats.station) % 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) 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)) 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 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 (usually called omega0) and the corner frequency assuming Aki's omega-square
source model. Has to be derived from instrument corrected displacement traces, source model. Has to be derived from instrument corrected displacement traces,
thus restitution and integration necessary! Integrated traces have to be rotated thus restitution and integration necessary! Integrated traces are rotated
into ray-coordinate system ZNE => LQT! 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 ....") print ("Calculating source spectrum ....")
@ -302,6 +346,7 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
Q = int(qu[0]) Q = int(qu[0])
# A, i.e. power of frequency # A, i.e. power of frequency
A = float(qu[1]) A = float(qu[1])
delta = delta * 1000 # hypocentral distance in [m]
fc = None fc = None
w0 = None w0 = None
@ -343,7 +388,11 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
LQT=cordat_copy.rotate('ZNE->LQT',azimuth, incidence) LQT=cordat_copy.rotate('ZNE->LQT',azimuth, incidence)
ldat = LQT.select(component="L") ldat = LQT.select(component="L")
if len(ldat) == 0: 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") ldat = LQT.select(component="Z")
# integrate to displacement # integrate to displacement
@ -407,8 +456,12 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
fi = np.where((f >= 1) & (f < 100)) fi = np.where((f >= 1) & (f < 100))
F = f[fi] F = f[fi]
YY = Y[fi] YY = Y[fi]
# correction for attenuation # 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 # get plateau (DC value) and corner frequency
# initial guess of plateau # initial guess of plateau
w0in = np.mean(YYcor[0:100]) w0in = np.mean(YYcor[0:100])
@ -428,7 +481,8 @@ def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, Qp, iplot):
# use of conventional fitting # use of conventional fitting
[w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot) [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]) w0 = np.median([w01, w02])
fc = np.median([fc1, fc2]) fc = np.median([fc1, fc2])
print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc))

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

@ -638,6 +638,11 @@ def autopickstation(wfstream, pickparam, verbose=False):
plt.legend([p1, p2], ['Data', 'CF1']) plt.legend([p1, p2], ['Data', 'CF1'])
plt.title('%s, P Weight=%d, SNR=None, ' plt.title('%s, P Weight=%d, SNR=None, '
'SNRdB=None' % (tr_filt.stats.channel, Pweight)) '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.yticks([])
plt.ylim([-1.5, 1.5]) plt.ylim([-1.5, 1.5])
plt.ylabel('Normalized Counts') 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) badpicks[i][1] = picks[badpicks[i][0]]['P']['mpp'] - float(res)
# get corresponding waveform stream # get corresponding waveform stream
msg = '#######################################################\n' \
'iteratepicker: Re-picking station {0}'.format(badpicks[i][0])
print(msg)
wf2pick = wf.select(station=badpicks[i][0]) wf2pick = wf.select(station=badpicks[i][0])
# modify some picking parameters # modify some picking parameters
@ -836,8 +844,8 @@ def iteratepicker(wf, NLLocfile, picks, badpicks, pickparameter):
Precalcwin_old = pickparameter.getParam('Precalcwin') Precalcwin_old = pickparameter.getParam('Precalcwin')
noisefactor_old = pickparameter.getParam('noisefactor') noisefactor_old = pickparameter.getParam('noisefactor')
zfac_old = pickparameter.getParam('zfac') zfac_old = pickparameter.getParam('zfac')
pickparameter.setParam(pstart=badpicks[i][1] - wf2pick[0].stats.starttime \ pickparameter.setParam(pstart=max([0, badpicks[i][1] - wf2pick[0].stats.starttime \
- pickparameter.getParam('tlta')) - pickparameter.getParam('tlta')]))
pickparameter.setParam(pstop=pickparameter.getParam('pstart') + \ pickparameter.setParam(pstop=pickparameter.getParam('pstart') + \
(3 * pickparameter.getParam('tlta'))) (3 * pickparameter.getParam('tlta')))
pickparameter.setParam(sstop=pickparameter.getParam('sstop') / 2) pickparameter.setParam(sstop=pickparameter.getParam('sstop') / 2)