From bfa7ffc960de7484ce642d6d2c878be47fe3358c Mon Sep 17 00:00:00 2001 From: Sebastian Wehling-Benatelli Date: Fri, 23 Sep 2016 15:12:04 +0200 Subject: [PATCH] [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