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

[change] do some major refactoring on Magnitude and subclasses to be more efficient and clean

Browse Source
This commit is contained in:
Sebastian Wehling-Benatelli 2016-09-26 10:49:02 +02:00
parent 51f4082e04
commit eaa3c2e75d
3 changed files with 69 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

10
autoPyLoT.py
View File

@ -10,7 +10,7 @@ import os
import numpy as np import numpy as np
from pylot.core.analysis.magnitude import M0Mw from pylot.core.analysis.magnitude import MomentMagnitude
from pylot.core.io.data import Data from pylot.core.io.data import Data
from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.inputs import AutoPickParameter
import pylot.core.loc.nll as nll import pylot.core.loc.nll as nll
@ -143,7 +143,7 @@ def autoPyLoT(inputfile):
# get latest NLLoc-location file if several are available # get latest NLLoc-location file if several are available
nllocfile = max(glob.glob(locsearch), key=os.path.getctime) nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \ nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \ parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir')) parameter.get('invdir'))
@ -185,7 +185,7 @@ def autoPyLoT(inputfile):
nlloccounter = maxnumit nlloccounter = maxnumit
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \ nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \ parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir')) parameter.get('invdir'))
@ -260,7 +260,7 @@ def autoPyLoT(inputfile):
# get latest NLLOc-location file if several are available # get latest NLLOc-location file if several are available
nllocfile = max(glob.glob(locsearch), key=os.path.getctime) nllocfile = max(glob.glob(locsearch), key=os.path.getctime)
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \ nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \ parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir')) parameter.get('invdir'))
@ -302,7 +302,7 @@ def autoPyLoT(inputfile):
nlloccounter = maxnumit nlloccounter = maxnumit
# calculating seismic moment Mo and moment magnitude Mw # calculating seismic moment Mo and moment magnitude Mw
finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \
nllocfile, picks, parameter.get('rho'), \ nllocfile, picks, parameter.get('rho'), \
parameter.get('vp'), parameter.get('Qp'), \ parameter.get('vp'), parameter.get('Qp'), \
parameter.get('invdir')) parameter.get('invdir'))

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

@ -9,7 +9,7 @@ import os
import matplotlib.pyplot as plt import matplotlib.pyplot as plt
import numpy as np import numpy as np
from obspy.core import Stream, UTCDateTime from obspy.core import Stream, UTCDateTime
from obspy.core.event import Magnitude as Obspymagnitude from obspy.core.event import Magnitude as ObspyMagnitude
from obspy.geodetics import degrees2kilometers from obspy.geodetics import degrees2kilometers
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
@ -33,16 +33,16 @@ class Magnitude(object):
and moment magnitudes. and moment magnitudes.
''' '''
def __init__(self, wfstream, To, pwin, iplot, NLLocfile=None, \ def __init__(self, wfstream, t0, pwin, iplot, NLLocfile=None, \
picks=None, rho=None, vp=None, Qp=None, invdir=None): picks=None, rho=None, vp=None, Qp=None, invdir=None):
''' '''
:param: wfstream :param: wfstream
:type: `~obspy.core.stream.Stream :type: `~obspy.core.stream.Stream
:param: To, onset time, P- or S phase :param: t0, onset time, P- or S phase
:type: float :type: float
:param: pwin, pick window [To To+pwin] to get maximum :param: pwin, pick window [t0 t0+pwin] to get maximum
peak-to-peak amplitude (WApp) or to calculate peak-to-peak amplitude (WApp) or to calculate
source spectrum (DCfc) around P onset source spectrum (DCfc) around P onset
:type: float :type: float
@ -69,37 +69,51 @@ class Magnitude(object):
assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream) assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream)
self.setwfstream(wfstream) self._stream = wfstream
self.setTo(To) self._invdir = invdir
self.setpwin(pwin) self._t0 = t0
self._pwin = pwin
self.setiplot(iplot) self.setiplot(iplot)
self.setNLLocfile(NLLocfile) self.setNLLocfile(NLLocfile)
self.setrho(rho) self.setrho(rho)
self.setpicks(picks) self.setpicks(picks)
self.setvp(vp) self.setvp(vp)
self.setQp(Qp) self.setQp(Qp)
self.setinvdir(invdir)
self.calcwapp() self.calcwapp()
self.calcsourcespec() self.calcsourcespec()
self.run_calcMoMw() self.run_calcMoMw()
def getwfstream(self): @property
return self.wfstream def stream(self):
return self._stream
def setwfstream(self, wfstream): @stream.setter
self.wfstream = wfstream def stream(self, wfstream):
self._stream = wfstream
def getTo(self): @property
return self.To def t0(self):
return self._t0
def setTo(self, To): @t0.setter
self.To = To def t0(self, value):
self._t0 = value
def getpwin(self): @property
return self.pwin def invdir(self):
return self._invdir
def setpwin(self, pwin): @invdir.setter
self.pwin = pwin def invdir(self, value):
self._invdir = value
@property
def pwin(self):
return self._pwin
@pwin.setter
def pwin(self, value):
self._pwin = value
def getiplot(self): def getiplot(self):
return self.iplot return self.iplot
@ -146,12 +160,12 @@ class Magnitude(object):
def getfc(self): def getfc(self):
return self.fc return self.fc
def setinvdir(self, invdir):
self.invdir = invdir
def get_metadata(self): def get_metadata(self):
return read_metadata(self.invdir) return read_metadata(self.invdir)
def plot(self):
pass
def getpicdic(self): def getpicdic(self):
return self.picdic return self.picdic
@ -165,7 +179,7 @@ class Magnitude(object):
self.pickdic = None self.pickdic = None
class WApp(Magnitude): class RichterMagnitude(Magnitude):
''' '''
Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- Method to derive peak-to-peak amplitude as seen on a Wood-Anderson-
seismograph. Has to be derived from instrument corrected traces! seismograph. Has to be derived from instrument corrected traces!
@ -176,7 +190,7 @@ class WApp(Magnitude):
print ("Simulating Wood-Anderson seismograph ...") print ("Simulating Wood-Anderson seismograph ...")
self.wapp = None self.wapp = None
stream = self.getwfstream() stream = self.stream
# poles, zeros and sensitivity of WA seismograph # poles, zeros and sensitivity of WA seismograph
# (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716)
@ -196,7 +210,7 @@ class WApp(Magnitude):
# get time array # get time array
th = np.arange(0, len(sqH) * stream[0].stats.delta, stream[0].stats.delta) th = np.arange(0, len(sqH) * stream[0].stats.delta, stream[0].stats.delta)
# get maximum peak within pick window # get maximum peak within pick window
iwin = getsignalwin(th, self.getTo(), self.getpwin()) iwin = getsignalwin(th, self.t0, self.pwin)
self.wapp = np.max(sqH[iwin]) 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
@ -205,7 +219,7 @@ class WApp(Magnitude):
f = plt.figure(2) f = plt.figure(2)
plt.plot(th, sqH) plt.plot(th, sqH)
plt.plot(th[iwin], sqH[iwin], 'g') plt.plot(th[iwin], sqH[iwin], 'g')
plt.plot([self.getTo(), self.getTo()], [0, max(sqH)], 'r', linewidth=2) plt.plot([self.t0, self.t0], [0, max(sqH)], 'r', linewidth=2)
plt.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \ plt.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \
% (stream[0].stats.station, self.wapp)) % (stream[0].stats.station, self.wapp))
plt.xlabel('Time [s]') plt.xlabel('Time [s]')
@ -215,7 +229,7 @@ class WApp(Magnitude):
plt.close(f) plt.close(f)
class M0Mw(Magnitude): class MomentMagnitude(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 calcsourcespec for calculating plateau w0 Requires results of class calcsourcespec for calculating plateau w0
@ -229,7 +243,7 @@ class M0Mw(Magnitude):
picks = self.getpicks() picks = self.getpicks()
nllocfile = self.getNLLocfile() nllocfile = self.getNLLocfile()
wfdat = self.getwfstream() wfdat = self.stream
self.picdic = None self.picdic = None
for key in picks: for key in picks:
@ -696,7 +710,7 @@ def calc_richter_magnitude(e, wf, metadata, amp_win):
print('stations used:\n') print('stations used:\n')
for s in mags.keys(): print('\t{0}'.format(s)) for s in mags.keys(): print('\t{0}'.format(s))
return Obspymagnitude(mag=mag, magnitude_type='ML') return ObspyMagnitude(mag=mag, magnitude_type='ML')
else: else:
return None return None
@ -725,6 +739,6 @@ def calc_moment_magnitude(e, wf, metadata, vp, Qp, rho):
print('number of stations used: {0}\n'.format(len(mags.values()))) print('number of stations used: {0}\n'.format(len(mags.values())))
print('stations used:\n') print('stations used:\n')
for s in mags.keys(): print('\t{0}'.format(s)) for s in mags.keys(): print('\t{0}'.format(s))
return Obspymagnitude(mag=mag, magnitude_type='Mw') return ObspyMagnitude(mag=mag, magnitude_type='Mw')
else: else:
return None return None

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

@ -20,7 +20,7 @@ from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.dataprocessing import restitute_data, read_metadata
from pylot.core.util.utils import getPatternLine from pylot.core.util.utils import getPatternLine
from pylot.core.io.data import Data from pylot.core.io.data import Data
from pylot.core.analysis.magnitude import WApp from pylot.core.analysis.magnitude import RichterMagnitude
def autopickevent(data, param): def autopickevent(data, param):
@ -571,11 +571,11 @@ def autopickstation(wfstream, pickparam, verbose=False):
# using subclass WApp of superclass Magnitude # using subclass WApp of superclass Magnitude
if restflag: if restflag:
if Sweight < 4: if Sweight < 4:
wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) wapp = RichterMagnitude(cordat, mpickS, mpickP + sstop, iplot)
else: else:
# use larger window for getting peak-to-peak amplitude # use larger window for getting peak-to-peak amplitude
# as the S pick is quite unsure # as the S pick is quite unsure
wapp = WApp(cordat, mpickP, mpickP + sstop + wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop +
(0.5 * (mpickP + sstop)), iplot) (0.5 * (mpickP + sstop)), iplot)
Ao = wapp.getwapp() Ao = wapp.getwapp()
@ -607,7 +607,7 @@ def autopickstation(wfstream, pickparam, verbose=False):
if restflag is True: if restflag is True:
# calculate WA-peak-to-peak amplitude # calculate WA-peak-to-peak amplitude
# using subclass WApp of superclass Magnitude # using subclass WApp of superclass Magnitude
wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP
+ sstop)), + sstop)),
iplot) iplot)
Ao = wapp.getwapp() Ao = wapp.getwapp()