[reformat] code reformatting with PyCharm
This commit is contained in:
@@ -6,27 +6,27 @@ Revised/extended summer 2017.
|
||||
|
||||
:author: Ludger Küperkoch / MAGS2 EP3 working group
|
||||
"""
|
||||
import os
|
||||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
import obspy.core.event as ope
|
||||
from obspy.geodetics import degrees2kilometers
|
||||
from scipy import integrate, signal
|
||||
from scipy.optimize import curve_fit
|
||||
from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, \
|
||||
select_for_phase
|
||||
from pylot.core.util.utils import common_range, fit_curve
|
||||
from scipy import integrate, signal
|
||||
from scipy.optimize import curve_fit
|
||||
|
||||
|
||||
def richter_magnitude_scaling(delta):
|
||||
distance = np.array([0, 10, 20, 25, 30, 35,40, 45, 50, 60, 70, 75, 85, 90, 100, 110,
|
||||
distance = np.array([0, 10, 20, 25, 30, 35, 40, 45, 50, 60, 70, 75, 85, 90, 100, 110,
|
||||
120, 130, 140, 150, 160, 170, 180, 190, 200, 210, 230, 240, 250,
|
||||
260, 270, 280, 290, 300, 310, 320, 330, 340, 350, 360, 370, 380,
|
||||
390, 400, 430, 470, 510, 560, 600, 700, 800, 900, 1000])
|
||||
richter_scaling = np.array([1.4, 1.5, 1.7, 1.9, 2.1, 2.3, 2.4, 2.5, 2.6, 2.8, 2.8, 2.9,
|
||||
2.9, 3.0, 3.1, 3.1, 3.2, 3.2, 3.3, 3.3, 3.4, 3.4, 3.5, 3.5,
|
||||
3.6, 3.7, 3.7, 3.8, 3.8, 3.9, 3.9, 4.0, 4.0, 4.1, 4.2, 4.2,
|
||||
4.2, 4.2, 4.3, 4.3, 4.3, 4.4, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9,
|
||||
5.1, 5.2, 5.4, 5.5, 5.7])
|
||||
richter_scaling = np.array([1.4, 1.5, 1.7, 1.9, 2.1, 2.3, 2.4, 2.5, 2.6, 2.8, 2.8, 2.9,
|
||||
2.9, 3.0, 3.1, 3.1, 3.2, 3.2, 3.3, 3.3, 3.4, 3.4, 3.5, 3.5,
|
||||
3.6, 3.7, 3.7, 3.8, 3.8, 3.9, 3.9, 4.0, 4.0, 4.1, 4.2, 4.2,
|
||||
4.2, 4.2, 4.3, 4.3, 4.3, 4.4, 4.4, 4.5, 4.6, 4.7, 4.8, 4.9,
|
||||
5.1, 5.2, 5.4, 5.5, 5.7])
|
||||
# prepare spline interpolation to calculate return value
|
||||
func, params = fit_curve(distance, richter_scaling)
|
||||
return func(delta, params)
|
||||
@@ -47,7 +47,7 @@ class Magnitude(object):
|
||||
|
||||
def __str__(self):
|
||||
print(
|
||||
'number of stations used: {0}\n'.format(len(self.magnitudes.values())))
|
||||
'number of stations used: {0}\n'.format(len(self.magnitudes.values())))
|
||||
print('\tstation\tmagnitude')
|
||||
for s, m in self.magnitudes.items(): print('\t{0}\t{1}'.format(s, m))
|
||||
|
||||
@@ -126,8 +126,8 @@ class Magnitude(object):
|
||||
# scaling necessary
|
||||
print("Scaling network magnitude ...")
|
||||
mag = ope.Magnitude(
|
||||
mag=np.median([M.mag for M in self.magnitudes.values()]) *\
|
||||
magscaling[0] + magscaling[1],
|
||||
mag=np.median([M.mag for M in self.magnitudes.values()]) * \
|
||||
magscaling[0] + magscaling[1],
|
||||
magnitude_type=self.type,
|
||||
origin_id=self.origin_id,
|
||||
station_count=len(self.magnitudes),
|
||||
@@ -215,7 +215,7 @@ class LocalMagnitude(Magnitude):
|
||||
th = np.arange(0, len(sqH) * dt, dt)
|
||||
# get maximum peak within pick window
|
||||
iwin = getsignalwin(th, t0 - stime, self.calc_win)
|
||||
ii = min([iwin[len(iwin)-1], len(th)])
|
||||
ii = min([iwin[len(iwin) - 1], len(th)])
|
||||
iwin = iwin[0:ii]
|
||||
wapp = np.max(sqH[iwin])
|
||||
if self.verbose:
|
||||
@@ -250,8 +250,8 @@ class LocalMagnitude(Magnitude):
|
||||
if not wf:
|
||||
if self.verbose:
|
||||
print(
|
||||
'WARNING: no waveform data found for station {0}'.format(
|
||||
station))
|
||||
'WARNING: no waveform data found for station {0}'.format(
|
||||
station))
|
||||
continue
|
||||
delta = degrees2kilometers(a.distance)
|
||||
onset = pick.time
|
||||
@@ -270,13 +270,14 @@ class LocalMagnitude(Magnitude):
|
||||
if str(self.wascaling) == '[0.0, 0.0, 0.0]':
|
||||
print("Calculating original Richter magnitude ...")
|
||||
magnitude = ope.StationMagnitude(mag=np.log10(a0) \
|
||||
+ richter_magnitude_scaling(delta))
|
||||
+ richter_magnitude_scaling(delta))
|
||||
else:
|
||||
print("Calculating scaled local magnitude ...")
|
||||
a0 = a0 * 1e03 # mm to nm (see Havskov & Ottemöller, 2010)
|
||||
a0 = a0 * 1e03 # mm to nm (see Havskov & Ottemöller, 2010)
|
||||
magnitude = ope.StationMagnitude(mag=np.log10(a0) \
|
||||
+ self.wascaling[0] * np.log10(delta) + self.wascaling[1]
|
||||
* delta + self.wascaling[2])
|
||||
+ self.wascaling[0] * np.log10(delta) + self.wascaling[1]
|
||||
* delta + self.wascaling[
|
||||
2])
|
||||
magnitude.origin_id = self.origin_id
|
||||
magnitude.waveform_id = pick.waveform_id
|
||||
magnitude.amplitude_id = amplitude.resource_id
|
||||
@@ -397,8 +398,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
|
||||
|
||||
if verbosity:
|
||||
print(
|
||||
"calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station {0} ...".format(
|
||||
tr.stats.station))
|
||||
"calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station {0} ...".format(
|
||||
tr.stats.station))
|
||||
|
||||
# additional common parameters for calculating Mo
|
||||
rP = 2 / np.sqrt(
|
||||
@@ -412,8 +413,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
|
||||
|
||||
if verbosity:
|
||||
print(
|
||||
"calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(
|
||||
Mo, Mw))
|
||||
"calcMoMw: Calculated seismic moment Mo = {0} Nm => Mw = {1:3.1f} ".format(
|
||||
Mo, Mw))
|
||||
|
||||
return Mo, Mw
|
||||
|
||||
@@ -452,7 +453,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
|
||||
:type: integer
|
||||
'''
|
||||
if verbosity:
|
||||
print ("Calculating source spectrum for station %s ...." % wfstream[0].stats.station)
|
||||
print("Calculating source spectrum for station %s ...." % wfstream[0].stats.station)
|
||||
|
||||
# get Q value
|
||||
Q, A = qp
|
||||
@@ -509,9 +510,9 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
|
||||
zc = crossings_nonzero_all(wfzc)
|
||||
if np.size(zc) == 0 or len(zc) <= 3:
|
||||
if verbosity:
|
||||
print ("calcsourcespec: Something is wrong with the waveform, "
|
||||
"no zero crossings derived!\n")
|
||||
print ("No calculation of source spectrum possible!")
|
||||
print("calcsourcespec: Something is wrong with the waveform, "
|
||||
"no zero crossings derived!\n")
|
||||
print("No calculation of source spectrum possible!")
|
||||
plotflag = 0
|
||||
else:
|
||||
plotflag = 1
|
||||
@@ -558,22 +559,22 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
|
||||
[optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin])
|
||||
w0 = optspecfit[0]
|
||||
fc = optspecfit[1]
|
||||
#w01 = optspecfit[0]
|
||||
#fc1 = optspecfit[1]
|
||||
# w01 = optspecfit[0]
|
||||
# fc1 = optspecfit[1]
|
||||
if verbosity:
|
||||
print ("calcsourcespec: Determined w0-value: %e m/Hz, \n"
|
||||
"calcsourcespec: Determined corner frequency: %f Hz" % (w0, fc))
|
||||
print("calcsourcespec: Determined w0-value: %e m/Hz, \n"
|
||||
"calcsourcespec: Determined corner frequency: %f Hz" % (w0, fc))
|
||||
|
||||
# use of conventional fitting
|
||||
# [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot, verbosity)
|
||||
# use of conventional fitting
|
||||
# [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot, verbosity)
|
||||
|
||||
# get w0 and fc as median of both
|
||||
# source spectrum fits
|
||||
#w0 = np.median([w01, w02])
|
||||
#fc = np.median([fc1, fc2])
|
||||
#if verbosity:
|
||||
# print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (
|
||||
# w0, fc))
|
||||
# get w0 and fc as median of both
|
||||
# source spectrum fits
|
||||
# w0 = np.median([w01, w02])
|
||||
# fc = np.median([fc1, fc2])
|
||||
# if verbosity:
|
||||
# print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (
|
||||
# w0, fc))
|
||||
|
||||
if iplot > 1:
|
||||
f1 = plt.figure()
|
||||
@@ -659,9 +660,9 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
|
||||
# left side of initial corner frequency
|
||||
fcstopl = max(f[0], fc0 - max(1, fc0 / 2))
|
||||
il = np.where(f <= fcstopl)
|
||||
il = il[0][np.size(il) - 1]
|
||||
il = il[0][np.size(il) - 1]
|
||||
# right side of initial corner frequency
|
||||
fcstopr = min(fc0 + (fc0 / 2), f[len(f) - 1])
|
||||
fcstopr = min(fc0 + (fc0 / 2), f[len(f) - 1])
|
||||
ir = np.where(f >= fcstopr)
|
||||
# check, if fcstopr is available
|
||||
if np.size(ir) == 0:
|
||||
@@ -672,16 +673,16 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
|
||||
|
||||
# vary corner frequency around initial point
|
||||
print("fitSourceModel: Varying corner frequency "
|
||||
"around initial corner frequency ...")
|
||||
"around initial corner frequency ...")
|
||||
# check difference of il and ir in order to
|
||||
# keep calculation time acceptable
|
||||
idiff = ir - il
|
||||
if idiff > 10000:
|
||||
increment = 100
|
||||
increment = 100
|
||||
elif idiff <= 20:
|
||||
increment = 1
|
||||
increment = 1
|
||||
else:
|
||||
increment = 10
|
||||
increment = 10
|
||||
|
||||
for i in range(il, ir, increment):
|
||||
FC = f[i]
|
||||
@@ -707,10 +708,10 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
|
||||
w0 = max(S)
|
||||
if verbosity:
|
||||
print(
|
||||
"fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
|
||||
"fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0))
|
||||
|
||||
if iplot > 1:
|
||||
plt.figure()#iplot)
|
||||
plt.figure() # iplot)
|
||||
plt.loglog(f, S, 'k')
|
||||
plt.loglog([f[0], fc], [w0, w0], 'g')
|
||||
plt.loglog([fc, fc], [w0 / 100, w0], 'g')
|
||||
@@ -719,7 +720,7 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False):
|
||||
plt.xlabel('Frequency [Hz]')
|
||||
plt.ylabel('Amplitude [m/Hz]')
|
||||
plt.grid()
|
||||
plt.figure()#iplot + 1)
|
||||
plt.figure() # iplot + 1)
|
||||
plt.subplot(311)
|
||||
plt.plot(f[il:ir], STD, '*')
|
||||
plt.title('Common Standard Deviations')
|
||||
|
||||
@@ -12,6 +12,7 @@ from obspy.core import Stream
|
||||
from pylot.core.pick.utils import getsignalwin
|
||||
from scipy.optimize import curve_fit
|
||||
|
||||
|
||||
class Magnitude(object):
|
||||
'''
|
||||
Superclass for calculating Wood-Anderson peak-to-peak
|
||||
@@ -45,7 +46,6 @@ class Magnitude(object):
|
||||
self.calcwapp()
|
||||
self.calcsourcespec()
|
||||
|
||||
|
||||
def getwfstream(self):
|
||||
return self.wfstream
|
||||
|
||||
@@ -85,6 +85,7 @@ class Magnitude(object):
|
||||
def calcsourcespec(self):
|
||||
self.sourcespek = None
|
||||
|
||||
|
||||
class WApp(Magnitude):
|
||||
'''
|
||||
Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
|
||||
@@ -92,8 +93,8 @@ class WApp(Magnitude):
|
||||
'''
|
||||
|
||||
def calcwapp(self):
|
||||
print ("Getting Wood-Anderson peak-to-peak amplitude ...")
|
||||
print ("Simulating Wood-Anderson seismograph ...")
|
||||
print("Getting Wood-Anderson peak-to-peak amplitude ...")
|
||||
print("Simulating Wood-Anderson seismograph ...")
|
||||
|
||||
self.wapp = None
|
||||
stream = self.getwfstream()
|
||||
@@ -118,7 +119,7 @@ class WApp(Magnitude):
|
||||
# get maximum peak within pick window
|
||||
iwin = getsignalwin(th, self.getTo(), self.getpwin())
|
||||
self.wapp = np.max(sqH[iwin])
|
||||
print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp
|
||||
print("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp
|
||||
|
||||
if self.getiplot() > 1:
|
||||
stream.plot()
|
||||
@@ -143,10 +144,10 @@ class DCfc(Magnitude):
|
||||
'''
|
||||
|
||||
def calcsourcespec(self):
|
||||
print ("Calculating source spectrum ....")
|
||||
print("Calculating source spectrum ....")
|
||||
|
||||
self.w0 = None # DC-value
|
||||
self.fc = None # corner frequency
|
||||
self.w0 = None # DC-value
|
||||
self.fc = None # corner frequency
|
||||
|
||||
stream = self.getwfstream()
|
||||
tr = stream[0]
|
||||
@@ -159,14 +160,14 @@ class DCfc(Magnitude):
|
||||
# fft
|
||||
fny = tr.stats.sampling_rate / 2
|
||||
l = len(xdat) / tr.stats.sampling_rate
|
||||
n = tr.stats.sampling_rate * l # number of fft bins after Bath
|
||||
n = tr.stats.sampling_rate * l # number of fft bins after Bath
|
||||
# find next power of 2 of data length
|
||||
m = pow(2, np.ceil(np.log(len(xdat)) / np.log(2)))
|
||||
N = int(np.power(m, 2))
|
||||
y = tr.stats.delta * np.fft.fft(xdat, N)
|
||||
Y = abs(y[: N/2])
|
||||
Y = abs(y[: N / 2])
|
||||
L = (N - 1) / tr.stats.sampling_rate
|
||||
f = np.arange(0, fny, 1/L)
|
||||
f = np.arange(0, fny, 1 / L)
|
||||
|
||||
# remove zero-frequency and frequencies above
|
||||
# corner frequency of seismometer (assumed
|
||||
@@ -185,20 +186,18 @@ class DCfc(Magnitude):
|
||||
[optspecfit, pcov] = curve_fit(synthsourcespec, F, YY.real, [DCin, Fcin])
|
||||
self.w0 = optspecfit[0]
|
||||
self.fc = optspecfit[1]
|
||||
print ("DCfc: Determined DC-value: %e m/Hz, \n" \
|
||||
"Determined corner frequency: %f Hz" % (self.w0, self.fc))
|
||||
|
||||
|
||||
#if self.getiplot() > 1:
|
||||
iplot=2
|
||||
if iplot > 1:
|
||||
print ("DCfc: Determined DC-value: %e m/Hz, \n"
|
||||
"Determined corner frequency: %f Hz" % (self.w0, self.fc))
|
||||
print("DCfc: Determined DC-value: %e m/Hz, \n" \
|
||||
"Determined corner frequency: %f Hz" % (self.w0, self.fc))
|
||||
|
||||
# if self.getiplot() > 1:
|
||||
iplot = 2
|
||||
if iplot > 1:
|
||||
print("DCfc: Determined DC-value: %e m/Hz, \n"
|
||||
"Determined corner frequency: %f Hz" % (self.w0, self.fc))
|
||||
|
||||
if self.getiplot() > 1:
|
||||
f1 = plt.figure()
|
||||
plt.subplot(2,1,1)
|
||||
plt.subplot(2, 1, 1)
|
||||
# show displacement in mm
|
||||
plt.plot(t, np.multiply(tr, 1000), 'k')
|
||||
plt.plot(t[iwin], np.multiply(xdat, 1000), 'g')
|
||||
@@ -206,12 +205,12 @@ class DCfc(Magnitude):
|
||||
plt.xlabel('Time since %s' % tr.stats.starttime)
|
||||
plt.ylabel('Displacement [mm]')
|
||||
|
||||
plt.subplot(2,1,2)
|
||||
plt.subplot(2, 1, 2)
|
||||
plt.loglog(f, Y.real, 'k')
|
||||
plt.loglog(F, YY.real)
|
||||
plt.loglog(F, fit, 'g')
|
||||
plt.title('Source Spectrum from P Pulse, DC=%e m/Hz, fc=%4.1f Hz' \
|
||||
% (self.w0, self.fc))
|
||||
% (self.w0, self.fc))
|
||||
plt.xlabel('Frequency [Hz]')
|
||||
plt.ylabel('Amplitude [m/Hz]')
|
||||
plt.grid()
|
||||
@@ -235,8 +234,7 @@ def synthsourcespec(f, omega0, fcorner):
|
||||
:type: float
|
||||
'''
|
||||
|
||||
#ssp = omega0 / (pow(2, (1 + f / fcorner)))
|
||||
# ssp = omega0 / (pow(2, (1 + f / fcorner)))
|
||||
ssp = omega0 / (1 + pow(2, (f / fcorner)))
|
||||
|
||||
return ssp
|
||||
|
||||
|
||||
Reference in New Issue
Block a user