Bugfix in class MomentMagnitude: calcsourcespec was feeded with 1 component only, thus rotation of traces was no more possible.

This commit is contained in:
Ludger Küperkoch 2016-10-27 09:07:08 +02:00
parent c055b52325
commit bcde390d42

View File

@ -6,14 +6,12 @@ Created autumn/winter 2015.
:author: Ludger Küperkoch / MAGS2 EP3 working group :author: Ludger Küperkoch / MAGS2 EP3 working group
""" """
import os import os
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
import obspy.core.event as ope import obspy.core.event as ope
from obspy.geodetics import degrees2kilometers from obspy.geodetics import degrees2kilometers
from scipy import integrate, signal from scipy import integrate, signal
from scipy.optimize import curve_fit from scipy.optimize import curve_fit
from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \ from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \
select_for_phase select_for_phase
from pylot.core.util.utils import common_range, fit_curve from pylot.core.util.utils import common_range, fit_curve
@ -317,8 +315,8 @@ class MomentMagnitude(Magnitude):
continue continue
pick = a.pick_id.get_referred_object() pick = a.pick_id.get_referred_object()
station = pick.waveform_id.station_code station = pick.waveform_id.station_code
wf = select_for_phase(self.stream.select( scopy = self.stream.copy()
station=station), a.phase) wf = scopy.select(station=station)
if not wf: if not wf:
continue continue
onset = pick.time onset = pick.time
@ -326,15 +324,16 @@ 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, w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance,
azimuth, azimuth, incidence, self.p_attenuation,
incidence, self.p_attenuation,
self.plot_flag, self.verbose) self.plot_flag, self.verbose)
if w0 is None or fc is None: if w0 is None or fc is None:
if self.verbose: 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(self.stream.select(
m0, mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, station=station), a.phase)
WF = select_for_phase(WF, "P")
m0, mw = calcMoMw(WF, w0, self.rock_density, self.p_velocity,
distance, self.verbose) distance, self.verbose)
self.moment_props = (station, dict(w0=w0, fc=fc, Mo=m0)) self.moment_props = (station, dict(w0=w0, fc=fc, Mo=m0))
magnitude = ope.StationMagnitude(mag=mw) magnitude = ope.StationMagnitude(mag=mw)
@ -426,7 +425,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
:type: integer :type: integer
''' '''
if verbosity: if verbosity:
print ("Calculating source spectrum ....") print ("Calculating source spectrum for station %s ...." % wfstream[0].stats.station)
# get Q value # get Q value
Q, A = qp Q, A = qp