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

[change] improved verbosity and plotting control for Magnitude objects

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2016-09-28 10:56:05 +02:00
parent 72d15e1fc5
commit c1bddd5c0b
1 changed files with 49 additions and 31 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -31,9 +31,9 @@ class Magnitude(object):
Base class object for Magnitude calculation within PyLoT. Base class object for Magnitude calculation within PyLoT.
""" """
def __init__(self, stream, event, verbosity=False): def __init__(self, stream, event, verbosity=False, iplot=0):
self._type = "M" self._type = "M"
self._plot_flag = False self._plot_flag = iplot
self._verbosity = verbosity self._verbosity = verbosity
self._event = event self._event = event
self._stream = stream self._stream = stream
@ -59,6 +59,17 @@ class Magnitude(object):
def plot_flag(self, value): def plot_flag(self, value):
self._plot_flag = value self._plot_flag = value
@property
def verbose(self):
return self._verbosity
@verbose.setter
def verbose(self, value):
if not isinstance(value, bool):
print('WARNING: only boolean values accepted...\n')
value = bool(value)
self._verbosity = value
@property @property
def stream(self): def stream(self):
return self._stream return self._stream
@ -115,8 +126,8 @@ class RichterMagnitude(Magnitude):
'sensitivity': 1 'sensitivity': 1
} }
def __init__(self, stream, event, calc_win, verbosity=False): def __init__(self, stream, event, calc_win, verbosity=False, iplot=0):
super(RichterMagnitude, self).__init__(stream, event, verbosity) super(RichterMagnitude, self).__init__(stream, event, verbosity, iplot)
self._calc_win = calc_win self._calc_win = calc_win
self._type = 'ML' self._type = 'ML'
@ -156,7 +167,7 @@ class RichterMagnitude(Magnitude):
# get maximum peak within pick window # get maximum peak within pick window
iwin = getsignalwin(th, t0 - stime, self.calc_win) iwin = getsignalwin(th, t0 - stime, self.calc_win)
wapp = np.max(sqH[iwin]) wapp = np.max(sqH[iwin])
if self._verbosity: if self.verbose:
print("Determined Wood-Anderson peak-to-peak amplitude: {0} " print("Determined Wood-Anderson peak-to-peak amplitude: {0} "
"mm".format(wapp)) "mm".format(wapp))
@ -186,6 +197,7 @@ class RichterMagnitude(Magnitude):
wf = select_for_phase(self.stream.select( wf = select_for_phase(self.stream.select(
station=station), a.phase) station=station), a.phase)
if not wf: if not wf:
if self.verbose:
print('WARNING: no waveform data found for station {0}'.format( print('WARNING: no waveform data found for station {0}'.format(
station)) station))
continue continue
@ -209,8 +221,8 @@ class MomentMagnitude(Magnitude):
corresponding moment magntiude Mw. corresponding moment magntiude Mw.
''' '''
def __init__(self, stream, event, vp, Qp, density, verbosity=False): def __init__(self, stream, event, vp, Qp, density, verbosity=False, iplot=False):
super(MomentMagnitude, self).__init__(stream, event) super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot)
self._vp = vp self._vp = vp
self._Qp = Qp self._Qp = Qp
@ -245,17 +257,18 @@ class MomentMagnitude(Magnitude):
azimuth = a.azimuth azimuth = a.azimuth
incidence = a.takeoff_angle incidence = a.takeoff_angle
w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance, azimuth, w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance, azimuth,
incidence, self.p_attenuation) incidence, self.p_attenuation, self.plot_flag, self.verbose)
if w0 is None or fc is None: if w0 is None or fc is None:
if self.verbose:
print("WARNING: insufficient frequency information") print("WARNING: insufficient frequency information")
continue continue
wf = select_for_phase(wf, "P") wf = select_for_phase(wf, "P")
M0, Mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, distance) M0, Mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, distance, self.verbose)
mag = dict(w0=w0, fc=fc, M0=M0, mag=Mw) mag = dict(w0=w0, fc=fc, M0=M0, mag=Mw)
self.magnitudes = (station, mag) self.magnitudes = (station, mag)
def calcMoMw(wfstream, w0, rho, vp, delta): def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
''' '''
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.
@ -279,8 +292,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta):
tr = wfstream[0] tr = wfstream[0]
delta = delta * 1000 # hypocentral distance in [m] delta = delta * 1000 # hypocentral distance in [m]
print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ if verbosity:
% tr.stats.station) print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station {0} ...".format(tr.stats.station))
# additional common parameters for calculating Mo # additional common parameters for calculating Mo
rP = 2 / np.sqrt(15) # average radiation pattern of P waves (Aki & Richards, 1980) rP = 2 / np.sqrt(15) # average radiation pattern of P waves (Aki & Richards, 1980)
@ -291,13 +304,14 @@ def calcMoMw(wfstream, w0, rho, vp, delta):
# 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 Mw = np.log10(Mo) * 2 / 3 - 6.7 # for metric units
print("calcMoMw: Calculated seismic moment Mo = %e Nm => Mw = %3.1f " % (Mo, Mw)) if verbosity:
print("calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(Mo, Mw))
return Mo, Mw return Mo, Mw
def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
qp, iplot=0): qp, iplot=0, verbosity=False):
''' '''
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
@ -329,6 +343,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
:param: iplot, show results (iplot>1) or not (iplot(<2) :param: iplot, show results (iplot>1) or not (iplot(<2)
:type: integer :type: integer
''' '''
if verbosity:
print ("Calculating source spectrum ....") print ("Calculating source spectrum ....")
# get Q value # get Q value
@ -359,6 +374,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
# if horizontal channels are 2 and 3 # if horizontal channels are 2 and 3
# no azimuth information is available and thus no # no azimuth information is available and thus no
# rotation is possible! # rotation is possible!
if verbosity:
print("calcsourcespec: Azimuth information is missing, " print("calcsourcespec: Azimuth information is missing, "
"no rotation of components possible!") "no rotation of components possible!")
ldat = LQT.select(component="Z") ldat = LQT.select(component="Z")
@ -381,8 +397,9 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
# 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:
if verbosity:
print ("calcsourcespec: Something is wrong with the waveform, " print ("calcsourcespec: Something is wrong with the waveform, "
"no zero crossings derived!") "no zero crossings derived!\n")
print ("No calculation of source spectrum possible!") print ("No calculation of source spectrum possible!")
plotflag = 0 plotflag = 0
else: else:
@ -430,16 +447,18 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
[optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin]) [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin])
w01 = optspecfit[0] w01 = optspecfit[0]
fc1 = optspecfit[1] fc1 = optspecfit[1]
if verbosity:
print ("calcsourcespec: 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
[w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot) [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot, verbosity)
# get w0 and fc as median of both # get w0 and fc as median of both
# source spectrum fits # source spectrum fits
w0 = np.median([w01, w02]) w0 = np.median([w01, w02])
fc = np.median([fc1, fc2]) fc = np.median([fc1, fc2])
if verbosity:
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))
if iplot > 1: if iplot > 1:
@ -504,7 +523,7 @@ def synthsourcespec(f, omega0, fcorner):
return ssp return ssp
def fitSourceModel(f, S, fc0, iplot): def fitSourceModel(f, S, fc0, iplot, verbosity=False):
''' '''
Calculates synthetic source spectrum by varying corner frequency fc. Calculates synthetic source spectrum by varying corner frequency fc.
Returns best approximated plateau omega0 and corner frequency, i.e. with least Returns best approximated plateau omega0 and corner frequency, i.e. with least
@ -558,9 +577,8 @@ def fitSourceModel(f, S, fc0, iplot):
elif len(STD) == 0: elif len(STD) == 0:
fc = fc0 fc = fc0
w0 = max(S) w0 = max(S)
if verbosity:
print("fitSourceModel: best fc: %fHz, best w0: %e m/Hz" \ print("fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
% (fc, w0))
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)