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-12-03 15:25:24 +01:00
commit 0d8324a4ae
3 changed files with 207 additions and 124 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

27
autoPyLoT.py
View File

@ -6,6 +6,7 @@ import argparse
import glob import glob
import subprocess import subprocess
import string import string
import numpy as np
from obspy.core import read, UTCDateTime from obspy.core import read, UTCDateTime
from pylot.core.read.data import Data from pylot.core.read.data import Data
from pylot.core.read.inputs import AutoPickParameter from pylot.core.read.inputs import AutoPickParameter
@ -161,6 +162,8 @@ def autoPyLoT(inputfile):
picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter) picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
# write phases to NLLoc-phase file # write phases to NLLoc-phase file
picksExport(picks, 'NLLoc', phasefile) picksExport(picks, 'NLLoc', phasefile)
# remove actual NLLoc-location file to keep only the last
os.remove(nllocfile)
# locate the event # locate the event
locate(nlloccall, ctrfile) locate(nlloccall, ctrfile)
print("autoPyLoT: Iteration No. %d finished." % nlloccounter) print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
@ -181,13 +184,23 @@ def autoPyLoT(inputfile):
finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \
nllocfile, picks, parameter.getParam('rho'), \ nllocfile, picks, parameter.getParam('rho'), \
parameter.getParam('vp'), parameter.getParam('invdir')) parameter.getParam('vp'), parameter.getParam('invdir'))
# get network moment magntiude
netMw = []
for key in finalpicks.getpicdic():
if finalpicks.getpicdic()[key]['P']['Mw'] is not None:
netMw.append(finalpicks.getpicdic()[key]['P']['Mw'])
netMw = np.median(netMw)
print("Network moment magnitude: %4.1f" % netMw)
else: else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!") print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
########################################################## ##########################################################
# write phase files for various location routines # write phase files for various location routines
# HYPO71 # HYPO71
hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, evID) hypo71file = '%s/autoPyLoT_HYPO71.pha' % event
if finalpicks.getpicdic() is not None:
writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
else:
writephases(picks, 'HYPO71', hypo71file)
endsplash = '''------------------------------------------\n' endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n' -----Finished event %s!-----\n'
@ -260,6 +273,8 @@ def autoPyLoT(inputfile):
picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter) picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter)
# write phases to NLLoc-phase file # write phases to NLLoc-phase file
picksExport(picks, 'NLLoc', phasefile) picksExport(picks, 'NLLoc', phasefile)
# remove actual NLLoc-location file to keep only the last
os.remove(nllocfile)
# locate the event # locate the event
locate(nlloccall, ctrfile) locate(nlloccall, ctrfile)
print("autoPyLoT: Iteration No. %d finished." % nlloccounter) print("autoPyLoT: Iteration No. %d finished." % nlloccounter)
@ -280,13 +295,23 @@ def autoPyLoT(inputfile):
finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \ finalpicks = M0Mw(wfdat, None, None, parameter.getParam('iplot'), \
nllocfile, picks, parameter.getParam('rho'), \ nllocfile, picks, parameter.getParam('rho'), \
parameter.getParam('vp'), parameter.getParam('invdir')) parameter.getParam('vp'), parameter.getParam('invdir'))
# get network moment magntiude
netMw = []
for key in finalpicks.getpicdic():
if finalpicks.getpicdic()[key]['P']['Mw'] is not None:
netMw.append(finalpicks.getpicdic()[key]['P']['Mw'])
netMw = np.median(netMw)
print("Network moment magnitude: %4.1f" % netMw)
else: else:
print("autoPyLoT: No NLLoc-location file available! Stop iteration!") print("autoPyLoT: No NLLoc-location file available! Stop iteration!")
########################################################## ##########################################################
# 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 finalpicks.getpicdic() is not None:
writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file) writephases(finalpicks.getpicdic(), 'HYPO71', hypo71file)
else:
writephases(picks, 'HYPO71', hypo71file)
endsplash = '''------------------------------------------\n' endsplash = '''------------------------------------------\n'
-----Finished event %s!-----\n' -----Finished event %s!-----\n'

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

@ -12,7 +12,7 @@ from obspy.core import Stream, UTCDateTime
from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all
from pylot.core.util.utils import getPatternLine from pylot.core.util.utils import getPatternLine
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
from scipy import integrate from scipy import integrate, signal
from pylot.core.read.data import Data from pylot.core.read.data import Data
class Magnitude(object): class Magnitude(object):
@ -200,7 +200,7 @@ class WApp(Magnitude):
class M0Mw(Magnitude): class M0Mw(Magnitude):
''' '''
Method to calculate seismic moment Mo and moment magnitude Mw. Method to calculate seismic moment Mo and moment magnitude Mw.
Requires results of class w0fc for calculating plateau w0 Requires results of class calcsourcespec for calculating plateau w0
and corner frequency fc of source spectrum, respectively. Uses and corner frequency fc of source spectrum, respectively. Uses
subfunction calcMoMw.py. Returns modified dictionary of picks including subfunction calcMoMw.py. Returns modified dictionary of picks including
Dc-value, corner frequency fc, seismic moment Mo and Dc-value, corner frequency fc, seismic moment Mo and
@ -212,17 +212,12 @@ class M0Mw(Magnitude):
picks = self.getpicks() picks = self.getpicks()
nllocfile = self.getNLLocfile() nllocfile = self.getNLLocfile()
wfdat = self.getwfstream() wfdat = self.getwfstream()
# get vertical component data only self.picdic = None
zdat = wfdat.select(component="Z")
if len(zdat) == 0: # check for other components
zdat = wfdat.select(component="3")
for key in picks: for key in picks:
if picks[key]['P']['weight'] < 4: if picks[key]['P']['weight'] < 4:
# select waveform # select waveform
selwf = zdat.select(station=key) selwf = wfdat.select(station=key)
# get hypocentral distance of station
# from NLLoc-location file
if len(key) > 4: if len(key) > 4:
Ppattern = '%s ? ? ? P' % key Ppattern = '%s ? ? ? P' % key
elif len(key) == 4: elif len(key) == 4:
@ -230,15 +225,22 @@ class M0Mw(Magnitude):
elif len(key) < 4: elif len(key) < 4:
Ppattern = '%s ? ? ? P' % key Ppattern = '%s ? ? ? P' % key
nllocline = getPatternLine(nllocfile, Ppattern) nllocline = getPatternLine(nllocfile, Ppattern)
# get hypocentral distance, station azimuth and
# angle of incidence from NLLoc-location file
delta = float(nllocline.split(None)[21]) delta = float(nllocline.split(None)[21])
az = float(nllocline.split(None)[22])
inc = float(nllocline.split(None)[24])
# 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.getiplot(), self.getinvdir()) self.getinvdir(), az, inc, 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
[Mo, Mw] = calcMoMw(selwf, w0, self.getrho(), self.getvp(), \ zdat = selwf.select(component="Z")
if len(zdat) == 0: # check for other components
zdat = selwf.select(component="3")
[Mo, Mw] = calcMoMw(zdat, w0, self.getrho(), self.getvp(), \
delta, self.getinvdir()) delta, self.getinvdir())
else: else:
Mo = None Mo = None
@ -276,70 +278,112 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv):
def calcsourcespec(wfstream, onset, iplot, inventory): def calcsourcespec(wfstream, onset, inventory, azimuth, incidence, 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! thus restitution and integration necessary! Integrated traces have to be rotated
into ray-coordinate system ZNE => LQT!
''' '''
print ("Calculating source spectrum ....") print ("Calculating source spectrum ....")
fc = None fc = None
w0 = None w0 = None
data = Data() data = Data()
z_copy = wfstream.copy() wf_copy = wfstream.copy()
[corzdat, restflag] = data.restituteWFData(inventory, z_copy)
[cordat, restflag] = data.restituteWFData(inventory, wf_copy)
if restflag == 1: if restflag == 1:
# integrate to displacment zdat = cordat.select(component="Z")
corintzdat = integrate.cumtrapz(corzdat[0], None, corzdat[0].stats.delta) if len(zdat) == 0:
z_copy[0].data = corintzdat zdat = cordat.select(component="3")
tr = z_copy[0] cordat_copy = cordat.copy()
# get equal time stamps and lengths of traces
# necessary for rotation of traces
tr0start = cordat_copy[0].stats.starttime
tr0start = tr0start.timestamp
tr0end = cordat_copy[0].stats.endtime
tr0end = tr0end.timestamp
tr1start = cordat_copy[1].stats.starttime
tr1start = tr1start.timestamp
tr1end = cordat_copy[1].stats.endtime
tr1end = tr1end.timestamp
tr2start = cordat_copy[2].stats.starttime
tr2start = tr2start.timestamp
tr2end = cordat_copy[0].stats.endtime
tr2end = tr2end.timestamp
trstart = UTCDateTime(max([tr0start, tr1start, tr2start]))
trend = UTCDateTime(min([tr0end, tr1end, tr2end]))
cordat_copy.trim(trstart, trend)
minlen = min([len(cordat_copy[0]), len(cordat_copy[1]), len(cordat_copy[2])])
cordat_copy[0].data = cordat_copy[0].data[0:minlen]
cordat_copy[1].data = cordat_copy[1].data[0:minlen]
cordat_copy[2].data = cordat_copy[2].data[0:minlen]
try:
# rotate into LQT (ray-coordindate-) system using Obspy's rotate
# L: P-wave direction
# Q: SV-wave direction
# T: SH-wave direction
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
ldat = LQT.select(component="Z")
# integrate to displacement
# unrotated vertical component (for copmarison)
inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None, \
zdat[0].stats.delta))
# rotated component Z => L
Ldat = signal.detrend(integrate.cumtrapz(ldat[0].data, None, \
ldat[0].stats.delta))
# get window after P pulse for # get window after P pulse for
# calculating source spectrum # calculating source spectrum
if tr.stats.sampling_rate <= 100: if zdat[0].stats.sampling_rate <= 100:
winzc = tr.stats.sampling_rate winzc = zdat[0].stats.sampling_rate
elif tr.stats.sampling_rate > 100 and \ elif zdat[0].stats.sampling_rate > 100 and \
tr.stats.sampling_rate <= 200: zdat[0].stats.sampling_rate <= 200:
winzc = 0.5 * tr.stats.sampling_rate winzc = 0.5 * zdat[0].stats.sampling_rate
elif tr.stats.sampling_rate > 200 and \ elif zdat[0].stats.sampling_rate > 200 and \
tr.stats.sampling_rate <= 400: zdat[0].stats.sampling_rate <= 400:
winzc = 0.2 * tr.stats.sampling_rate winzc = 0.2 * zdat[0].stats.sampling_rate
elif tr.stats.sampling_rate > 400: elif zdat[0].stats.sampling_rate > 400:
winzc = tr.stats.sampling_rate winzc = zdat[0].stats.sampling_rate
tstart = UTCDateTime(tr.stats.starttime) tstart = UTCDateTime(zdat[0].stats.starttime)
tonset = onset.timestamp -tstart.timestamp tonset = onset.timestamp -tstart.timestamp
impickP = tonset * tr.stats.sampling_rate impickP = tonset * zdat[0].stats.sampling_rate
wfzc = tr.data[impickP : impickP + winzc] wfzc = Ldat[impickP : impickP + winzc]
# get time array # get time array
t = np.arange(0, len(tr) * tr.stats.delta, tr.stats.delta) t = np.arange(0, len(inttrz) * zdat[0].stats.delta, \
zdat[0].stats.delta)
# calculate spectrum using only first cycles of # calculate spectrum using only first cycles of
# waveform after P onset! # waveform after P onset!
zc = crossings_nonzero_all(wfzc) zc = crossings_nonzero_all(wfzc)
if np.size(zc) == 0 or len(zc) <= 3: if np.size(zc) == 0 or len(zc) <= 3:
print ("Something is wrong with the waveform, " print ("calcsourcespec: Something is wrong with the waveform, "
"no zero crossings derived!") "no zero crossings derived!")
print ("No calculation of source spectrum possible!") print ("No calculation of source spectrum possible!")
plotflag = 0 plotflag = 0
else: else:
plotflag = 1 plotflag = 1
index = min([3, len(zc) - 1]) index = min([3, len(zc) - 1])
calcwin = (zc[index] - zc[0]) * z_copy[0].stats.delta calcwin = (zc[index] - zc[0]) * zdat[0].stats.delta
iwin = getsignalwin(t, tonset, calcwin) iwin = getsignalwin(t, tonset, calcwin)
xdat = tr.data[iwin] xdat = Ldat[iwin]
# fft # fft
fny = tr.stats.sampling_rate / 2 fny = zdat[0].stats.sampling_rate / 2
l = len(xdat) / tr.stats.sampling_rate l = len(xdat) / zdat[0].stats.sampling_rate
n = tr.stats.sampling_rate * l # number of fft bins after Bath # number of fft bins after Bath
n = zdat[0].stats.sampling_rate * l
# find next power of 2 of data length # find next power of 2 of data length
m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2))) m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
N = int(np.power(m, 2)) N = int(np.power(m, 2))
y = tr.stats.delta * np.fft.fft(xdat, N) y = zdat[0].stats.delta * np.fft.fft(xdat, N)
Y = abs(y[: N/2]) Y = abs(y[: N/2])
L = (N - 1) / tr.stats.sampling_rate L = (N - 1) / zdat[0].stats.sampling_rate
f = np.arange(0, fny, 1/L) f = np.arange(0, fny, 1/L)
# remove zero-frequency and frequencies above # remove zero-frequency and frequencies above
@ -361,7 +405,7 @@ def calcsourcespec(wfstream, onset, iplot, inventory):
[optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [w0in, Fcin]) [optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [w0in, Fcin])
w01 = optspecfit[0] w01 = optspecfit[0]
fc1 = optspecfit[1] fc1 = optspecfit[1]
print ("w0fc: Determined w0-value: %e m/Hz, \n" print ("calcsourcespec: Determined w0-value: %e m/Hz, \n"
"Determined corner frequency: %f Hz" % (w01, fc1)) "Determined corner frequency: %f Hz" % (w01, fc1))
# use of conventional fitting # use of conventional fitting
@ -370,21 +414,28 @@ def calcsourcespec(wfstream, onset, iplot, inventory):
# get w0 and fc as median # get w0 and fc as median
w0 = np.median([w01, w02]) w0 = np.median([w01, w02])
fc = np.median([fc1, fc2]) fc = np.median([fc1, fc2])
print("w0fc: 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))
except TypeError as er:
raise TypeError('''{0}'''.format(er))
if iplot > 1: if iplot > 1:
f1 = plt.figure() f1 = plt.figure()
tLdat = np.arange(0, len(Ldat) * zdat[0].stats.delta, \
zdat[0].stats.delta)
plt.subplot(2,1,1) plt.subplot(2,1,1)
# show displacement in mm # show displacement in mm
plt.plot(t, np.multiply(tr, 1000), 'k') p1, = plt.plot(t, np.multiply(inttrz, 1000), 'k')
p2, = plt.plot(tLdat, np.multiply(Ldat, 1000))
plt.legend([p1, p2], ['Displacement', 'Rotated Displacement'])
if plotflag == 1: if plotflag == 1:
plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') plt.plot(t[iwin], np.multiply(xdat, 1000), 'g')
plt.title('Seismogram and P Pulse, Station %s-%s' \ plt.title('Seismogram and P Pulse, Station %s-%s' \
% (tr.stats.station, tr.stats.channel)) % (zdat[0].stats.station, zdat[0].stats.channel))
else: else:
plt.title('Seismogram, Station %s-%s' \ plt.title('Seismogram, Station %s-%s' \
% (tr.stats.station, tr.stats.channel)) % (zdat[0].stats.station, zdat[0].stats.channel))
plt.xlabel('Time since %s' % tr.stats.starttime) plt.xlabel('Time since %s' % zdat[0].stats.starttime)
plt.ylabel('Displacement [mm]') plt.ylabel('Displacement [mm]')
if plotflag == 1: if plotflag == 1:
@ -474,8 +525,13 @@ def fitSourceModel(f, S, fc0, iplot):
STD.append(stddc + stdFC) STD.append(stddc + stdFC)
# get best found w0 anf fc from minimum # get best found w0 anf fc from minimum
if len(STD) > 0:
fc = fc[np.argmin(STD)] fc = fc[np.argmin(STD)]
w0 = w0[np.argmin(STD)] w0 = w0[np.argmin(STD)]
elif len(STD) == 0:
fc = fc0
w0 = max(S)
print("fitSourceModel: best fc: %fHz, best w0: %e m/Hz" \ print("fitSourceModel: best fc: %fHz, best w0: %e m/Hz" \
% (fc, w0)) % (fc, w0))

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

@ -204,7 +204,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
z_copy[0].data = tr_filt.data z_copy[0].data = tr_filt.data
zne = z_copy zne = z_copy
if len(ndat) == 0 or len(edat) == 0: if len(ndat) == 0 or len(edat) == 0:
msg = 'One or more horizontal components missing!\nSignal ' \ msg = 'One or more horizontal component(s) missing!\nSignal ' \
'length only checked on vertical component!\n' \ 'length only checked on vertical component!\n' \
'Decreasing minsiglengh from {0} to ' \ 'Decreasing minsiglengh from {0} to ' \
'{1}'.format(minsiglength, minsiglength / 2) '{1}'.format(minsiglength, minsiglength / 2)
@ -424,6 +424,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
'SNR: {1}\nGo on with refined picking ...\n' \ 'SNR: {1}\nGo on with refined picking ...\n' \
'autopickstation: re-filtering horizontal traces ' \ 'autopickstation: re-filtering horizontal traces ' \
'...'.format(aicarhpick.getSlope(), aicarhpick.getSNR()) '...'.format(aicarhpick.getSlope(), aicarhpick.getSNR())
if verbose: print(msg)
# re-calculate CF from re-filtered trace in vicinity of initial # re-calculate CF from re-filtered trace in vicinity of initial
# onset # onset
cuttimesh2 = [round(aicarhpick.getpick() - Srecalcwin), cuttimesh2 = [round(aicarhpick.getpick() - Srecalcwin),
@ -774,7 +775,8 @@ def autopickstation(wfstream, pickparam, verbose=False):
# for P phase # for P phase
phase = 'P' phase = 'P'
phasepick = {'lpp': lpickP, 'epp': epickP, 'mpp': mpickP, 'spe': Perror, phasepick = {'lpp': lpickP, 'epp': epickP, 'mpp': mpickP, 'spe': Perror,
'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM} 'snr': SNRP, 'snrdb': SNRPdB, 'weight': Pweight, 'fm': FM,
'w0': None, 'fc': None, 'Mo': None, 'Mw': None}
picks = {phase: phasepick} picks = {phase: phasepick}
# add P marker # add P marker
picks[phase]['marked'] = Pmarker picks[phase]['marked'] = Pmarker