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

[fix] reformatted code and fixed magnitude_type bug

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2016-09-29 12:44:37 +02:00
parent dfefd8af87
commit 2e840cdfeb
1 changed files with 37 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -14,7 +14,8 @@ 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, select_for_phase from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \
select_for_phase
from pylot.core.util.utils import common_range, fit_curve from pylot.core.util.utils import common_range, fit_curve
@ -40,7 +41,8 @@ class Magnitude(object):
self._magnitudes = dict() self._magnitudes = dict()
def __str__(self): def __str__(self):
print('number of stations used: {0}\n'.format(len(self.magnitudes.values()))) print(
'number of stations used: {0}\n'.format(len(self.magnitudes.values())))
print('\tstation\tmagnitude') print('\tstation\tmagnitude')
for s, m in self.magnitudes.items(): print('\t{0}\t{1}'.format(s, m)) for s, m in self.magnitudes.items(): print('\t{0}\t{1}'.format(s, m))
@ -120,10 +122,12 @@ class Magnitude(object):
# Magnitude object # Magnitude object
# mag_error => weights (magnitude error estimate from peak_to_peak, calcsourcespec?) # mag_error => weights (magnitude error estimate from peak_to_peak, calcsourcespec?)
# weights => StationMagnitdeContribution # weights => StationMagnitdeContribution
mag = ope.Magnitude(mag=np.median([M.mag for M in self.magnitudes.values()]), mag = ope.Magnitude(
type=self.type, origin_id=self.origin_id, mag=np.median([M.mag for M in self.magnitudes.values()]),
station_count=len(self.magnitudes), magnitude_type=self.type,
azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap) origin_id=self.origin_id,
station_count=len(self.magnitudes),
azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap)
return mag return mag
return None return None
@ -206,8 +210,9 @@ class RichterMagnitude(Magnitude):
plt.plot(th, sqH) plt.plot(th, sqH)
plt.plot(th[iwin], sqH[iwin], 'g') plt.plot(th[iwin], sqH[iwin], 'g')
plt.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2) plt.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2)
plt.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \ plt.title(
% (st[0].stats.station, wapp)) 'Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \
% (st[0].stats.station, wapp))
plt.xlabel('Time [s]') plt.xlabel('Time [s]')
plt.ylabel('Displacement [mm]') plt.ylabel('Displacement [mm]')
plt.show() plt.show()
@ -226,7 +231,8 @@ class RichterMagnitude(Magnitude):
station=station), a.phase) station=station), a.phase)
if not wf: if not wf:
if self.verbose: 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
delta = degrees2kilometers(a.distance) delta = degrees2kilometers(a.distance)
@ -244,7 +250,8 @@ class RichterMagnitude(Magnitude):
# using standard Gutenberg-Richter relation # using standard Gutenberg-Richter relation
# TODO make the ML calculation more flexible by allowing # TODO make the ML calculation more flexible by allowing
# use of custom relation functions # use of custom relation functions
magnitude = ope.StationMagnitude(mag=np.log10(a0) + richter_magnitude_scaling(delta)) magnitude = ope.StationMagnitude(
mag=np.log10(a0) + richter_magnitude_scaling(delta))
magnitude.origin_id = self.origin_id magnitude.origin_id = self.origin_id
magnitude.waveform_id = pick.waveform_id magnitude.waveform_id = pick.waveform_id
magnitude.amplitude_id = amplitude.resource_id magnitude.amplitude_id = amplitude.resource_id
@ -265,7 +272,8 @@ class MomentMagnitude(Magnitude):
_props = dict() _props = dict()
def __init__(self, stream, event, vp, Qp, density, verbosity=False, iplot=False): def __init__(self, stream, event, vp, Qp, density, verbosity=False,
iplot=False):
super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot) super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot)
self._vp = vp self._vp = vp
@ -317,14 +325,17 @@ class MomentMagnitude(Magnitude):
distance = degrees2kilometers(a.distance) distance = degrees2kilometers(a.distance)
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,
incidence, self.p_attenuation, self.plot_flag, self.verbose) azimuth,
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: 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, self.verbose) m0, mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity,
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)
magnitude.origin_id = self.origin_id magnitude.origin_id = self.origin_id
@ -359,10 +370,13 @@ def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
delta = delta * 1000 # hypocentral distance in [m] delta = delta * 1000 # hypocentral distance in [m]
if verbosity: if verbosity:
print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station {0} ...".format(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)
freesurf = 2.0 # free surface correction, assuming vertical incidence freesurf = 2.0 # free surface correction, assuming vertical incidence
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)
@ -371,7 +385,9 @@ def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
Mw = np.log10(Mo) * 2 / 3 - 6.7 # for metric units Mw = np.log10(Mo) * 2 / 3 - 6.7 # for metric units
if verbosity: if verbosity:
print("calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(Mo, Mw)) print(
"calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(
Mo, Mw))
return Mo, Mw return Mo, Mw
@ -525,7 +541,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
w0 = np.median([w01, w02]) w0 = np.median([w01, w02])
fc = np.median([fc1, fc2]) fc = np.median([fc1, fc2])
if verbosity: 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:
f1 = plt.figure() f1 = plt.figure()
@ -644,7 +661,8 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
fc = fc0 fc = fc0
w0 = max(S) w0 = max(S)
if verbosity: if verbosity:
print("fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0)) print(
"fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
if iplot > 1: if iplot > 1:
plt.figure(iplot) plt.figure(iplot)