From 80e0ca99d7ad90128ab6d88989e408a8062fe46f Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Wed, 21 Sep 2016 14:12:58 +0200 Subject: [PATCH 01/32] [new] added function to calculate symmetrized pickerror on the fly --- QtPyLoT.py | 3 +++ pylot/core/pick/utils.py | 13 ++++++++++++- 2 files changed, 15 insertions(+), 1 deletion(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index ef494399..66560ab8 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -47,6 +47,7 @@ from pylot.core.io.data import Data from pylot.core.io.inputs import FilterOptions, AutoPickParameter from pylot.core.pick.autopick import autopickevent from pylot.core.pick.compare import Comparison +from pylot.core.pick.utils import symmetrize_error from pylot.core.io.phases import picksdict_from_picks import pylot.core.loc.nll as nll from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP, \ @@ -894,6 +895,8 @@ class MainWindow(QMainWindow): epp = picks['epp'] - stime lpp = picks['lpp'] - stime spe = picks['spe'] + if not spe: + spe = symmetrize_error(mpp - epp, lpp - mpp) if picktype == 'manual': ax.fill_between([epp, lpp], ylims[0], ylims[1], diff --git a/pylot/core/pick/utils.py b/pylot/core/pick/utils.py index 09405fb1..20c2227c 100644 --- a/pylot/core/pick/utils.py +++ b/pylot/core/pick/utils.py @@ -103,7 +103,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=None, stealth_mode=False): # by weighting latest possible pick two times earliest possible pick diffti_tl = LPick - Pick1 diffti_te = Pick1 - EPick - PickError = (diffti_te + 2 * diffti_tl) / 3 + PickError = symmetrize_error(diffti_te, diffti_tl) if iplot > 1: p = plt.figure(iplot) @@ -325,6 +325,17 @@ def crossings_nonzero_all(data): return ((pos[:-1] & npos[1:]) | (npos[:-1] & pos[1:])).nonzero()[0] +def symmetrize_error(dte, dtl): + """ + takes earliest and latest possible pick and returns the symmetrized pick + uncertainty value + :param dte: relative lower uncertainty + :param dtl: relative upper uncertainty + :return: symmetrized error + """ + return (dte + 2 * dtl) / 3 + + def getSNR(X, TSNR, t1, tracenum=0): ''' Function to calculate SNR of certain part of seismogram relative to From 180cd25b512c499af71bf7a45346b062b36e8604 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 20 Sep 2016 14:30:24 +0200 Subject: [PATCH 02/32] Fixed bug in read_metadata.py: path to inventory file was not taken into account. --- pylot/core/util/dataprocessing.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index 497f4762..7ec89e7d 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -183,7 +183,8 @@ def read_metadata(path_to_inventory): "RESP-file found!") elif invtype == 'dless': # prevent multiple read of large dlsv if len(inv[invtype]) == 1: - robj = Parser(inv[invtype][0]) + fullpath_inv = os.path.join(path_to_inventory, inv[invtype][0]) + robj = Parser(fullpath_inv) else: robj = inv[invtype] else: From f35559e7c022d8b4e064fc64aaecce1bb7695494 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 22 Sep 2016 10:53:09 +0200 Subject: [PATCH 03/32] [new] added data file and a function that evaluates the Gutenberg-Richter relation for a given distance --- inputs/gutenberg_richter.data | 53 ++++++++++++++++++++++++++++++++ pylot/core/analysis/magnitude.py | 10 ++++-- 2 files changed, 61 insertions(+), 2 deletions(-) create mode 100644 inputs/gutenberg_richter.data diff --git a/inputs/gutenberg_richter.data b/inputs/gutenberg_richter.data new file mode 100644 index 00000000..74909881 --- /dev/null +++ b/inputs/gutenberg_richter.data @@ -0,0 +1,53 @@ + 0 1.4 + 10 1.5 + 20 1.7 + 25 1.9 + 30 2.1 + 35 2.3 + 40 2.4 + 45 2.5 + 50 2.6 + 60 2.8 + 70 2.8 + 75 2.9 + 85 2.9 + 90 3.0 +100 3.0 +110 3.1 +120 3.1 +130 3.2 +140 3.2 +150 3.3 +160 3.3 +170 3.4 +180 3.4 +190 3.5 +200 3.5 +210 3.6 +230 3.7 +240 3.7 +250 3.8 +260 3.8 +270 3.9 +280 3.9 +290 4.0 +300 4.0 +310 4.1 +320 4.2 +330 4.2 +340 4.2 +350 4.3 +360 4.3 +370 4.3 +380 4.4 +390 4.4 +400 4.5 +430 4.6 +470 4.7 +510 4.8 +560 4.9 +600 5.1 +700 5.2 +800 5.4 +900 5.5 +1000 5.7 \ No newline at end of file diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index b65c9ab2..14af46a0 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -5,7 +5,7 @@ Created autumn/winter 2015. :author: Ludger Küperkoch / MAGS2 EP3 working group """ - +import os import matplotlib.pyplot as plt import numpy as np from obspy.core import Stream, UTCDateTime @@ -15,8 +15,14 @@ from scipy.optimize import curve_fit from scipy import integrate, signal from pylot.core.io.data import Data from pylot.core.util.dataprocessing import restitute_data, read_metadata -from pylot.core.util.utils import common_range +from pylot.core.util.utils import common_range, fit_curve +def gutenberg_richter_relation(delta): + relation = np.loadtxt(os.path.join(os.path.expanduser('~'), + '.pylot', 'gutenberg_richter.data')) + # prepare spline interpolation to calculate return value + func, params = fit_curve(relation[:,0], relation[:, 1]) + return func(delta, params) class Magnitude(object): ''' From 8307974edf8f0665567c6b37a421e9d7eb204ca5 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 22 Sep 2016 11:39:07 +0200 Subject: [PATCH 04/32] [new] added richter magnitude calculation (to be tested) --- QtPyLoT.py | 47 ++++++++++++++++++++++++++++++-- pylot/core/analysis/magnitude.py | 36 +++++++++++++++++++++++- 2 files changed, 79 insertions(+), 4 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 66560ab8..f794ec7a 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -42,12 +42,13 @@ from obspy import UTCDateTime from obspy.geodetics import degrees2kilometers from obspy.core.event import Magnitude -from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw +from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw, \ + calc_woodanderson_pp, gutenberg_richter_relation from pylot.core.io.data import Data from pylot.core.io.inputs import FilterOptions, AutoPickParameter from pylot.core.pick.autopick import autopickevent from pylot.core.pick.compare import Comparison -from pylot.core.pick.utils import symmetrize_error +from pylot.core.pick.utils import symmetrize_error, select_for_phase from pylot.core.io.phases import picksdict_from_picks import pylot.core.loc.nll as nll from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP, \ @@ -972,7 +973,47 @@ class MainWindow(QMainWindow): self.get_data().applyEVTData(lt.read_location(locpath), type='event') self.get_data().get_evt_data().magnitudes.append(self.calc_magnitude()) - def calc_magnitude(self): + + def calc_magnitude(self, type='ML'): + if type == 'ML': + return self.calc_richter_magnitude() + elif type == 'Mw': + return self.calc_moment_magnitude() + else: + return None + + + def calc_richter_magnitude(self): + e = self.get_data().get_evt_data() + if e.origins: + o = e.origins[0] + mags = dict() + for a in o.arrivals: + if a.phase not in 'sS': + continue + pick = a.pick_id.get_referred_object() + station = pick.waveform_id.station_code + wf = select_for_phase(self.get_data().getWFData().select( + station=station), a.phase) + delta = degrees2kilometers(a.distance) + wapp = calc_woodanderson_pp(wf, pick.time, self.inputs.get( + 'sstop')) + # using standard Gutenberg-Richter relation + # TODO make the ML calculation more flexible by allowing + # use of custom relation functions + mag = np.log10(wapp) + gutenberg_richter_relation(delta) + mags[station] = mag + mag = np.median([M[1] for M in mags.values()]) + # give some information on the processing + print('number of stations used: {0}\n'.format(len(mags.values()))) + print('stations used:\n') + for s in mags.keys(): print('\t{0}'.format(s)) + + return Magnitude(mag=mag, magnitude_type='ML') + else: + return None + + def calc_moment_magnitude(self): e = self.get_data().get_evt_data() settings = QSettings() if e.origins: diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 14af46a0..d913c75d 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -271,6 +271,41 @@ class M0Mw(Magnitude): self.picdic = picks +def calc_woodanderson_pp(st, T0, win=10., verbosity=False): + if verbosity: + print ("Getting Wood-Anderson peak-to-peak amplitude ...") + print ("Simulating Wood-Anderson seismograph ...") + + # poles, zeros and sensitivity of WA seismograph + # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) + paz_wa = { + 'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j], + 'zeros': [0j, 0j], + 'gain': 2080, + 'sensitivity': 1 + } + + stime, etime = common_range(st) + st.trim(stime, etime) + + dt = st[0].stats.delta + + st.simulate(paz_remove=None, paz_simulate=paz_wa) + + # get RMS of both horizontal components + sqH = np.sqrt(np.sum([np.power(tr.data, 2) for tr in st if + tr.stats.channel[-1] not in 'Z3'])) + # get time array + th = np.arange(0, len(sqH) * dt, dt) + # get maximum peak within pick window + iwin = getsignalwin(th, T0, win) + wapp = np.max(sqH[iwin]) + if verbosity: + print("Determined Wood-Anderson peak-to-peak amplitude: {0} " + "mm".format(wapp)) + return wapp + + def calcMoMw(wfstream, w0, rho, vp, delta): ''' Subfunction of run_calcMoMw to calculate individual @@ -358,7 +393,6 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, fc = None w0 = None - data = Data() wf_copy = wfstream.copy() invtype, inventory = metadata From 04ec43c699bf3d155c1ffaffae13aa800e83c173 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 22 Sep 2016 14:12:24 +0200 Subject: [PATCH 05/32] [fix] restitute waveform data prior to Wood-Anderson simulation --- QtPyLoT.py | 26 +++++++++++++++++++++----- pylot/core/analysis/magnitude.py | 18 ++++++++++++++---- pylot/core/pick/autopick.py | 2 +- 3 files changed, 36 insertions(+), 10 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index f794ec7a..df879554 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -126,6 +126,11 @@ class MainWindow(QMainWindow): # load and display waveform data self.dirty = False self.load_data() + finv = settings.value("inventoryFile", None) + if finv is not None: + self._metadata = read_metadata(finv) + else: + self._metadata = None if self.loadWaveformData(): self.updateFilterOptions() else: @@ -370,6 +375,17 @@ class MainWindow(QMainWindow): self.setCentralWidget(_widget) + + @property + def metadata(self): + return self._metadata + + + @metadata.setter + def metadata(self, value): + self._metadata = value + + def updateFileMenu(self): self.fileMenu.clear() @@ -996,14 +1012,14 @@ class MainWindow(QMainWindow): wf = select_for_phase(self.get_data().getWFData().select( station=station), a.phase) delta = degrees2kilometers(a.distance) - wapp = calc_woodanderson_pp(wf, pick.time, self.inputs.get( - 'sstop')) + wapp = calc_woodanderson_pp(wf, self.metadata, pick.time, + self.inputs.get('sstop')) # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # use of custom relation functions mag = np.log10(wapp) + gutenberg_richter_relation(delta) mags[station] = mag - mag = np.median([M[1] for M in mags.values()]) + mag = np.median([M for M in mags.values()]) # give some information on the processing print('number of stations used: {0}\n'.format(len(mags.values()))) print('stations used:\n') @@ -1020,7 +1036,7 @@ class MainWindow(QMainWindow): o = e.origins[0] mags = dict() fninv = settings.value("inventoryFile", None) - if fninv is None: + if fninv is None and not self.metadata: fninv, _ = QFileDialog.getOpenFileName(self, self.tr( "Select inventory..."), self.tr("Select file")) ans = QMessageBox.question(self, self.tr("Make default..."), @@ -1032,7 +1048,7 @@ class MainWindow(QMainWindow): if ans == QMessageBox.Yes: settings.setValue("inventoryFile", fninv) settings.sync() - metadata = read_metadata(fninv) + self.metadata = read_metadata(fninv) for a in o.arrivals: if a.phase in 'sS': continue diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index d913c75d..f987a4ae 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -271,7 +271,7 @@ class M0Mw(Magnitude): self.picdic = picks -def calc_woodanderson_pp(st, T0, win=10., verbosity=False): +def calc_woodanderson_pp(st, metadata, T0, win=10., verbosity=False): if verbosity: print ("Getting Wood-Anderson peak-to-peak amplitude ...") print ("Simulating Wood-Anderson seismograph ...") @@ -288,17 +288,27 @@ def calc_woodanderson_pp(st, T0, win=10., verbosity=False): stime, etime = common_range(st) st.trim(stime, etime) + invtype, inventory = metadata + + st = restitute_data(st, invtype, inventory) + dt = st[0].stats.delta st.simulate(paz_remove=None, paz_simulate=paz_wa) # get RMS of both horizontal components - sqH = np.sqrt(np.sum([np.power(tr.data, 2) for tr in st if - tr.stats.channel[-1] not in 'Z3'])) + power = [np.power(tr.data, 2) for tr in st if tr.stats.channel[-1] not + in 'Z3'] + if len(power) != 2: + raise ValueError('Wood-Anderson amplitude defintion only valid for ' + 'two horizontals: {0} given'.format(len(power))) + power_sum = power[0] + power[1] + sqH = np.sqrt(power_sum) + # get time array th = np.arange(0, len(sqH) * dt, dt) # get maximum peak within pick window - iwin = getsignalwin(th, T0, win) + iwin = getsignalwin(th, T0 - stime, win) wapp = np.max(sqH[iwin]) if verbosity: print("Determined Wood-Anderson peak-to-peak amplitude: {0} " diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 3b42eb1a..553a7196 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -604,7 +604,7 @@ def autopickstation(wfstream, pickparam, verbose=False): hdat += ndat h_copy = hdat.copy() [cordat, restflag] = restitute_data(h_copy, invtype, inventory) - if restflag == 1: + if restflag is True: # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP From bfa7ffc960de7484ce642d6d2c878be47fe3358c Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Fri, 23 Sep 2016 15:12:04 +0200 Subject: [PATCH 06/32] [move] moving functions for Richter and moment magnitude calculation to magnitude module for re-use in autoPyLoT --- QtPyLoT.py | 112 ++++++------------------------- pylot/core/analysis/magnitude.py | 67 +++++++++++++++++- 2 files changed, 86 insertions(+), 93 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index df879554..7a64fca2 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -23,10 +23,11 @@ https://www.iconfinder.com/iconsets/flavour (http://www.gnu.org/copyleft/lesser.html) """ -import matplotlib import os import sys +import matplotlib + matplotlib.use('Qt4Agg') matplotlib.rcParams['backend.qt4'] = 'PySide' @@ -37,18 +38,14 @@ from PySide.QtGui import QMainWindow, QInputDialog, QIcon, QFileDialog, \ QDialog, QErrorMessage, QApplication, QPixmap, QMessageBox, QSplashScreen, \ QActionGroup, QListWidget, QDockWidget, QLineEdit import numpy as np -import subprocess from obspy import UTCDateTime -from obspy.geodetics import degrees2kilometers -from obspy.core.event import Magnitude -from pylot.core.analysis.magnitude import calcsourcespec, calcMoMw, \ - calc_woodanderson_pp, gutenberg_richter_relation +from pylot.core.analysis.magnitude import calc_richter_magnitude, calc_moment_magnitude from pylot.core.io.data import Data from pylot.core.io.inputs import FilterOptions, AutoPickParameter from pylot.core.pick.autopick import autopickevent from pylot.core.pick.compare import Comparison -from pylot.core.pick.utils import symmetrize_error, select_for_phase +from pylot.core.pick.utils import symmetrize_error from pylot.core.io.phases import picksdict_from_picks import pylot.core.loc.nll as nll from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP, \ @@ -66,7 +63,6 @@ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.thread import AutoPickThread from pylot.core.util.version import get_git_version as _getVersionString -import icons_rc locateTool = dict(nll=nll) @@ -991,89 +987,25 @@ class MainWindow(QMainWindow): def calc_magnitude(self, type='ML'): - if type == 'ML': - return self.calc_richter_magnitude() - elif type == 'Mw': - return self.calc_moment_magnitude() - else: - return None - - - def calc_richter_magnitude(self): - e = self.get_data().get_evt_data() - if e.origins: - o = e.origins[0] - mags = dict() - for a in o.arrivals: - if a.phase not in 'sS': - continue - pick = a.pick_id.get_referred_object() - station = pick.waveform_id.station_code - wf = select_for_phase(self.get_data().getWFData().select( - station=station), a.phase) - delta = degrees2kilometers(a.distance) - wapp = calc_woodanderson_pp(wf, self.metadata, pick.time, - self.inputs.get('sstop')) - # using standard Gutenberg-Richter relation - # TODO make the ML calculation more flexible by allowing - # use of custom relation functions - mag = np.log10(wapp) + gutenberg_richter_relation(delta) - mags[station] = mag - mag = np.median([M for M in mags.values()]) - # give some information on the processing - print('number of stations used: {0}\n'.format(len(mags.values()))) - print('stations used:\n') - for s in mags.keys(): print('\t{0}'.format(s)) - - return Magnitude(mag=mag, magnitude_type='ML') - else: - return None - - def calc_moment_magnitude(self): - e = self.get_data().get_evt_data() settings = QSettings() - if e.origins: - o = e.origins[0] - mags = dict() - fninv = settings.value("inventoryFile", None) - if fninv is None and not self.metadata: - fninv, _ = QFileDialog.getOpenFileName(self, self.tr( - "Select inventory..."), self.tr("Select file")) - ans = QMessageBox.question(self, self.tr("Make default..."), - self.tr( - "New inventory filename set.\n" + \ - "Do you want to make it the default value?"), - QMessageBox.Yes | QMessageBox.No, - QMessageBox.No) - if ans == QMessageBox.Yes: - settings.setValue("inventoryFile", fninv) - settings.sync() - self.metadata = read_metadata(fninv) - for a in o.arrivals: - if a.phase in 'sS': - continue - pick = a.pick_id.get_referred_object() - station = pick.waveform_id.station_code - wf = self.get_data().getWFData().select(station=station) - if not wf: - continue - onset = pick.time - dist = degrees2kilometers(a.distance) - w0, fc = calcsourcespec(wf, onset, metadata, self.inputs.get('vp'), dist, - a.azimuth, a.takeoff_angle, - self.inputs.get('Qp'), 0) - if w0 is None or fc is None: - continue - station_mag = calcMoMw(wf, w0, self.inputs.get('rho'), - self.inputs.get('vp'), dist) - mags[station] = station_mag - mag = np.median([M[1] for M in mags.values()]) - # give some information on the processing - print('number of stations used: {0}\n'.format(len(mags.values()))) - print('stations used:\n') - for s in mags.keys(): print('\t{0}'.format(s)) - - return Magnitude(mag=mag, magnitude_type='Mw') + fninv = settings.value("inventoryFile", None) + if fninv is None and not self.metadata: + fninv, _ = QFileDialog.getOpenFileName(self, self.tr( + "Select inventory..."), self.tr("Select file")) + ans = QMessageBox.question(self, self.tr("Make default..."), + self.tr( + "New inventory filename set.\n" + \ + "Do you want to make it the default value?"), + QMessageBox.Yes | QMessageBox.No, + QMessageBox.No) + if ans == QMessageBox.Yes: + settings.setValue("inventoryFile", fninv) + settings.sync() + self.metadata = read_metadata(fninv) + if type == 'ML': + return calc_richter_magnitude(self.get_data().get_evt_data(), self.get_data().getWFData(), self.metadata, self.inputs.get('sstop')) + elif type == 'Mw': + return calc_moment_magnitude(self.get_data().get_evt_data(), self.get_data().getWFData(), self.metadata, self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho')) else: return None diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index f987a4ae..73dd422c 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -9,11 +9,13 @@ import os import matplotlib.pyplot as plt import numpy as np from obspy.core import Stream, UTCDateTime -from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all +from obspy.core.event import Magnitude +from obspy.geodetics import degrees2kilometers + +from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, select_for_phase from pylot.core.util.utils import getPatternLine from scipy.optimize import curve_fit from scipy import integrate, signal -from pylot.core.io.data import Data from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.utils import common_range, fit_curve @@ -441,7 +443,7 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, ldat[0].stats.delta)) # get window after P pulse for - # calculating source spectrum + # calculating source spectrum tstart = UTCDateTime(zdat[0].stats.starttime) tonset = onset.timestamp - tstart.timestamp impickP = tonset * zdat[0].stats.sampling_rate @@ -667,3 +669,62 @@ def fitSourceModel(f, S, fc0, iplot): plt.close() return w0, fc + + +def calc_richter_magnitude(e, wf, metadata, amp_win): + if e.origins: + o = e.origins[0] + mags = dict() + for a in o.arrivals: + if a.phase not in 'sS': + continue + pick = a.pick_id.get_referred_object() + station = pick.waveform_id.station_code + wf = select_for_phase(wf.select( + station=station), a.phase) + delta = degrees2kilometers(a.distance) + wapp = calc_woodanderson_pp(wf, metadata, pick.time, + amp_win) + # using standard Gutenberg-Richter relation + # TODO make the ML calculation more flexible by allowing + # use of custom relation functions + mag = np.log10(wapp) + gutenberg_richter_relation(delta) + mags[station] = mag + mag = np.median([M for M in mags.values()]) + # give some information on the processing + print('number of stations used: {0}\n'.format(len(mags.values()))) + print('stations used:\n') + for s in mags.keys(): print('\t{0}'.format(s)) + + return Magnitude(mag=mag, magnitude_type='ML') + else: + return None + + +def calc_moment_magnitude(e, wf, metadata, vp, Qp, rho): + if e.origins: + o = e.origins[0] + mags = dict() + for a in o.arrivals: + if a.phase in 'sS': + continue + pick = a.pick_id.get_referred_object() + station = pick.waveform_id.station_code + wf = wf.select(station=station) + if not wf: + continue + onset = pick.time + dist = degrees2kilometers(a.distance) + w0, fc = calcsourcespec(wf, onset, metadata, vp, dist, a.azimuth, a.takeoff_angle, Qp, 0) + if w0 is None or fc is None: + continue + station_mag = calcMoMw(wf, w0, rho, vp, dist) + mags[station] = station_mag + mag = np.median([M[1] for M in mags.values()]) + # give some information on the processing + print('number of stations used: {0}\n'.format(len(mags.values()))) + print('stations used:\n') + for s in mags.keys(): print('\t{0}'.format(s)) + return Magnitude(mag=mag, magnitude_type='Mw') + else: + return None From 0dffe37d3b86217c1685c3d6393928d598412681 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Fri, 23 Sep 2016 14:28:02 +0200 Subject: [PATCH 07/32] [namefix] rename data file and corresponding function --- inputs/{gutenberg_richter.data => richter_scaling.data} | 0 pylot/core/analysis/magnitude.py | 4 ++-- 2 files changed, 2 insertions(+), 2 deletions(-) rename inputs/{gutenberg_richter.data => richter_scaling.data} (100%) diff --git a/inputs/gutenberg_richter.data b/inputs/richter_scaling.data similarity index 100% rename from inputs/gutenberg_richter.data rename to inputs/richter_scaling.data diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 73dd422c..8c3612a0 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -19,9 +19,9 @@ from scipy import integrate, signal from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.utils import common_range, fit_curve -def gutenberg_richter_relation(delta): +def richter_magnitude_scaling(delta): relation = np.loadtxt(os.path.join(os.path.expanduser('~'), - '.pylot', 'gutenberg_richter.data')) + '.pylot', 'richter_scaling.data')) # prepare spline interpolation to calculate return value func, params = fit_curve(relation[:,0], relation[:, 1]) return func(delta, params) From 51f4082e048d71bf893ef84d39869c9ffa33dc5a Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Fri, 23 Sep 2016 15:21:34 +0200 Subject: [PATCH 08/32] [fix] imported Magnitude overwrite prevented by renamed import; changed wrong function call --- pylot/core/analysis/magnitude.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 8c3612a0..cb1ace0f 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -9,7 +9,7 @@ import os import matplotlib.pyplot as plt import numpy as np from obspy.core import Stream, UTCDateTime -from obspy.core.event import Magnitude +from obspy.core.event import Magnitude as Obspymagnitude from obspy.geodetics import degrees2kilometers from pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, select_for_phase @@ -688,7 +688,7 @@ def calc_richter_magnitude(e, wf, metadata, amp_win): # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # use of custom relation functions - mag = np.log10(wapp) + gutenberg_richter_relation(delta) + mag = np.log10(wapp) + richter_magnitude_scaling(delta) mags[station] = mag mag = np.median([M for M in mags.values()]) # give some information on the processing @@ -696,7 +696,7 @@ def calc_richter_magnitude(e, wf, metadata, amp_win): print('stations used:\n') for s in mags.keys(): print('\t{0}'.format(s)) - return Magnitude(mag=mag, magnitude_type='ML') + return Obspymagnitude(mag=mag, magnitude_type='ML') else: return None @@ -725,6 +725,6 @@ def calc_moment_magnitude(e, wf, metadata, vp, Qp, rho): print('number of stations used: {0}\n'.format(len(mags.values()))) print('stations used:\n') for s in mags.keys(): print('\t{0}'.format(s)) - return Magnitude(mag=mag, magnitude_type='Mw') + return Obspymagnitude(mag=mag, magnitude_type='Mw') else: return None From eaa3c2e75dc2953e360e16ee7d03fbf2d217d9f5 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 26 Sep 2016 10:49:02 +0200 Subject: [PATCH 09/32] [change] do some major refactoring on Magnitude and subclasses to be more efficient and clean --- autoPyLoT.py | 34 +++++++------- pylot/core/analysis/magnitude.py | 76 +++++++++++++++++++------------- pylot/core/pick/autopick.py | 14 +++--- 3 files changed, 69 insertions(+), 55 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 01b430ab..21d7eb52 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -10,7 +10,7 @@ import os 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.inputs import AutoPickParameter import pylot.core.loc.nll as nll @@ -143,10 +143,10 @@ def autoPyLoT(inputfile): # get latest NLLoc-location file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # calculating seismic moment Mo and moment magnitude Mw - finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ + nllocfile, picks, parameter.get('rho'), \ + parameter.get('vp'), parameter.get('Qp'), \ + parameter.get('invdir')) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -185,10 +185,10 @@ def autoPyLoT(inputfile): nlloccounter = maxnumit # calculating seismic moment Mo and moment magnitude Mw - finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ + nllocfile, picks, parameter.get('rho'), \ + parameter.get('vp'), parameter.get('Qp'), \ + parameter.get('invdir')) # get network moment magntiude netMw = [] for key in finalpicks.getpicdic(): @@ -260,10 +260,10 @@ def autoPyLoT(inputfile): # get latest NLLOc-location file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) # calculating seismic moment Mo and moment magnitude Mw - finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ + nllocfile, picks, parameter.get('rho'), \ + parameter.get('vp'), parameter.get('Qp'), \ + parameter.get('invdir')) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -302,10 +302,10 @@ def autoPyLoT(inputfile): nlloccounter = maxnumit # calculating seismic moment Mo and moment magnitude Mw - finalpicks = M0Mw(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ + nllocfile, picks, parameter.get('rho'), \ + parameter.get('vp'), parameter.get('Qp'), \ + parameter.get('invdir')) # get network moment magntiude netMw = [] for key in finalpicks.getpicdic(): diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index cb1ace0f..c6084a08 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -9,7 +9,7 @@ import os import matplotlib.pyplot as plt import numpy as np 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 pylot.core.pick.utils import getsignalwin, crossings_nonzero_all, select_for_phase @@ -33,16 +33,16 @@ class Magnitude(object): 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): ''' :param: wfstream :type: `~obspy.core.stream.Stream - :param: To, onset time, P- or S phase + :param: t0, onset time, P- or S phase :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 source spectrum (DCfc) around P onset :type: float @@ -69,37 +69,51 @@ class Magnitude(object): assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream) - self.setwfstream(wfstream) - self.setTo(To) - self.setpwin(pwin) + self._stream = wfstream + self._invdir = invdir + self._t0 = t0 + self._pwin = pwin self.setiplot(iplot) self.setNLLocfile(NLLocfile) self.setrho(rho) self.setpicks(picks) self.setvp(vp) self.setQp(Qp) - self.setinvdir(invdir) self.calcwapp() self.calcsourcespec() self.run_calcMoMw() - def getwfstream(self): - return self.wfstream + @property + def stream(self): + return self._stream - def setwfstream(self, wfstream): - self.wfstream = wfstream + @stream.setter + def stream(self, wfstream): + self._stream = wfstream - def getTo(self): - return self.To + @property + def t0(self): + return self._t0 - def setTo(self, To): - self.To = To + @t0.setter + def t0(self, value): + self._t0 = value - def getpwin(self): - return self.pwin + @property + def invdir(self): + return self._invdir - def setpwin(self, pwin): - self.pwin = pwin + @invdir.setter + 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): return self.iplot @@ -146,12 +160,12 @@ class Magnitude(object): def getfc(self): return self.fc - def setinvdir(self, invdir): - self.invdir = invdir - def get_metadata(self): return read_metadata(self.invdir) + def plot(self): + pass + def getpicdic(self): return self.picdic @@ -165,7 +179,7 @@ class Magnitude(object): self.pickdic = None -class WApp(Magnitude): +class RichterMagnitude(Magnitude): ''' Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- seismograph. Has to be derived from instrument corrected traces! @@ -176,7 +190,7 @@ class WApp(Magnitude): print ("Simulating Wood-Anderson seismograph ...") self.wapp = None - stream = self.getwfstream() + stream = self.stream # poles, zeros and sensitivity of WA seismograph # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) @@ -196,7 +210,7 @@ class WApp(Magnitude): # get time array th = np.arange(0, len(sqH) * stream[0].stats.delta, stream[0].stats.delta) # 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]) print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp @@ -205,7 +219,7 @@ class WApp(Magnitude): f = plt.figure(2) plt.plot(th, sqH) 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' \ % (stream[0].stats.station, self.wapp)) plt.xlabel('Time [s]') @@ -215,7 +229,7 @@ class WApp(Magnitude): plt.close(f) -class M0Mw(Magnitude): +class MomentMagnitude(Magnitude): ''' Method to calculate seismic moment Mo and moment magnitude Mw. Requires results of class calcsourcespec for calculating plateau w0 @@ -229,7 +243,7 @@ class M0Mw(Magnitude): picks = self.getpicks() nllocfile = self.getNLLocfile() - wfdat = self.getwfstream() + wfdat = self.stream self.picdic = None for key in picks: @@ -696,7 +710,7 @@ def calc_richter_magnitude(e, wf, metadata, amp_win): print('stations used:\n') 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: 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('stations used:\n') 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: return None diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 553a7196..5a017adc 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -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.utils import getPatternLine 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): @@ -571,12 +571,12 @@ def autopickstation(wfstream, pickparam, verbose=False): # using subclass WApp of superclass Magnitude if restflag: if Sweight < 4: - wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) + wapp = RichterMagnitude(cordat, mpickS, mpickP + sstop, iplot) else: # use larger window for getting peak-to-peak amplitude # as the S pick is quite unsure - wapp = WApp(cordat, mpickP, mpickP + sstop + - (0.5 * (mpickP + sstop)), iplot) + wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop + + (0.5 * (mpickP + sstop)), iplot) Ao = wapp.getwapp() else: @@ -607,9 +607,9 @@ def autopickstation(wfstream, pickparam, verbose=False): if restflag is True: # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP - + sstop)), - iplot) + wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP + + sstop)), + iplot) Ao = wapp.getwapp() else: From dc38bd6e79fc53019e1cd660dfa54a5c2989544b Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 26 Sep 2016 14:47:50 +0200 Subject: [PATCH 10/32] [fix, refactor] started major refactoring of magnitude.py and fixed some smaller bugs --- QtPyLoT.py | 1 + pylot/core/analysis/magnitude.py | 115 +++++++++++++++++++++++++++---- pylot/core/util/widgets.py | 2 +- 3 files changed, 103 insertions(+), 15 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 7a64fca2..84d32732 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -63,6 +63,7 @@ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.thread import AutoPickThread from pylot.core.util.version import get_git_version as _getVersionString +import icons_rc locateTool = dict(nll=nll) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index c6084a08..7d6a29e4 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -26,6 +26,43 @@ def richter_magnitude_scaling(delta): func, params = fit_curve(relation[:,0], relation[:, 1]) return func(delta, params) + +class Magnitude(object): + """ + Base class object for Magnitude calculation within PyLoT. + """ + + def __init__(self, stream): + self._stream = stream + self._magnitudes = dict() + + @property + def stream(self): + return self._stream + + @stream.setter + def stream(self, value): + self._stream = value + + @property + def magnitudes(self): + return self._magnitudes + + @magnitudes.setter + def magnitudes(self, value): + """ + takes a tuple and saves the key value pair to private + attribute _magnitudes + :param value: station, magnitude value pair + :type value: tuple or list + :return: + """ + station, magnitude = value + self._magnitudes[station] = magnitude + + def get(self): + return self.magnitudes + class Magnitude(object): ''' Superclass for calculating Wood-Anderson peak-to-peak @@ -73,7 +110,7 @@ class Magnitude(object): self._invdir = invdir self._t0 = t0 self._pwin = pwin - self.setiplot(iplot) + self._iplot = iplot self.setNLLocfile(NLLocfile) self.setrho(rho) self.setpicks(picks) @@ -115,11 +152,13 @@ class Magnitude(object): def pwin(self, value): self._pwin = value - def getiplot(self): + @property + def plot_flag(self): return self.iplot - def setiplot(self, iplot): - self.iplot = iplot + @plot_flag.setter + def plot_flag(self, value): + self._iplot = value def setNLLocfile(self, NLLocfile): self.NLLocfile = NLLocfile @@ -185,6 +224,62 @@ class RichterMagnitude(Magnitude): seismograph. Has to be derived from instrument corrected traces! ''' + # poles, zeros and sensitivity of WA seismograph + # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) + _paz = { + 'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j], + 'zeros': [0j, 0j], + 'gain': 2080, + 'sensitivity': 1 + } + + def __init__(self, stream, t0, calc_win): + super(RichterMagnitude, self).__init__(stream) + + self._t0 = t0 + self._calc_win = calc_win + + @property + def t0(self): + return self._t0 + + @t0.setter + def t0(self, value): + self._t0 = value + + @property + def calc_win(self): + return self._calc_win + + @calc_win.setter + def calc_win(self, value): + self._calc_win = value + + def get(self): + + st = self.stream + + # simulate Wood-Anderson response + st.simulate(paz_remove=None, paz_simulate=self._paz) + + # trim waveform to common range + stime, etime = common_range(st) + st.trim(stime, etime) + + # get time delta from waveform data + dt = st[0].stats.delta + + power = [np.power(tr.data, 2) for tr in st if tr.stats.channel[-1] not + in 'Z3'] + if len(power) != 2: + raise ValueError('Wood-Anderson amplitude defintion only valid for ' + 'two horizontals: {0} given'.format(len(power))) + power_sum = power[0] + power[1] + # + sqH = np.sqrt(power_sum) + + + def calcwapp(self): print ("Getting Wood-Anderson peak-to-peak amplitude ...") print ("Simulating Wood-Anderson seismograph ...") @@ -192,14 +287,6 @@ class RichterMagnitude(Magnitude): self.wapp = None stream = self.stream - # poles, zeros and sensitivity of WA seismograph - # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) - paz_wa = { - 'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j], - 'zeros': [0j, 0j], - 'gain': 2080, - 'sensitivity': 1} - stream.simulate(paz_remove=None, paz_simulate=paz_wa) trH1 = stream[0].data @@ -214,7 +301,7 @@ class RichterMagnitude(Magnitude): self.wapp = np.max(sqH[iwin]) print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp - if self.getiplot() > 1: + if self.plot_flag > 1: stream.plot() f = plt.figure(2) plt.plot(th, sqH) @@ -266,7 +353,7 @@ class MomentMagnitude(Magnitude): # and to derive w0 and fc [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \ self.get_metadata(), self.getvp(), delta, az, \ - inc, self.getQp(), self.getiplot()) + inc, self.getQp(), self.plot_flag) if w0 is not None: # call subfunction to calculate Mo and Mw diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 4d3e3d16..c81a94f3 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -34,7 +34,7 @@ from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS, LOCTOOLS, \ COMPPOSITION_MAP from pylot.core.util.utils import prepTimeAxis, full_range, scaleWFData, \ demeanTrace, isSorted, findComboBoxIndex, clims - +import icons_rc def getDataType(parent): type = QInputDialog().getItem(parent, "Select phases type", "Type:", From 9288a169a4aa2e16eb8d7fec069b9b5189b55e5e Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 26 Sep 2016 14:45:36 +0200 Subject: [PATCH 11/32] [change] if folder selection is canceled do not empty editable text --- pylot/core/util/widgets.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index c81a94f3..da5740ba 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -1381,7 +1381,9 @@ class LocalisationTab(PropTab): def selectDirectory(self, edit): selected_directory = QFileDialog.getExistingDirectory() - edit.setText(selected_directory) + # check if string is empty + if selected_directory: + edit.setText(selected_directory) def getValues(self): loctool = self.locToolComboBox.currentText() From c52277e4a2271ea653b25daee84167c0b6235611 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 26 Sep 2016 15:55:56 +0200 Subject: [PATCH 12/32] [new] added attributes, properties and special method __str__ to the Magnitude superclass -> improves significantly convenience of sub-class programming --- pylot/core/analysis/magnitude.py | 30 +++++++++++++++++++++++++++--- 1 file changed, 27 insertions(+), 3 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 7d6a29e4..4a3cef8f 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -32,10 +32,26 @@ class Magnitude(object): Base class object for Magnitude calculation within PyLoT. """ - def __init__(self, stream): + def __init__(self, stream, event, verbosity=False): + self._plot_flag = False + self._verbosity = verbosity + self._event = event self._stream = stream self._magnitudes = dict() + def __str__(self): + print('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)) + + @property + def plot_flag(self): + return self._plot_flag + + @plot_flag.setter + def plot_flag(self, value): + self._plot_flag = value + @property def stream(self): return self._stream @@ -44,6 +60,14 @@ class Magnitude(object): def stream(self, value): self._stream = value + @property + def event(self): + return self._event + + @property + def arrivals(self): + return self._event.origins[0].arrivals + @property def magnitudes(self): return self._magnitudes @@ -233,8 +257,8 @@ class RichterMagnitude(Magnitude): 'sensitivity': 1 } - def __init__(self, stream, t0, calc_win): - super(RichterMagnitude, self).__init__(stream) + def __init__(self, stream, event, t0, calc_win, verbosity=False): + super(RichterMagnitude, self).__init__(stream, event, verbosity) self._t0 = t0 self._calc_win = calc_win From d4481e4acdf5d7d5e08ee116a14fbeaf4ab15b37 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Mon, 26 Sep 2016 16:04:09 +0200 Subject: [PATCH 13/32] [new] added peak_to_peak, get and net_magnitude giving Wood-Anderson simulated peak amplitude, single station magnitudes and network magnitude for a given event, respectively --- pylot/core/analysis/magnitude.py | 57 ++++++++++++++++++-------------- 1 file changed, 32 insertions(+), 25 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 4a3cef8f..0d975ef1 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -279,9 +279,7 @@ class RichterMagnitude(Magnitude): def calc_win(self, value): self._calc_win = value - def get(self): - - st = self.stream + def peak_to_peak(self, st): # simulate Wood-Anderson response st.simulate(paz_remove=None, paz_simulate=self._paz) @@ -302,43 +300,52 @@ class RichterMagnitude(Magnitude): # sqH = np.sqrt(power_sum) - - - def calcwapp(self): - print ("Getting Wood-Anderson peak-to-peak amplitude ...") - print ("Simulating Wood-Anderson seismograph ...") - - self.wapp = None - stream = self.stream - - stream.simulate(paz_remove=None, paz_simulate=paz_wa) - - trH1 = stream[0].data - trH2 = stream[1].data - ilen = min([len(trH1), len(trH2)]) - # get RMS of both horizontal components - sqH = np.sqrt(np.power(trH1[0:ilen], 2) + np.power(trH2[0:ilen], 2)) # get time array - th = np.arange(0, len(sqH) * stream[0].stats.delta, stream[0].stats.delta) + th = np.arange(0, len(sqH) * dt, dt) # get maximum peak within pick window - iwin = getsignalwin(th, self.t0, self.pwin) - self.wapp = np.max(sqH[iwin]) - print ("Determined Wood-Anderson peak-to-peak amplitude: %f mm") % self.wapp + iwin = getsignalwin(th, self.t0 - stime, self.calc_win) + wapp = np.max(sqH[iwin]) + if self._verbosity: + print("Determined Wood-Anderson peak-to-peak amplitude: {0} " + "mm".format(wapp)) + # check for plot flag (for debugging only) if self.plot_flag > 1: - stream.plot() + st.plot() f = plt.figure(2) plt.plot(th, sqH) plt.plot(th[iwin], sqH[iwin], 'g') 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' \ - % (stream[0].stats.station, self.wapp)) + % (st[0].stats.station, wapp)) plt.xlabel('Time [s]') plt.ylabel('Displacement [mm]') plt.show() raw_input() plt.close(f) + return wapp + + def get(self): + for a in self.arrivals: + if a.phase not in 'sS': + continue + pick = a.pick_id.get_referred_object() + station = pick.waveform_id.station_code + wf = select_for_phase(self.stream.select( + station=station), a.phase) + delta = degrees2kilometers(a.distance) + wapp = self.peak_to_peak(wf) + # using standard Gutenberg-Richter relation + # TODO make the ML calculation more flexible by allowing + # use of custom relation functions + mag = np.log10(wapp) + richter_magnitude_scaling(delta) + self.magnitudes = (station, mag) + return self.magnitudes + + def net_magnitude(self): + return np.median([M for M in self.magnitudes.values()]) + class MomentMagnitude(Magnitude): ''' From 405402ffdc414fde7193a00edded6e17fd16900e Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Tue, 27 Sep 2016 13:57:14 +0200 Subject: [PATCH 14/32] [refactor] major refactoring of Magnitude objects finished now the changed usage of the Magnitude object has to be implemented into autoPyLoT and QtPyLoT (pending) --- pylot/core/analysis/magnitude.py | 542 +++++++++++++------------------ 1 file changed, 223 insertions(+), 319 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 0d975ef1..48417557 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -39,39 +39,52 @@ class Magnitude(object): self._stream = stream self._magnitudes = dict() + def __str__(self): print('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)) + + def __nonzero__(self): + return bool(self.magnitudes) + + @property def plot_flag(self): return self._plot_flag + @plot_flag.setter def plot_flag(self, value): self._plot_flag = value + @property def stream(self): return self._stream + @stream.setter def stream(self, value): self._stream = value + @property def event(self): return self._event + @property def arrivals(self): return self._event.origins[0].arrivals + @property def magnitudes(self): return self._magnitudes + @magnitudes.setter def magnitudes(self, value): """ @@ -84,169 +97,22 @@ class Magnitude(object): station, magnitude = value self._magnitudes[station] = magnitude + def get(self): return self.magnitudes -class Magnitude(object): - ''' - Superclass for calculating Wood-Anderson peak-to-peak - amplitudes, local magnitudes, source spectra, seismic moments - and moment magnitudes. - ''' - def __init__(self, wfstream, t0, pwin, iplot, NLLocfile=None, \ - picks=None, rho=None, vp=None, Qp=None, invdir=None): - ''' - :param: wfstream - :type: `~obspy.core.stream.Stream - - :param: t0, onset time, P- or S phase - :type: float - - :param: pwin, pick window [t0 t0+pwin] to get maximum - peak-to-peak amplitude (WApp) or to calculate - source spectrum (DCfc) around P onset - :type: float - - :param: iplot, no. of figure window for plotting interims results - :type: integer - - :param: NLLocfile, name and full path to NLLoc-location file - needed when calling class MoMw - :type: string - - :param: picks, dictionary containing picking results - :type: dictionary - - :param: rho [kg/m³], rock density, parameter from autoPyLoT.in - :type: integer - - :param: vp [m/s], P-velocity - :param: integer - - :param: invdir, name and path to inventory or dataless-SEED file - :type: string - ''' - - assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream) - - self._stream = wfstream - self._invdir = invdir - self._t0 = t0 - self._pwin = pwin - self._iplot = iplot - self.setNLLocfile(NLLocfile) - self.setrho(rho) - self.setpicks(picks) - self.setvp(vp) - self.setQp(Qp) - self.calcwapp() - self.calcsourcespec() - self.run_calcMoMw() - - @property - def stream(self): - return self._stream - - @stream.setter - def stream(self, wfstream): - self._stream = wfstream - - @property - def t0(self): - return self._t0 - - @t0.setter - def t0(self, value): - self._t0 = value - - @property - def invdir(self): - return self._invdir - - @invdir.setter - def invdir(self, value): - self._invdir = value - - @property - def pwin(self): - return self._pwin - - @pwin.setter - def pwin(self, value): - self._pwin = value - - @property - def plot_flag(self): - return self.iplot - - @plot_flag.setter - def plot_flag(self, value): - self._iplot = value - - def setNLLocfile(self, NLLocfile): - self.NLLocfile = NLLocfile - - def getNLLocfile(self): - return self.NLLocfile - - def setrho(self, rho): - self.rho = rho - - def getrho(self): - return self.rho - - def setvp(self, vp): - self.vp = vp - - def getvp(self): - return self.vp - - def setQp(self, Qp): - self.Qp = Qp - - def getQp(self): - return self.Qp - - def setpicks(self, picks): - self.picks = picks - - def getpicks(self): - return self.picks - - def getwapp(self): - return self.wapp - - def getw0(self): - return self.w0 - - def getfc(self): - return self.fc - - def get_metadata(self): - return read_metadata(self.invdir) - - def plot(self): - pass - - def getpicdic(self): - return self.picdic - - def calcwapp(self): - self.wapp = None - - def calcsourcespec(self): - self.sourcespek = None - - def run_calcMoMw(self): - self.pickdic = None + def net_magnitude(self): + if self: + return np.median([M["mag"] for M in self.magnitudes.values()]) + return None class RichterMagnitude(Magnitude): - ''' + """ Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- seismograph. Has to be derived from instrument corrected traces! - ''' + """ # poles, zeros and sensitivity of WA seismograph # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) @@ -257,29 +123,23 @@ class RichterMagnitude(Magnitude): 'sensitivity': 1 } - def __init__(self, stream, event, t0, calc_win, verbosity=False): + def __init__(self, stream, event, calc_win, verbosity=False): super(RichterMagnitude, self).__init__(stream, event, verbosity) - self._t0 = t0 self._calc_win = calc_win - @property - def t0(self): - return self._t0 - - @t0.setter - def t0(self, value): - self._t0 = value @property def calc_win(self): return self._calc_win + @calc_win.setter def calc_win(self, value): self._calc_win = value - def peak_to_peak(self, st): + + def peak_to_peak(self, st, t0): # simulate Wood-Anderson response st.simulate(paz_remove=None, paz_simulate=self._paz) @@ -303,7 +163,7 @@ class RichterMagnitude(Magnitude): # get time array th = np.arange(0, len(sqH) * dt, dt) # get maximum peak within pick window - iwin = getsignalwin(th, self.t0 - stime, self.calc_win) + iwin = getsignalwin(th, t0 - stime, self.calc_win) wapp = np.max(sqH[iwin]) if self._verbosity: print("Determined Wood-Anderson peak-to-peak amplitude: {0} " @@ -315,7 +175,7 @@ class RichterMagnitude(Magnitude): f = plt.figure(2) plt.plot(th, sqH) plt.plot(th[iwin], sqH[iwin], 'g') - plt.plot([self.t0, self.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' \ % (st[0].stats.station, wapp)) plt.xlabel('Time [s]') @@ -326,6 +186,7 @@ class RichterMagnitude(Magnitude): return wapp + def get(self): for a in self.arrivals: if a.phase not in 'sS': @@ -334,18 +195,20 @@ class RichterMagnitude(Magnitude): station = pick.waveform_id.station_code wf = select_for_phase(self.stream.select( station=station), a.phase) + if not wf: + print('WARNING: no waveform data found for station {0}'.format( + station)) + continue delta = degrees2kilometers(a.distance) - wapp = self.peak_to_peak(wf) + onset = pick.time + wapp = self.peak_to_peak(wf, onset) # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # use of custom relation functions - mag = np.log10(wapp) + richter_magnitude_scaling(delta) + mag = dict(mag=np.log10(wapp) + richter_magnitude_scaling(delta)) self.magnitudes = (station, mag) return self.magnitudes - def net_magnitude(self): - return np.median([M for M in self.magnitudes.values()]) - class MomentMagnitude(Magnitude): ''' @@ -357,6 +220,29 @@ class MomentMagnitude(Magnitude): corresponding moment magntiude Mw. ''' + def __init__(self, stream, event, vp, Qp, density, verbosity=False): + super(MomentMagnitude, self).__init__(stream, event) + + self._vp = vp + self._Qp = Qp + self._density = density + + + @property + def p_velocity(self): + return self._vp + + + @property + def p_attenuation(self): + return self._Qp + + + @property + def rock_density(self): + return self._density + + def run_calcMoMw(self): picks = self.getpicks() @@ -405,6 +291,32 @@ class MomentMagnitude(Magnitude): self.picdic = picks + def get(self): + for a in self.arrivals: + if a.phase not in 'pP': + continue + pick = a.pick_id.get_referred_object() + station = pick.waveform_id.station_code + wf = select_for_phase(self.stream.select( + station=station), a.phase) + if not wf: + continue + onset = pick.time + distance = degrees2kilometers(a.distance) + azimuth = a.azimuth + incidence = a.takeoff_angle + w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance, azimuth, + incidence, self.p_attenuation) + if w0 is None or fc is None: + print("WARNING: insufficient frequency information") + continue + wf = select_for_phase(wf, "P") + M0, Mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, distance) + mag = dict(w0=w0, fc=fc, M0=M0, mag=Mw) + self.magnitudes = (station, mag) + return self.magnitudes + + def calc_woodanderson_pp(st, metadata, T0, win=10., verbosity=False): if verbosity: print ("Getting Wood-Anderson peak-to-peak amplitude ...") @@ -491,8 +403,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta): return Mo, Mw -def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, - qp, iplot): +def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, + qp, iplot=0): ''' 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 @@ -500,16 +412,12 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, thus restitution and integration necessary! Integrated traces are rotated into ray-coordinate system ZNE => LQT using Obspy's rotate modul! - :param: wfstream + :param: wfstream (corrected for instrument) :type: `~obspy.core.stream.Stream` :param: onset, P-phase onset time :type: float - :param: metadata, tuple or list containing type of inventory and either - list of files or inventory object - :type: tuple or list - :param: vp, Vp-wave velocity :type: float @@ -533,163 +441,151 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, # get Q value Q, A = qp - delta = delta * 1000 # hypocentral distance in [m] + dist = delta * 1000 # hypocentral distance in [m] fc = None w0 = None - wf_copy = wfstream.copy() - invtype, inventory = metadata + zdat = select_for_phase(wfstream, "P") - [cordat, restflag] = restitute_data(wf_copy, invtype, inventory) - if restflag is True: - zdat = cordat.select(component="Z") - if len(zdat) == 0: - zdat = cordat.select(component="3") - cordat_copy = cordat.copy() - # get equal time stamps and lengths of traces - # necessary for rotation of traces - trstart, trend = common_range(cordat_copy) - cordat_copy.trim(trstart, trend) - try: - # rotate into LQT (ray-coordindate-) system using Obspy's rotate - # L: P-wave direction - # Q: SV-wave direction - # T: SH-wave direction - LQT = cordat_copy.rotate('ZNE->LQT', azimuth, incidence) - ldat = LQT.select(component="L") - if len(ldat) == 0: - # if horizontal channels are 2 and 3 - # no azimuth information is available and thus no - # rotation is possible! - print("calcsourcespec: Azimuth information is missing, " - "no rotation of components possible!") - ldat = LQT.select(component="Z") + dt = zdat[0].stats.delta - # integrate to displacement - # unrotated vertical component (for copmarison) - inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None, - zdat[0].stats.delta)) - # rotated component Z => L - Ldat = signal.detrend(integrate.cumtrapz(ldat[0].data, None, - ldat[0].stats.delta)) + freq = zdat[0].stats.sampling_rate - # get window after P pulse for - # calculating source spectrum - tstart = UTCDateTime(zdat[0].stats.starttime) - tonset = onset.timestamp - tstart.timestamp - impickP = tonset * zdat[0].stats.sampling_rate - wfzc = Ldat[impickP: len(Ldat) - 1] - # get time array - t = np.arange(0, len(inttrz) * zdat[0].stats.delta, \ - zdat[0].stats.delta) - # calculate spectrum using only first cycles of - # waveform after P onset! - zc = crossings_nonzero_all(wfzc) - if np.size(zc) == 0 or len(zc) <= 3: - print ("calcsourcespec: Something is wrong with the waveform, " - "no zero crossings derived!") - print ("No calculation of source spectrum possible!") - plotflag = 0 - else: - plotflag = 1 - index = min([3, len(zc) - 1]) - calcwin = (zc[index] - zc[0]) * zdat[0].stats.delta - iwin = getsignalwin(t, tonset, calcwin) - xdat = Ldat[iwin] + # trim traces to common range (for rotation) + trstart, trend = common_range(wfstream) + wfstream.trim(trstart, trend) - # fft - fny = zdat[0].stats.sampling_rate / 2 - l = len(xdat) / zdat[0].stats.sampling_rate - # number of fft bins after Bath - n = zdat[0].stats.sampling_rate * l - # 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 = zdat[0].stats.delta * np.fft.fft(xdat, N) - Y = abs(y[: N / 2]) - L = (N - 1) / zdat[0].stats.sampling_rate - f = np.arange(0, fny, 1 / L) + # rotate into LQT (ray-coordindate-) system using Obspy's rotate + # L: P-wave direction + # Q: SV-wave direction + # T: SH-wave direction + LQT = wfstream.rotate('ZNE->LQT', azimuth, incidence) + ldat = LQT.select(component="L") + if len(ldat) == 0: + # if horizontal channels are 2 and 3 + # no azimuth information is available and thus no + # rotation is possible! + print("calcsourcespec: Azimuth information is missing, " + "no rotation of components possible!") + ldat = LQT.select(component="Z") - # remove zero-frequency and frequencies above - # corner frequency of seismometer (assumed - # to be 100 Hz) - fi = np.where((f >= 1) & (f < 100)) - F = f[fi] - YY = Y[fi] + # integrate to displacement + # unrotated vertical component (for comparison) + inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None, dt)) - # correction for attenuation - wa = 2 * np.pi * F # angular frequency - D = np.exp((wa * delta) / (2 * vp * Q * F ** A)) - YYcor = YY.real * D + # rotated component Z => L + Ldat = signal.detrend(integrate.cumtrapz(ldat[0].data, None, dt)) - # get plateau (DC value) and corner frequency - # initial guess of plateau - w0in = np.mean(YYcor[0:100]) - # initial guess of corner frequency - # where spectral level reached 50% of flat level - iin = np.where(YYcor >= 0.5 * w0in) - Fcin = F[iin[0][np.size(iin) - 1]] + # get window after P pulse for + # calculating source spectrum + rel_onset = onset - trstart + impickP = int(rel_onset * freq) + wfzc = Ldat[impickP: len(Ldat) - 1] + # get time array + t = np.arange(0, len(inttrz) * dt, dt) + # calculate spectrum using only first cycles of + # waveform after P onset! + zc = crossings_nonzero_all(wfzc) + if np.size(zc) == 0 or len(zc) <= 3: + print ("calcsourcespec: Something is wrong with the waveform, " + "no zero crossings derived!") + print ("No calculation of source spectrum possible!") + plotflag = 0 + else: + plotflag = 1 + index = min([3, len(zc) - 1]) + calcwin = (zc[index] - zc[0]) * dt + iwin = getsignalwin(t, rel_onset, calcwin) + xdat = Ldat[iwin] - # use of implicit scipy otimization function - fit = synthsourcespec(F, w0in, Fcin) - [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, - Fcin]) - w01 = optspecfit[0] - fc1 = optspecfit[1] - print ("calcsourcespec: Determined w0-value: %e m/Hz, \n" - "Determined corner frequency: %f Hz" % (w01, fc1)) + # fft + fny = freq / 2 + l = len(xdat) / freq + # number of fft bins after Bath + n = freq * l + # 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 = dt * np.fft.fft(xdat, N) + Y = abs(y[: N / 2]) + L = (N - 1) / freq + f = np.arange(0, fny, 1 / L) - # use of conventional fitting - [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot) + # remove zero-frequency and frequencies above + # corner frequency of seismometer (assumed + # to be 100 Hz) + fi = np.where((f >= 1) & (f < 100)) + F = f[fi] + YY = Y[fi] - # get w0 and fc as median of both - # source spectrum fits - w0 = np.median([w01, w02]) - fc = np.median([fc1, fc2]) - print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) + # correction for attenuation + wa = 2 * np.pi * F # angular frequency + D = np.exp((wa * dist) / (2 * vp * Q * F ** A)) + YYcor = YY.real * D - except TypeError as er: - raise TypeError('''{0}'''.format(er)) + # get plateau (DC value) and corner frequency + # initial guess of plateau + w0in = np.mean(YYcor[0:100]) + # initial guess of corner frequency + # where spectral level reached 50% of flat level + iin = np.where(YYcor >= 0.5 * w0in) + Fcin = F[iin[0][np.size(iin) - 1]] - if iplot > 1: - f1 = plt.figure() - tLdat = np.arange(0, len(Ldat) * zdat[0].stats.delta, \ - zdat[0].stats.delta) - plt.subplot(2, 1, 1) - # show displacement in mm - p1, = plt.plot(t, np.multiply(inttrz, 1000), 'k') - p2, = plt.plot(tLdat, np.multiply(Ldat, 1000)) - plt.legend([p1, p2], ['Displacement', 'Rotated Displacement']) - if plotflag == 1: - plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') - plt.title('Seismogram and P Pulse, Station %s-%s' \ - % (zdat[0].stats.station, zdat[0].stats.channel)) - else: - plt.title('Seismogram, Station %s-%s' \ - % (zdat[0].stats.station, zdat[0].stats.channel)) - plt.xlabel('Time since %s' % zdat[0].stats.starttime) - plt.ylabel('Displacement [mm]') + # use of implicit scipy otimization function + fit = synthsourcespec(F, w0in, Fcin) + [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin]) + w01 = optspecfit[0] + fc1 = optspecfit[1] + print ("calcsourcespec: Determined w0-value: %e m/Hz, \n" + "Determined corner frequency: %f Hz" % (w01, fc1)) - if plotflag == 1: - plt.subplot(2, 1, 2) - p1, = plt.loglog(f, Y.real, 'k') - p2, = plt.loglog(F, YY.real) - p3, = plt.loglog(F, YYcor, 'r') - p4, = plt.loglog(F, fit, 'g') - plt.loglog([fc, fc], [w0 / 100, w0], 'g') - plt.legend([p1, p2, p3, p4], ['Raw Spectrum', \ - 'Used Raw Spectrum', \ - 'Q-Corrected Spectrum', \ - 'Fit to Spectrum']) - plt.title('Source Spectrum from P Pulse, w0=%e m/Hz, fc=%6.2f Hz' \ - % (w0, fc)) - plt.xlabel('Frequency [Hz]') - plt.ylabel('Amplitude [m/Hz]') - plt.grid() - plt.show() - raw_input() - plt.close(f1) + # use of conventional fitting + [w02, fc2] = fitSourceModel(F, YYcor, Fcin, iplot) + + # get w0 and fc as median of both + # source spectrum fits + w0 = np.median([w01, w02]) + fc = np.median([fc1, fc2]) + print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) + + if iplot > 1: + f1 = plt.figure() + tLdat = np.arange(0, len(Ldat) * dt, dt) + plt.subplot(2, 1, 1) + # show displacement in mm + p1, = plt.plot(t, np.multiply(inttrz, 1000), 'k') + p2, = plt.plot(tLdat, np.multiply(Ldat, 1000)) + plt.legend([p1, p2], ['Displacement', 'Rotated Displacement']) + if plotflag == 1: + plt.plot(t[iwin], np.multiply(xdat, 1000), 'g') + plt.title('Seismogram and P Pulse, Station %s-%s' \ + % (zdat[0].stats.station, zdat[0].stats.channel)) + else: + plt.title('Seismogram, Station %s-%s' \ + % (zdat[0].stats.station, zdat[0].stats.channel)) + plt.xlabel('Time since %s' % zdat[0].stats.starttime) + plt.ylabel('Displacement [mm]') + + if plotflag == 1: + plt.subplot(2, 1, 2) + p1, = plt.loglog(f, Y.real, 'k') + p2, = plt.loglog(F, YY.real) + p3, = plt.loglog(F, YYcor, 'r') + p4, = plt.loglog(F, fit, 'g') + plt.loglog([fc, fc], [w0 / 100, w0], 'g') + plt.legend([p1, p2, p3, p4], ['Raw Spectrum', \ + 'Used Raw Spectrum', \ + 'Q-Corrected Spectrum', \ + 'Fit to Spectrum']) + plt.title('Source Spectrum from P Pulse, w0=%e m/Hz, fc=%6.2f Hz' \ + % (w0, fc)) + plt.xlabel('Frequency [Hz]') + plt.ylabel('Amplitude [m/Hz]') + plt.grid() + plt.show() + raw_input() + plt.close(f1) return w0, fc @@ -847,9 +743,17 @@ def calc_moment_magnitude(e, wf, metadata, vp, Qp, rho): continue onset = pick.time dist = degrees2kilometers(a.distance) - w0, fc = calcsourcespec(wf, onset, metadata, vp, dist, a.azimuth, a.takeoff_angle, Qp, 0) + invtype, inventory = metadata + [corr_wf, rest_flag] = restitute_data(wf, invtype, inventory) + if not rest_flag: + print("WARNING: data for {0} could not be corrected".format( + station)) + continue + w0, fc = calcsourcespec(corr_wf, onset, vp, dist, a.azimuth, + a.takeoff_angle, Qp, 0) if w0 is None or fc is None: continue + wf = select_for_phase(corr_wf, "P") station_mag = calcMoMw(wf, w0, rho, vp, dist) mags[station] = station_mag mag = np.median([M[1] for M in mags.values()]) From 28a5cedbc6bf67e412bdb037f75f0e42e80176ce Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Tue, 27 Sep 2016 15:12:14 +0200 Subject: [PATCH 15/32] [refactor] further refactoring done -> obsolete functions deleted, imports optimized, output suppressed and calculation done in __init__ --- pylot/core/analysis/magnitude.py | 212 +++---------------------------- 1 file changed, 19 insertions(+), 193 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 48417557..ed338a41 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -6,24 +6,23 @@ Created autumn/winter 2015. :author: Ludger Küperkoch / MAGS2 EP3 working group """ import os + import matplotlib.pyplot as plt import numpy as np -from obspy.core import Stream, UTCDateTime from obspy.core.event import Magnitude as ObspyMagnitude 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 getPatternLine -from scipy.optimize import curve_fit -from scipy import integrate, signal -from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.utils import common_range, fit_curve + def richter_magnitude_scaling(delta): relation = np.loadtxt(os.path.join(os.path.expanduser('~'), - '.pylot', 'richter_scaling.data')) + '.pylot', 'richter_scaling.data')) # prepare spline interpolation to calculate return value - func, params = fit_curve(relation[:,0], relation[:, 1]) + func, params = fit_curve(relation[:, 0], relation[:, 1]) return func(delta, params) @@ -33,58 +32,53 @@ class Magnitude(object): """ def __init__(self, stream, event, verbosity=False): + self._type = "M" self._plot_flag = False self._verbosity = verbosity self._event = event self._stream = stream self._magnitudes = dict() - def __str__(self): print('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)) - def __nonzero__(self): return bool(self.magnitudes) + @property + def type(self): + return self._type @property def plot_flag(self): return self._plot_flag - @plot_flag.setter def plot_flag(self, value): self._plot_flag = value - @property def stream(self): return self._stream - @stream.setter def stream(self, value): self._stream = value - @property def event(self): return self._event - @property def arrivals(self): return self._event.origins[0].arrivals - @property def magnitudes(self): return self._magnitudes - @magnitudes.setter def magnitudes(self, value): """ @@ -97,14 +91,12 @@ class Magnitude(object): station, magnitude = value self._magnitudes[station] = magnitude - - def get(self): - return self.magnitudes - + def calc(self): + pass def net_magnitude(self): if self: - return np.median([M["mag"] for M in self.magnitudes.values()]) + return ObspyMagnitude(mag=np.median([M["mag"] for M in self.magnitudes.values()]), type=self.type) return None @@ -127,18 +119,17 @@ class RichterMagnitude(Magnitude): super(RichterMagnitude, self).__init__(stream, event, verbosity) self._calc_win = calc_win - + self._type = 'ML' + self.calc() @property def calc_win(self): return self._calc_win - @calc_win.setter def calc_win(self, value): self._calc_win = value - def peak_to_peak(self, st, t0): # simulate Wood-Anderson response @@ -186,8 +177,7 @@ class RichterMagnitude(Magnitude): return wapp - - def get(self): + def calc(self): for a in self.arrivals: if a.phase not in 'sS': continue @@ -207,7 +197,6 @@ class RichterMagnitude(Magnitude): # use of custom relation functions mag = dict(mag=np.log10(wapp) + richter_magnitude_scaling(delta)) self.magnitudes = (station, mag) - return self.magnitudes class MomentMagnitude(Magnitude): @@ -226,72 +215,22 @@ class MomentMagnitude(Magnitude): self._vp = vp self._Qp = Qp self._density = density - + self._type = 'Mw' + self.calc() @property def p_velocity(self): return self._vp - @property def p_attenuation(self): return self._Qp - @property def rock_density(self): return self._density - - def run_calcMoMw(self): - - picks = self.getpicks() - nllocfile = self.getNLLocfile() - wfdat = self.stream - self.picdic = None - - for key in picks: - if picks[key]['P']['weight'] < 4: - # select waveform - selwf = wfdat.select(station=key) - if len(key) > 4: - Ppattern = '%s ? ? ? P' % key - elif len(key) == 4: - Ppattern = '%s ? ? ? P' % key - elif len(key) < 4: - Ppattern = '%s ? ? ? P' % key - nllocline = getPatternLine(nllocfile, Ppattern) - # get hypocentral distance, station azimuth and - # angle of incidence from NLLoc-location file - delta = float(nllocline.split(None)[21]) - az = float(nllocline.split(None)[22]) - inc = float(nllocline.split(None)[24]) - # call subfunction to estimate source spectrum - # and to derive w0 and fc - [w0, fc] = calcsourcespec(selwf, picks[key]['P']['mpp'], \ - self.get_metadata(), self.getvp(), delta, az, \ - inc, self.getQp(), self.plot_flag) - - if w0 is not None: - # call subfunction to calculate Mo and Mw - zdat = selwf.select(component="Z") - if len(zdat) == 0: # check for other components - zdat = selwf.select(component="3") - [Mo, Mw] = calcMoMw(zdat, w0, self.getrho(), self.getvp(), - delta) - else: - Mo = None - Mw = None - - # add w0, fc, Mo and Mw to dictionary - picks[key]['P']['w0'] = w0 - picks[key]['P']['fc'] = fc - picks[key]['P']['Mo'] = Mo - picks[key]['P']['Mw'] = Mw - self.picdic = picks - - - def get(self): + def calc(self): for a in self.arrivals: if a.phase not in 'pP': continue @@ -314,52 +253,6 @@ class MomentMagnitude(Magnitude): M0, Mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, distance) mag = dict(w0=w0, fc=fc, M0=M0, mag=Mw) self.magnitudes = (station, mag) - return self.magnitudes - - -def calc_woodanderson_pp(st, metadata, T0, win=10., verbosity=False): - if verbosity: - print ("Getting Wood-Anderson peak-to-peak amplitude ...") - print ("Simulating Wood-Anderson seismograph ...") - - # poles, zeros and sensitivity of WA seismograph - # (see Uhrhammer & Collins, 1990, BSSA, pp. 702-716) - paz_wa = { - 'poles': [5.6089 - 5.4978j, -5.6089 - 5.4978j], - 'zeros': [0j, 0j], - 'gain': 2080, - 'sensitivity': 1 - } - - stime, etime = common_range(st) - st.trim(stime, etime) - - invtype, inventory = metadata - - st = restitute_data(st, invtype, inventory) - - dt = st[0].stats.delta - - st.simulate(paz_remove=None, paz_simulate=paz_wa) - - # get RMS of both horizontal components - power = [np.power(tr.data, 2) for tr in st if tr.stats.channel[-1] not - in 'Z3'] - if len(power) != 2: - raise ValueError('Wood-Anderson amplitude defintion only valid for ' - 'two horizontals: {0} given'.format(len(power))) - power_sum = power[0] + power[1] - sqH = np.sqrt(power_sum) - - # get time array - th = np.arange(0, len(sqH) * dt, dt) - # get maximum peak within pick window - iwin = getsignalwin(th, T0 - stime, win) - wapp = np.max(sqH[iwin]) - if verbosity: - print("Determined Wood-Anderson peak-to-peak amplitude: {0} " - "mm".format(wapp)) - return wapp def calcMoMw(wfstream, w0, rho, vp, delta): @@ -697,70 +590,3 @@ def fitSourceModel(f, S, fc0, iplot): plt.close() return w0, fc - - -def calc_richter_magnitude(e, wf, metadata, amp_win): - if e.origins: - o = e.origins[0] - mags = dict() - for a in o.arrivals: - if a.phase not in 'sS': - continue - pick = a.pick_id.get_referred_object() - station = pick.waveform_id.station_code - wf = select_for_phase(wf.select( - station=station), a.phase) - delta = degrees2kilometers(a.distance) - wapp = calc_woodanderson_pp(wf, metadata, pick.time, - amp_win) - # using standard Gutenberg-Richter relation - # TODO make the ML calculation more flexible by allowing - # use of custom relation functions - mag = np.log10(wapp) + richter_magnitude_scaling(delta) - mags[station] = mag - mag = np.median([M for M in mags.values()]) - # give some information on the processing - print('number of stations used: {0}\n'.format(len(mags.values()))) - print('stations used:\n') - for s in mags.keys(): print('\t{0}'.format(s)) - - return ObspyMagnitude(mag=mag, magnitude_type='ML') - else: - return None - - -def calc_moment_magnitude(e, wf, metadata, vp, Qp, rho): - if e.origins: - o = e.origins[0] - mags = dict() - for a in o.arrivals: - if a.phase in 'sS': - continue - pick = a.pick_id.get_referred_object() - station = pick.waveform_id.station_code - wf = wf.select(station=station) - if not wf: - continue - onset = pick.time - dist = degrees2kilometers(a.distance) - invtype, inventory = metadata - [corr_wf, rest_flag] = restitute_data(wf, invtype, inventory) - if not rest_flag: - print("WARNING: data for {0} could not be corrected".format( - station)) - continue - w0, fc = calcsourcespec(corr_wf, onset, vp, dist, a.azimuth, - a.takeoff_angle, Qp, 0) - if w0 is None or fc is None: - continue - wf = select_for_phase(corr_wf, "P") - station_mag = calcMoMw(wf, w0, rho, vp, dist) - mags[station] = station_mag - mag = np.median([M[1] for M in mags.values()]) - # give some information on the processing - print('number of stations used: {0}\n'.format(len(mags.values()))) - print('stations used:\n') - for s in mags.keys(): print('\t{0}'.format(s)) - return ObspyMagnitude(mag=mag, magnitude_type='Mw') - else: - return None From cf514ae024b596f852c2004e3c8700d7e9e7d841 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Tue, 27 Sep 2016 15:13:51 +0200 Subject: [PATCH 16/32] [change] traces that could not be restituted are now removed from trace --- pylot/core/util/dataprocessing.py | 12 +++++------- 1 file changed, 5 insertions(+), 7 deletions(-) diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index 7ec89e7d..a0ba4bcf 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -214,10 +214,6 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False): for tr in data: seed_id = tr.get_id() # check, whether this trace has already been corrected - # TODO read actual value of processing key in Trace's stats instead - # of just looking for thr key: if processing is setit doesn't - # necessarily mean that the trace has been corrected so far but only - # processed in some way, e.g. normalized if 'processing' in tr.stats.keys() \ and np.all(['remove' in p for p in tr.stats.processing]) \ and not force: @@ -252,7 +248,7 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False): finv = invlist[0] inventory = read_inventory(finv, format='STATIONXML') else: - restflag.append(False) + data.remove(tr) continue # apply restitution to data try: @@ -268,7 +264,9 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False): if msg0 not in e.message or msg1 not in e.message: raise else: - restflag.append(False) + # restitution done to copies of data thus deleting traces + # that failed should not be a problem + data.remove(tr) continue restflag.append(True) # check if ALL traces could be restituted, take care of large datasets @@ -319,4 +317,4 @@ def get_prefilt(trace, tlow=(0.5, 0.9), thi=(5., 2.), verbosity=0): if __name__ == "__main__": import doctest - doctest.testmod() \ No newline at end of file + doctest.testmod() From 699ba6f12298e40f8981976dfc8641cc9d12b872 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Tue, 27 Sep 2016 15:14:48 +0200 Subject: [PATCH 17/32] [new] added a new Error -> ProcessingError raised in case of failed restitution --- pylot/core/util/errors.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/pylot/core/util/errors.py b/pylot/core/util/errors.py index a6015e55..8805b488 100644 --- a/pylot/core/util/errors.py +++ b/pylot/core/util/errors.py @@ -24,3 +24,6 @@ class OverwriteError(IOError): class ParameterError(Exception): pass + +class ProcessingError(RuntimeError): + pass From 72d15e1fc5cafbaf6d9f9f4de93e848ac166ee42 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Tue, 27 Sep 2016 15:15:53 +0200 Subject: [PATCH 18/32] [new] implemented new magnitude concept into QtPyLoT --- QtPyLoT.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 84d32732..a2148fe3 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -40,7 +40,7 @@ from PySide.QtGui import QMainWindow, QInputDialog, QIcon, QFileDialog, \ import numpy as np from obspy import UTCDateTime -from pylot.core.analysis.magnitude import calc_richter_magnitude, calc_moment_magnitude +from pylot.core.analysis.magnitude import RichterMagnitude, MomentMagnitude from pylot.core.io.data import Data from pylot.core.io.inputs import FilterOptions, AutoPickParameter from pylot.core.pick.autopick import autopickevent @@ -51,9 +51,9 @@ import pylot.core.loc.nll as nll from pylot.core.util.defaults import FILTERDEFAULTS, COMPNAME_MAP, \ AUTOMATIC_DEFAULTS from pylot.core.util.errors import FormatError, DatastructureError, \ - OverwriteError + OverwriteError, ProcessingError from pylot.core.util.connection import checkurl -from pylot.core.util.dataprocessing import read_metadata +from pylot.core.util.dataprocessing import read_metadata, restitute_data from pylot.core.util.utils import fnConstructor, getLogin, \ full_range from pylot.core.io.location import create_creation_info, create_event @@ -1003,10 +1003,16 @@ class MainWindow(QMainWindow): settings.setValue("inventoryFile", fninv) settings.sync() self.metadata = read_metadata(fninv) + wf_copy = self.get_data().getWFData().copy() + [corr_wf, rest_flag] = restitute_data(wf_copy, *self.metadata) + if not rest_flag: + raise ProcessingError('Restitution of waveform data failed!') if type == 'ML': - return calc_richter_magnitude(self.get_data().get_evt_data(), self.get_data().getWFData(), self.metadata, self.inputs.get('sstop')) + local_mag = RichterMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('sstop')) + return local_mag.net_magnitude() elif type == 'Mw': - return calc_moment_magnitude(self.get_data().get_evt_data(), self.get_data().getWFData(), self.metadata, self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho')) + moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho')) + return moment_mag.net_magnitude() else: return None From c1bddd5c0be3de49bf020101704f5dcb457d6738 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 10:56:05 +0200 Subject: [PATCH 19/32] [change] improved verbosity and plotting control for Magnitude objects --- pylot/core/analysis/magnitude.py | 80 +++++++++++++++++++------------- 1 file changed, 49 insertions(+), 31 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index ed338a41..085cbe8c 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -31,9 +31,9 @@ class Magnitude(object): 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._plot_flag = False + self._plot_flag = iplot self._verbosity = verbosity self._event = event self._stream = stream @@ -59,6 +59,17 @@ class Magnitude(object): def plot_flag(self, 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 def stream(self): return self._stream @@ -115,8 +126,8 @@ class RichterMagnitude(Magnitude): 'sensitivity': 1 } - def __init__(self, stream, event, calc_win, verbosity=False): - super(RichterMagnitude, self).__init__(stream, event, verbosity) + def __init__(self, stream, event, calc_win, verbosity=False, iplot=0): + super(RichterMagnitude, self).__init__(stream, event, verbosity, iplot) self._calc_win = calc_win self._type = 'ML' @@ -156,7 +167,7 @@ class RichterMagnitude(Magnitude): # get maximum peak within pick window iwin = getsignalwin(th, t0 - stime, self.calc_win) wapp = np.max(sqH[iwin]) - if self._verbosity: + if self.verbose: print("Determined Wood-Anderson peak-to-peak amplitude: {0} " "mm".format(wapp)) @@ -186,8 +197,9 @@ class RichterMagnitude(Magnitude): wf = select_for_phase(self.stream.select( station=station), a.phase) if not wf: - print('WARNING: no waveform data found for station {0}'.format( - station)) + if self.verbose: + print('WARNING: no waveform data found for station {0}'.format( + station)) continue delta = degrees2kilometers(a.distance) onset = pick.time @@ -209,8 +221,8 @@ class MomentMagnitude(Magnitude): corresponding moment magntiude Mw. ''' - def __init__(self, stream, event, vp, Qp, density, verbosity=False): - super(MomentMagnitude, self).__init__(stream, event) + def __init__(self, stream, event, vp, Qp, density, verbosity=False, iplot=False): + super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot) self._vp = vp self._Qp = Qp @@ -245,17 +257,18 @@ class MomentMagnitude(Magnitude): azimuth = a.azimuth incidence = a.takeoff_angle 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: - print("WARNING: insufficient frequency information") + if self.verbose: + print("WARNING: insufficient frequency information") continue 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) 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 seismic moments and corresponding moment magnitudes. @@ -279,8 +292,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta): tr = wfstream[0] delta = delta * 1000 # hypocentral distance in [m] - print("calcMoMw: Calculating seismic moment Mo and moment magnitude Mw for station %s ..." \ - % tr.stats.station) + if verbosity: + print("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(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) * 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 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 (usually called omega0) and the corner frequency assuming Aki's omega-square @@ -329,7 +343,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, :param: iplot, show results (iplot>1) or not (iplot(<2) :type: integer ''' - print ("Calculating source spectrum ....") + if verbosity: + print ("Calculating source spectrum ....") # get Q value Q, A = qp @@ -359,8 +374,9 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, # if horizontal channels are 2 and 3 # no azimuth information is available and thus no # rotation is possible! - print("calcsourcespec: Azimuth information is missing, " - "no rotation of components possible!") + if verbosity: + print("calcsourcespec: Azimuth information is missing, " + "no rotation of components possible!") ldat = LQT.select(component="Z") # integrate to displacement @@ -381,9 +397,10 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, # waveform after P onset! zc = crossings_nonzero_all(wfzc) if np.size(zc) == 0 or len(zc) <= 3: - print ("calcsourcespec: Something is wrong with the waveform, " - "no zero crossings derived!") - print ("No calculation of source spectrum possible!") + if verbosity: + 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 @@ -430,17 +447,19 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, [optspecfit, _] = curve_fit(synthsourcespec, F, YYcor, [w0in, Fcin]) w01 = optspecfit[0] fc1 = optspecfit[1] - print ("calcsourcespec: Determined w0-value: %e m/Hz, \n" - "Determined corner frequency: %f Hz" % (w01, fc1)) + if verbosity: + print ("calcsourcespec: Determined w0-value: %e m/Hz, \n" + "Determined corner frequency: %f Hz" % (w01, fc1)) # 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 # source spectrum fits w0 = np.median([w01, w02]) fc = np.median([fc1, fc2]) - print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) + if verbosity: + print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % (w0, fc)) if iplot > 1: f1 = plt.figure() @@ -504,7 +523,7 @@ def synthsourcespec(f, omega0, fcorner): 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. 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: fc = fc0 w0 = max(S) - - print("fitSourceModel: best fc: %fHz, best w0: %e m/Hz" \ - % (fc, w0)) + if verbosity: + print("fitSourceModel: best fc: {0} Hz, best w0: {1} m/Hz".format(fc, w0)) if iplot > 1: plt.figure(iplot) From ae967b3429d14769c8f183e4330fd0d42e4c8650 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 10:57:52 +0200 Subject: [PATCH 20/32] [remove] removed Wood-Anderson peak-to-peak amplitude reading from autopick.py; newly introduced in autoPyLoT in a future commit --- pylot/core/pick/autopick.py | 32 -------------------------------- 1 file changed, 32 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 5a017adc..35111071 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -17,10 +17,8 @@ from pylot.core.pick.charfuns import CharacteristicFunction from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \ getSNR, fmpicker, checkPonsets, wadaticheck -from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.utils import getPatternLine from pylot.core.io.data import Data -from pylot.core.analysis.magnitude import RichterMagnitude def autopickevent(data, param): @@ -121,8 +119,6 @@ def autopickstation(wfstream, pickparam, verbose=False): # parameter to check for spuriously picked S onset zfac = pickparam.get('zfac') # path to inventory-, dataless- or resp-files - invdir = pickparam.get('invdir') - invtype, inventory = read_metadata(invdir) # initialize output Pweight = 4 # weight for P onset @@ -565,24 +561,6 @@ def autopickstation(wfstream, pickparam, verbose=False): # re-create stream object including both horizontal components hdat = edat.copy() hdat += ndat - h_copy = hdat.copy() - [cordat, restflag] = restitute_data(h_copy, invtype, inventory) - # calculate WA-peak-to-peak amplitude - # using subclass WApp of superclass Magnitude - if restflag: - if Sweight < 4: - wapp = RichterMagnitude(cordat, mpickS, mpickP + sstop, iplot) - else: - # use larger window for getting peak-to-peak amplitude - # as the S pick is quite unsure - wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop + - (0.5 * (mpickP + sstop)), iplot) - - Ao = wapp.getwapp() - else: - print("Go on processing data without source parameter " - "determination!") - else: msg = 'Bad initial (AIC) S-pick, skipping this onset!\n' \ 'AIC-SNR={0}, AIC-Slope={1}counts/s\n' \ @@ -602,16 +580,6 @@ def autopickstation(wfstream, pickparam, verbose=False): # re-create stream object including both horizontal components hdat = edat.copy() hdat += ndat - h_copy = hdat.copy() - [cordat, restflag] = restitute_data(h_copy, invtype, inventory) - if restflag is True: - # calculate WA-peak-to-peak amplitude - # using subclass WApp of superclass Magnitude - wapp = RichterMagnitude(cordat, mpickP, mpickP + sstop + (0.5 * (mpickP - + sstop)), - iplot) - Ao = wapp.getwapp() - else: print('autopickstation: No horizontal component data available or ' \ 'bad P onset, skipping S picking!') From 4e520df145e20af55c6e3e033b501ea67cc5bb80 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 10:59:50 +0200 Subject: [PATCH 21/32] [new] added Wood-Anderson amplitude output for further analysis --- pylot/core/analysis/magnitude.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 085cbe8c..994c2230 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -126,6 +126,8 @@ class RichterMagnitude(Magnitude): 'sensitivity': 1 } + _amplitudes = dict() + def __init__(self, stream, event, calc_win, verbosity=False, iplot=0): super(RichterMagnitude, self).__init__(stream, event, verbosity, iplot) @@ -141,6 +143,15 @@ class RichterMagnitude(Magnitude): def calc_win(self, value): self._calc_win = value + @property + def amplitudes(self): + return self._amplitudes + + @amplitudes.setter + def amplitudes(self, value): + station, a0 = value + self._amplitudes[station] = a0 + def peak_to_peak(self, st, t0): # simulate Wood-Anderson response @@ -203,11 +214,12 @@ class RichterMagnitude(Magnitude): continue delta = degrees2kilometers(a.distance) onset = pick.time - wapp = self.peak_to_peak(wf, onset) + a0 = self.peak_to_peak(wf, onset) + self.amplitudes = (station, a0) # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # use of custom relation functions - mag = dict(mag=np.log10(wapp) + richter_magnitude_scaling(delta)) + mag = dict(mag=np.log10(a0) + richter_magnitude_scaling(delta)) self.magnitudes = (station, mag) From 231e7dafa93d00ade36a222109c604db79a2cca5 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 11:01:09 +0200 Subject: [PATCH 22/32] [new] added a function to easily add amplitude information to a given Obspy event object --- pylot/core/io/phases.py | 20 +++++++++++++++++++- 1 file changed, 19 insertions(+), 1 deletion(-) diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py index d3064f6f..7871f4bb 100644 --- a/pylot/core/io/phases.py +++ b/pylot/core/io/phases.py @@ -15,6 +15,24 @@ from pylot.core.pick.utils import select_for_phase from pylot.core.util.utils import getOwner, full_range, four_digits +def add_amplitudes(event, amplitudes): + amplitude_list = [] + for pick in event.picks: + try: + a0 = amplitudes[pick.waveform_id.station_code] + amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3) + amplitude.unit = 'm' + amplitude.category = 'point' + amplitude.waveform_id = pick.waveform_id + amplitude.magnitude_hint = 'ML' + amplitude.pick_id = pick.resource_id + amplitude.type = 'AML' + amplitude_list.append(amplitude) + except KeyError: + continue + event.amplitudes = amplitude_list + return event + def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs): """ readPILOTEvent - function @@ -549,4 +567,4 @@ def merge_picks(event, picks): if p.waveform_id.station_code == station and p.phase_hint == phase: p.time, p.time_errors = time, err del time, err, phase, station - return event \ No newline at end of file + return event From 5f92d1f0db0038c3df4d8312c465c5fcf1f570f6 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 11:03:53 +0200 Subject: [PATCH 23/32] [refactor] major refactoring of autoPyLoT and making use of the newly introduced Magnitude objects --- autoPyLoT.py | 213 ++++++++++++++------------------------------------- 1 file changed, 56 insertions(+), 157 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 21d7eb52..2bf14b76 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -5,17 +5,18 @@ from __future__ import print_function import argparse import glob -import string import os -import numpy as np +from obspy import read_events -from pylot.core.analysis.magnitude import MomentMagnitude +import pylot.core.loc.hsat as hsat +import pylot.core.loc.nll as nll +from pylot.core.analysis.magnitude import MomentMagnitude, RichterMagnitude from pylot.core.io.data import Data from pylot.core.io.inputs import AutoPickParameter -import pylot.core.loc.nll as nll -import pylot.core.loc.hsat as hsat +from pylot.core.io.phases import add_amplitudes from pylot.core.pick.autopick import autopickevent, iteratepicker +from pylot.core.util.dataprocessing import restitute_data, read_metadata from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.version import get_git_version as _getVersionString @@ -54,6 +55,7 @@ def autoPyLoT(inputfile): data = Data() + evt = None # getting information on data structure if parameter.hasParam('datastructure'): @@ -97,139 +99,25 @@ def autoPyLoT(inputfile): print("!!No source parameter estimation possible!!") print(" !!! ") - # multiple event processing - # read each event in database datapath = datastructure.expandDataPath() if not parameter.hasParam('eventID'): - for event in [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)]: - data.setWFData(glob.glob(os.path.join(datapath, event, '*'))) - print('Working on event %s' % event) - print(data) - wfdat = data.getWFData() # all available streams - ########################################################## - # !automated picking starts here! - picks = autopickevent(wfdat, parameter) - finalpicks = picks - ########################################################## - # locating - if locflag == 1: - # write phases to NLLoc-phase file - nll.export(picks, phasefile) - - # For locating the event the NLLoc-control file has to be modified! - evID = event[string.rfind(event, "/") + 1: len(events) - 1] - nllocout = '%s_%s' % (evID, nllocoutpatter) - # create comment line for NLLoc-control file - nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, - ttpat) - - # locate the event - nll.locate(ctrfile) - - # !iterative picking if traces remained unpicked or occupied with bad picks! - # get theoretical onset times for picks with weights >= 4 - # in order to reprocess them using smaller time windows around theoretical onset - # get stations with bad onsets - badpicks = [] - for key in picks: - if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: - badpicks.append([key, picks[key]['P']['mpp']]) - - if len(badpicks) == 0: - print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") - # get NLLoc-location file - locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) - if len(glob.glob(locsearch)) > 0: - # get latest NLLoc-location file if several are available - nllocfile = max(glob.glob(locsearch), key=os.path.getctime) - # calculating seismic moment Mo and moment magnitude Mw - finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) - else: - print("autoPyLoT: No NLLoc-location file available!") - print("No source parameter estimation possible!") - else: - # get theoretical P-onset times from NLLoc-location file - locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) - if len(glob.glob(locsearch)) > 0: - # get latest file if several are available - nllocfile = max(glob.glob(locsearch), key=os.path.getctime) - nlloccounter = 0 - while len(badpicks) > 0 and nlloccounter <= maxnumit: - nlloccounter += 1 - if nlloccounter > maxnumit: - print("autoPyLoT: Number of maximum iterations reached, stop iterative picking!") - break - print("autoPyLoT: Starting with iteration No. %d ..." % nlloccounter) - picks = iteratepicker(wfdat, nllocfile, picks, badpicks, parameter) - # write phases to NLLoc-phase file - nll.export(picks, phasefile) - # remove actual NLLoc-location file to keep only the last - os.remove(nllocfile) - # locate the event - nll.locate(ctrfile) - print("autoPyLoT: Iteration No. %d finished." % nlloccounter) - # get updated NLLoc-location file - nllocfile = max(glob.glob(locsearch), key=os.path.getctime) - # check for bad picks - badpicks = [] - for key in picks: - if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: - badpicks.append([key, picks[key]['P']['mpp']]) - print("autoPyLoT: After iteration No. %d: %d bad onsets found ..." % (nlloccounter, \ - len(badpicks))) - if len(badpicks) == 0: - print("autoPyLoT: No more bad onsets found, stop iterative picking!") - nlloccounter = maxnumit - - # calculating seismic moment Mo and moment magnitude Mw - finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) - # get network moment magntiude - netMw = [] - for key in finalpicks.getpicdic(): - if finalpicks.getpicdic()[key]['P']['Mw'] is not None: - netMw.append(finalpicks.getpicdic()[key]['P']['Mw']) - netMw = np.median(netMw) - print("Network moment magnitude: %4.1f" % netMw) - else: - print("autoPyLoT: No NLLoc-location file available! Stop iteration!") - ########################################################## - # write phase files for various location routines - # HYPO71 - hypo71file = '%s/autoPyLoT_HYPO71.pha' % event - if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None: - hsat.export(finalpicks.getpicdic(), hypo71file) - data.applyEVTData(finalpicks.getpicdic()) - else: - hsat.export(picks, hypo71file) - data.applyEVTData(picks) - fnqml = '%s/autoPyLoT' % event - data.exportEvent(fnqml) - - endsplash = '''------------------------------------------\n' - -----Finished event %s!-----\n' - ------------------------------------------'''.format \ - (version=_getVersionString()) % evID - print(endsplash) - if locflag == 0: - print("autoPyLoT was running in non-location mode!") - - # single event processing + # multiple event processing + # read each event in database + events = [events for events in glob.glob(os.path.join(datapath, '*')) if os.path.isdir(events)] else: - data.setWFData(glob.glob(os.path.join(datapath, parameter.get('eventID'), '*'))) - print("Working on event {0}".format(parameter.get('eventID'))) + # single event processing + events = glob.glob(os.path.join(datapath, parameter.get('eventID'))) + for event in events: + data.setWFData(glob.glob(os.path.join(datapath, event, '*'))) + evID = os.path.split(event)[-1] + print('Working on event %s' % event) print(data) - wfdat = data.getWFData() # all available streams + metadata = read_metadata(parameter.get('invdir')) + corr_dat, rest_flag = restitute_data(wfdat.copy(), *metadata) ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) - finalpicks = picks ########################################################## # locating if locflag == 1: @@ -237,12 +125,14 @@ def autoPyLoT(inputfile): nll.export(picks, phasefile) # For locating the event the NLLoc-control file has to be modified! - nllocout = '%s_%s' % (parameter.get('eventID'), nllocoutpatter) + nllocout = '%s_%s' % (evID, nllocoutpatter) # create comment line for NLLoc-control file - nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, ttpat) + nll.modify_inputs(ctrf, nllocroot, nllocout, phasef, + ttpat) # locate the event nll.locate(ctrfile) + # !iterative picking if traces remained unpicked or occupied with bad picks! # get theoretical onset times for picks with weights >= 4 # in order to reprocess them using smaller time windows around theoretical onset @@ -252,18 +142,28 @@ def autoPyLoT(inputfile): if picks[key]['P']['weight'] >= 4 or picks[key]['S']['weight'] >= 4: badpicks.append([key, picks[key]['P']['mpp']]) + # TODO keep code DRY (Don't Repeat Yourself) the following part is written twice + # suggestion: delete block and modify the later similar block to work properly + if len(badpicks) == 0: print("autoPyLoT: No bad onsets found, thus no iterative picking necessary!") # get NLLoc-location file locsearch = '%s/loc/%s.????????.??????.grid?.loc.hyp' % (nllocroot, nllocout) if len(glob.glob(locsearch)) > 0: - # 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) + evt = read_events(nllocfile) # calculating seismic moment Mo and moment magnitude Mw - finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), + parameter.get('Qp'), parameter.get('rho'), True, 2) + local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2) + for station, amplitude in local_mag.amplitudes: + picks[station]['S']['Ao'] = amplitude + evt = add_amplitudes(evt, local_mag.amplitudes) + netML = local_mag.net_magnitude() + netMw = moment_mag.net_magnitude() + evt.origins[0].magnitudes.append(netMw) + evt.origins[0].magnitudes.append(netML) else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -300,38 +200,37 @@ def autoPyLoT(inputfile): if len(badpicks) == 0: print("autoPyLoT: No more bad onsets found, stop iterative picking!") nlloccounter = maxnumit - + evt = read_events(nllocfile) # calculating seismic moment Mo and moment magnitude Mw - finalpicks = MomentMagnitude(wfdat, None, None, parameter.get('iplot'), \ - nllocfile, picks, parameter.get('rho'), \ - parameter.get('vp'), parameter.get('Qp'), \ - parameter.get('invdir')) + moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), + parameter.get('Qp'), parameter.get('rho'), True, 2) + local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2) + for station, amplitude in local_mag.amplitudes: + picks[station]['S']['Ao'] = amplitude + evt = add_amplitudes(evt, local_mag.amplitudes) + netML = local_mag.net_magnitude() # get network moment magntiude - netMw = [] - for key in finalpicks.getpicdic(): - if finalpicks.getpicdic()[key]['P']['Mw'] is not None: - netMw.append(finalpicks.getpicdic()[key]['P']['Mw']) - netMw = np.median(netMw) - print("Network moment magnitude: %4.1f" % netMw) + netMw = moment_mag.net_magnitude() + evt.origins[0].magnitudes.append(netMw) + evt.origins[0].magnitudes.append(netML) + print("Network moment magnitude: %4.1f" % netMw.mag) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") ########################################################## # write phase files for various location routines # HYPO71 - hypo71file = '%s/%s/autoPyLoT_HYPO71.pha' % (datapath, parameter.get('eventID')) - if hasattr(finalpicks, 'getpicdic') and finalpicks.getpicdic() is not None: - hsat.export(finalpicks.getpicdic(), hypo71file) - data.applyEVTData(finalpicks.getpicdic()) - else: - hsat.export(picks, hypo71file) - data.applyEVTData(picks) - fnqml = '%s/%s/autoPyLoT' % (datapath, parameter.get('eventID')) + hypo71file = '%s/autoPyLoT_HYPO71.pha' % event + hsat.export(picks, hypo71file) + data.applyEVTData(picks) + if evt is not None: + data.applyEVTData(evt, 'event') + fnqml = '%s/autoPyLoT' % event data.exportEvent(fnqml) endsplash = '''------------------------------------------\n' -----Finished event %s!-----\n' ------------------------------------------'''.format \ - (version=_getVersionString()) % parameter.get('eventID') + (version=_getVersionString()) % evID print(endsplash) if locflag == 0: print("autoPyLoT was running in non-location mode!") From be2bacf5e8cd1b3ee51313e779a88457ed0e22ab Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 28 Sep 2016 14:37:24 +0200 Subject: [PATCH 24/32] bugfix: metadata not read from default file --- QtPyLoT.py | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index a2148fe3..e63be200 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -988,11 +988,11 @@ class MainWindow(QMainWindow): def calc_magnitude(self, type='ML'): - settings = QSettings() - fninv = settings.value("inventoryFile", None) - if fninv is None and not self.metadata: + def set_inv(settings): fninv, _ = QFileDialog.getOpenFileName(self, self.tr( "Select inventory..."), self.tr("Select file")) + if not fninv: + return False ans = QMessageBox.question(self, self.tr("Make default..."), self.tr( "New inventory filename set.\n" + \ @@ -1003,15 +1003,35 @@ class MainWindow(QMainWindow): settings.setValue("inventoryFile", fninv) settings.sync() self.metadata = read_metadata(fninv) + return True + + settings = QSettings() + fninv = settings.value("inventoryFile", None) + + if fninv is None and not self.metadata: + if not set_inv(settings): + return None + elif fninv is not None and not self.metadata: + ans = QMessageBox.question(self, self.tr("Use default..."), + self.tr( + "Do you want to use the default value?"), + QMessageBox.Yes | QMessageBox.No, + QMessageBox.Yes) + if ans == QMessageBox.No: + if not set_inv(settings): + return None + else: + self.metadata = read_metadata(fninv) + wf_copy = self.get_data().getWFData().copy() [corr_wf, rest_flag] = restitute_data(wf_copy, *self.metadata) if not rest_flag: raise ProcessingError('Restitution of waveform data failed!') if type == 'ML': - local_mag = RichterMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('sstop')) + local_mag = RichterMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('sstop'), verbosity = True) return local_mag.net_magnitude() elif type == 'Mw': - moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho')) + moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho'), verbosity = True) return moment_mag.net_magnitude() else: return None From d093349b50c3091aeb072e01bcb69371f79e033f Mon Sep 17 00:00:00 2001 From: Marcel Paffrath Date: Wed, 28 Sep 2016 14:38:16 +0200 Subject: [PATCH 25/32] changed function position --- pylot/RELEASE-VERSION | 2 +- pylot/core/active/seismicArrayPreparation.py | 70 ++++++++++---------- 2 files changed, 36 insertions(+), 36 deletions(-) diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION index 30ada269..4c41ac0c 100644 --- a/pylot/RELEASE-VERSION +++ b/pylot/RELEASE-VERSION @@ -1 +1 @@ -55a58-dirty +5f92-dirty diff --git a/pylot/core/active/seismicArrayPreparation.py b/pylot/core/active/seismicArrayPreparation.py index 976c37f3..907e9161 100644 --- a/pylot/core/active/seismicArrayPreparation.py +++ b/pylot/core/active/seismicArrayPreparation.py @@ -3,6 +3,41 @@ import sys import numpy as np from scipy.interpolate import griddata +def readMygridNlayers(filename): + infile = open(filename, 'r') + nlayers = len(infile.readlines()) / 2 + infile.close() + + return nlayers + +def readMygrid(filename): + ztop = []; + zbot = []; + vtop = []; + vbot = [] + infile = open(filename, 'r') + nlayers = readMygridNlayers(filename) + + print('\nreadMygrid: Reading file %s.' % filename) + for index in range(nlayers): + line1 = infile.readline() + line2 = infile.readline() + ztop.append(float(line1.split()[0])) + vtop.append(float(line1.split()[1])) + zbot.append(float(line2.split()[0])) + vbot.append(float(line2.split()[1])) + print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]' + '\n[Bot: v = %s [km/s], z = %s [m]]' + % (index + 1, vtop[index], ztop[index], + vbot[index], zbot[index])) + + if not ztop[0] == 0: + print('ERROR: there must be a velocity set for z = 0 in the file %s' % filename) + print('e.g.:\n0 0.33\n-5 1.0\netc.') + + infile.close() + return ztop, zbot, vtop, vbot + class SeisArray(object): ''' @@ -714,41 +749,6 @@ class SeisArray(object): # rad = angle / 180 * PI # return rad - def readMygridNlayers(filename): - infile = open(filename, 'r') - nlayers = len(infile.readlines()) / 2 - infile.close() - - return nlayers - - def readMygrid(filename): - ztop = []; - zbot = []; - vtop = []; - vbot = [] - infile = open(filename, 'r') - nlayers = readMygridNlayers(filename) - - print('\nreadMygrid: Reading file %s.' % filename) - for index in range(nlayers): - line1 = infile.readline() - line2 = infile.readline() - ztop.append(float(line1.split()[0])) - vtop.append(float(line1.split()[1])) - zbot.append(float(line2.split()[0])) - vbot.append(float(line2.split()[1])) - print('Layer %s:\n[Top: v = %s [km/s], z = %s [m]]' - '\n[Bot: v = %s [km/s], z = %s [m]]' - % (index + 1, vtop[index], ztop[index], - vbot[index], zbot[index])) - - if not ztop[0] == 0: - print('ERROR: there must be a velocity set for z = 0 in the file %s' % filename) - print('e.g.:\n0 0.33\n-5 1.0\netc.') - - infile.close() - return ztop, zbot, vtop, vbot - R = 6371. vmin = 0.34 decm = 0.3 # diagonal elements of the covariance matrix (grid3dg's default value is 0.3) From 019a3ae0f3fcc6e0c1aacabc123aae9752496473 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Wed, 28 Sep 2016 14:53:57 +0200 Subject: [PATCH 26/32] [new] added origin information to the net_magnitude --- pylot/core/analysis/magnitude.py | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 994c2230..4cf02535 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -82,6 +82,10 @@ class Magnitude(object): def event(self): return self._event + @property + def origin_id(self): + return self._event.origins[0].resource_id + @property def arrivals(self): return self._event.origins[0].arrivals @@ -107,7 +111,9 @@ class Magnitude(object): def net_magnitude(self): if self: - return ObspyMagnitude(mag=np.median([M["mag"] for M in self.magnitudes.values()]), type=self.type) + return ObspyMagnitude(mag=np.median([M["mag"] for M in + self.magnitudes.values()]), + type=self.type, origin_id=self.origin_id) return None From d68a1bcf0e4425c17a2d41cd734ef93661bd3852 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Wed, 28 Sep 2016 14:56:01 +0200 Subject: [PATCH 27/32] [fix] some bugs found and fixed --- autoPyLoT.py | 32 +++++++++++++++++++------------- 1 file changed, 19 insertions(+), 13 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 2bf14b76..0c727278 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -16,7 +16,8 @@ from pylot.core.io.data import Data from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.phases import add_amplitudes from pylot.core.pick.autopick import autopickevent, iteratepicker -from pylot.core.util.dataprocessing import restitute_data, read_metadata +from pylot.core.util.dataprocessing import restitute_data, read_metadata, \ + remove_underscores from pylot.core.util.structure import DATASTRUCTURE from pylot.core.util.version import get_git_version as _getVersionString @@ -113,6 +114,7 @@ def autoPyLoT(inputfile): print('Working on event %s' % event) print(data) wfdat = data.getWFData() # all available streams + wfdat = remove_underscores(wfdat) metadata = read_metadata(parameter.get('invdir')) corr_dat, rest_flag = restitute_data(wfdat.copy(), *metadata) ########################################################## @@ -152,18 +154,20 @@ def autoPyLoT(inputfile): if len(glob.glob(locsearch)) > 0: # get latest NLLoc-location file if several are available nllocfile = max(glob.glob(locsearch), key=os.path.getctime) - evt = read_events(nllocfile) + evt = read_events(nllocfile)[0] # calculating seismic moment Mo and moment magnitude Mw moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), - parameter.get('Qp'), parameter.get('rho'), True, 2) - local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2) - for station, amplitude in local_mag.amplitudes: + parameter.get('Qp'), + parameter.get('rho'), True, 0) + local_mag = RichterMagnitude(corr_dat, evt, + parameter.get('sstop'), True, 0) + for station, amplitude in local_mag.amplitudes.items(): picks[station]['S']['Ao'] = amplitude evt = add_amplitudes(evt, local_mag.amplitudes) netML = local_mag.net_magnitude() netMw = moment_mag.net_magnitude() - evt.origins[0].magnitudes.append(netMw) - evt.origins[0].magnitudes.append(netML) + mags = [netML, netMw] + evt.magnitudes += mags else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -200,19 +204,21 @@ def autoPyLoT(inputfile): if len(badpicks) == 0: print("autoPyLoT: No more bad onsets found, stop iterative picking!") nlloccounter = maxnumit - evt = read_events(nllocfile) + evt = read_events(nllocfile)[0] # calculating seismic moment Mo and moment magnitude Mw moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), - parameter.get('Qp'), parameter.get('rho'), True, 2) - local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 2) - for station, amplitude in local_mag.amplitudes: + parameter.get('Qp'), + parameter.get('rho'), True, 0) + local_mag = RichterMagnitude(corr_dat, evt, + parameter.get('sstop'), True, 0) + for station, amplitude in local_mag.amplitudes.items(): picks[station]['S']['Ao'] = amplitude evt = add_amplitudes(evt, local_mag.amplitudes) netML = local_mag.net_magnitude() # get network moment magntiude netMw = moment_mag.net_magnitude() - evt.origins[0].magnitudes.append(netMw) - evt.origins[0].magnitudes.append(netML) + mags = [netML, netMw] + evt.magnitudes += mags print("Network moment magnitude: %4.1f" % netMw.mag) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") From 010963dcd18841b97c91548d4e2470efdefdfc6c Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Wed, 28 Sep 2016 15:07:49 +0200 Subject: [PATCH 28/32] [bugfix] not all processing entries have to contain remove but at least one of them --- pylot/core/util/dataprocessing.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index a0ba4bcf..040c8ff3 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -215,7 +215,7 @@ def restitute_data(data, invtype, inobj, unit='VEL', force=False): seed_id = tr.get_id() # check, whether this trace has already been corrected if 'processing' in tr.stats.keys() \ - and np.all(['remove' in p for p in tr.stats.processing]) \ + and np.any(['remove' in p for p in tr.stats.processing]) \ and not force: print("Trace {0} has already been corrected!".format(seed_id)) continue From 900c7af9313e1300ef7dae72bfa082599fe6dcb2 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Thu, 29 Sep 2016 11:53:25 +0200 Subject: [PATCH 29/32] [new] added referenced information on Magnitude properties to the recently introduced Magnitude objects --- pylot/core/analysis/magnitude.py | 68 +++++++++++++++++++++++++++----- 1 file changed, 58 insertions(+), 10 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 4cf02535..a5477f2a 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -9,7 +9,7 @@ import os import matplotlib.pyplot as plt import numpy as np -from obspy.core.event import Magnitude as ObspyMagnitude +import obspy.core.event as ope from obspy.geodetics import degrees2kilometers from scipy import integrate, signal from scipy.optimize import curve_fit @@ -109,11 +109,22 @@ class Magnitude(object): def calc(self): pass + def updated_event(self): + self.event.magnitudes.append(self.net_magnitude()) + return self.event + def net_magnitude(self): if self: - return ObspyMagnitude(mag=np.median([M["mag"] for M in - self.magnitudes.values()]), - type=self.type, origin_id=self.origin_id) + # TODO if an average Magnitude instead of the median is calculated + # StationMagnitudeContributions should be added to the returned + # Magnitude object + # mag_error => weights (magnitude error estimate from peak_to_peak, calcsourcespec?) + # weights => StationMagnitdeContribution + mag = ope.Magnitude(mag=np.median([M.mag for M in self.magnitudes.values()]), + type=self.type, origin_id=self.origin_id, + station_count=len(self.magnitudes), + azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap) + return mag return None @@ -221,12 +232,25 @@ class RichterMagnitude(Magnitude): delta = degrees2kilometers(a.distance) onset = pick.time a0 = self.peak_to_peak(wf, onset) - self.amplitudes = (station, a0) + amplitude = ope.Amplitude(generic_amplitude=a0 * 1e-3) + amplitude.unit = 'm' + amplitude.category = 'point' + amplitude.waveform_id = pick.waveform_id + amplitude.magnitude_hint = self.type + amplitude.pick_id = pick.resource_id + amplitude.type = 'AML' + self.event.amplitudes.append(amplitude) + self.amplitudes = (station, amplitude) # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # use of custom relation functions - mag = dict(mag=np.log10(a0) + richter_magnitude_scaling(delta)) - self.magnitudes = (station, mag) + magnitude = ope.StationMagnitude(mag=np.log10(a0) + richter_magnitude_scaling(delta)) + magnitude.origin_id = self.origin_id + magnitude.waveform_id = pick.waveform_id + magnitude.amplitude_id = amplitude.resource_id + magnitude.station_magnitude_type = self.type + self.event.station_magnitudes.append(magnitude) + self.magnitudes = (station, magnitude) class MomentMagnitude(Magnitude): @@ -239,6 +263,8 @@ class MomentMagnitude(Magnitude): corresponding moment magntiude Mw. ''' + _props = dict() + def __init__(self, stream, event, vp, Qp, density, verbosity=False, iplot=False): super(MomentMagnitude, self).__init__(stream, event, verbosity, iplot) @@ -260,6 +286,23 @@ class MomentMagnitude(Magnitude): def rock_density(self): return self._density + @property + def moment_props(self): + return self._props + + @moment_props.setter + def moment_props(self, value): + station, props = value + self._props[station] = props + + @property + def seismic_moment(self): + return self._m0 + + @seismic_moment.setter + def seismic_moment(self, value): + self._m0 = value + def calc(self): for a in self.arrivals: if a.phase not in 'pP': @@ -281,9 +324,14 @@ class MomentMagnitude(Magnitude): print("WARNING: insufficient frequency information") continue wf = select_for_phase(wf, "P") - M0, Mw = calcMoMw(wf, w0, self.rock_density, self.p_velocity, distance, self.verbose) - mag = dict(w0=w0, fc=fc, M0=M0, mag=Mw) - self.magnitudes = (station, mag) + 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)) + magnitude = ope.StationMagnitude(mag=mw) + magnitude.origin_id = self.origin_id + magnitude.waveform_id = pick.waveform_id + magnitude.station_magnitude_type = self.type + self.event.station_magnitudes.append(magnitude) + self.magnitudes = (station, magnitude) def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False): From 34c1e80e78abefef65dc83c2bfcc14e4dea7d8f2 Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Thu, 29 Sep 2016 11:55:08 +0200 Subject: [PATCH 30/32] [enhancement] make use of the new methods introduced in [900c7af9313e1300ef7dae72bfa082599fe6dcb2] --- autoPyLoT.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index 0c727278..8c81b51d 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -14,7 +14,6 @@ import pylot.core.loc.nll as nll from pylot.core.analysis.magnitude import MomentMagnitude, RichterMagnitude from pylot.core.io.data import Data from pylot.core.io.inputs import AutoPickParameter -from pylot.core.io.phases import add_amplitudes from pylot.core.pick.autopick import autopickevent, iteratepicker from pylot.core.util.dataprocessing import restitute_data, read_metadata, \ remove_underscores @@ -159,15 +158,15 @@ def autoPyLoT(inputfile): moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), parameter.get('Qp'), parameter.get('rho'), True, 0) + # update pick with moment property values (w0, fc, Mo) + for station, props in moment_mag.moment_props.items(): + picks[station]['P'].update(props) + evt = moment_mag.updated_event() local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 0) for station, amplitude in local_mag.amplitudes.items(): - picks[station]['S']['Ao'] = amplitude - evt = add_amplitudes(evt, local_mag.amplitudes) - netML = local_mag.net_magnitude() - netMw = moment_mag.net_magnitude() - mags = [netML, netMw] - evt.magnitudes += mags + picks[station]['S']['Ao'] = amplitude.generic_amplitude + evt = local_mag.updated_event() else: print("autoPyLoT: No NLLoc-location file available!") print("No source parameter estimation possible!") @@ -209,17 +208,17 @@ def autoPyLoT(inputfile): moment_mag = MomentMagnitude(corr_dat, evt, parameter.get('vp'), parameter.get('Qp'), parameter.get('rho'), True, 0) + # update pick with moment property values (w0, fc, Mo) + for station, props in moment_mag.moment_props.items(): + picks[station]['P'].update(props) + evt = moment_mag.updated_event() local_mag = RichterMagnitude(corr_dat, evt, parameter.get('sstop'), True, 0) for station, amplitude in local_mag.amplitudes.items(): - picks[station]['S']['Ao'] = amplitude - evt = add_amplitudes(evt, local_mag.amplitudes) - netML = local_mag.net_magnitude() - # get network moment magntiude - netMw = moment_mag.net_magnitude() - mags = [netML, netMw] - evt.magnitudes += mags - print("Network moment magnitude: %4.1f" % netMw.mag) + picks[station]['S']['Ao'] = amplitude.generic_amplitude + evt = local_mag.updated_event() + net_mw = moment_mag.net_magnitude() + print("Network moment magnitude: %4.1f" % net_mw.mag) else: print("autoPyLoT: No NLLoc-location file available! Stop iteration!") ########################################################## From dfefd8af876261e542eba5375a1e04d0b7235761 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 29 Sep 2016 12:08:59 +0200 Subject: [PATCH 31/32] [enhancement] make use of new Magnitude method in QtPyLoT --- QtPyLoT.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index e63be200..87eadfd2 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -984,7 +984,7 @@ class MainWindow(QMainWindow): os.remove(phasepath) self.get_data().applyEVTData(lt.read_location(locpath), type='event') - self.get_data().get_evt_data().magnitudes.append(self.calc_magnitude()) + self.get_data().applyEVTData(self.calc_magnitude(), type='event') def calc_magnitude(self, type='ML'): @@ -1029,10 +1029,10 @@ class MainWindow(QMainWindow): raise ProcessingError('Restitution of waveform data failed!') if type == 'ML': local_mag = RichterMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('sstop'), verbosity = True) - return local_mag.net_magnitude() + return local_mag.updated_event() elif type == 'Mw': moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'), self.inputs.get('Qp'), self.inputs.get('rho'), verbosity = True) - return moment_mag.net_magnitude() + return moment_mag.updated_event() else: return None From 2e840cdfeb862ff372abdce744075f351f24b177 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 29 Sep 2016 12:44:37 +0200 Subject: [PATCH 32/32] [fix] reformatted code and fixed magnitude_type bug --- pylot/core/analysis/magnitude.py | 56 +++++++++++++++++++++----------- 1 file changed, 37 insertions(+), 19 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index a5477f2a..c75dacec 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -14,7 +14,8 @@ 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.pick.utils import getsignalwin, crossings_nonzero_all, \ + select_for_phase from pylot.core.util.utils import common_range, fit_curve @@ -40,7 +41,8 @@ class Magnitude(object): self._magnitudes = dict() 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') for s, m in self.magnitudes.items(): print('\t{0}\t{1}'.format(s, m)) @@ -120,10 +122,12 @@ class Magnitude(object): # Magnitude object # mag_error => weights (magnitude error estimate from peak_to_peak, calcsourcespec?) # weights => StationMagnitdeContribution - mag = ope.Magnitude(mag=np.median([M.mag for M in self.magnitudes.values()]), - type=self.type, origin_id=self.origin_id, - station_count=len(self.magnitudes), - azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap) + mag = ope.Magnitude( + mag=np.median([M.mag for M in self.magnitudes.values()]), + magnitude_type=self.type, + origin_id=self.origin_id, + station_count=len(self.magnitudes), + azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap) return mag return None @@ -206,8 +210,9 @@ class RichterMagnitude(Magnitude): plt.plot(th, sqH) plt.plot(th[iwin], sqH[iwin], 'g') plt.plot([t0, t0], [0, max(sqH)], 'r', linewidth=2) - plt.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \ - % (st[0].stats.station, wapp)) + plt.title( + 'Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \ + % (st[0].stats.station, wapp)) plt.xlabel('Time [s]') plt.ylabel('Displacement [mm]') plt.show() @@ -226,7 +231,8 @@ class RichterMagnitude(Magnitude): station=station), a.phase) 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)) continue delta = degrees2kilometers(a.distance) @@ -244,7 +250,8 @@ class RichterMagnitude(Magnitude): # using standard Gutenberg-Richter relation # TODO make the ML calculation more flexible by allowing # 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.waveform_id = pick.waveform_id magnitude.amplitude_id = amplitude.resource_id @@ -265,7 +272,8 @@ class MomentMagnitude(Magnitude): _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) self._vp = vp @@ -317,14 +325,17 @@ class MomentMagnitude(Magnitude): distance = degrees2kilometers(a.distance) azimuth = a.azimuth incidence = a.takeoff_angle - w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance, azimuth, - incidence, self.p_attenuation, self.plot_flag, self.verbose) + w0, fc = calcsourcespec(wf, onset, self.p_velocity, distance, + azimuth, + incidence, self.p_attenuation, + self.plot_flag, self.verbose) if w0 is None or fc is None: if self.verbose: print("WARNING: insufficient frequency information") continue 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)) magnitude = ope.StationMagnitude(mag=mw) 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] 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 - 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 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 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 @@ -525,7 +541,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence, 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)) + print("calcsourcespec: Using w0-value = %e m/Hz and fc = %f Hz" % ( + w0, fc)) if iplot > 1: f1 = plt.figure() @@ -644,7 +661,8 @@ def fitSourceModel(f, S, fc0, iplot, verbosity=False): fc = fc0 w0 = max(S) 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: plt.figure(iplot)