From ef8ebc300e1fd6261405d5b4d1617070c70336cd Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 28 Aug 2015 11:27:09 +0200 Subject: [PATCH 1/7] New object to calculate magnitude. Finished class wapp to calculate amplitude as seen on Wood-Anderson seismograph. --- pylot/core/analysis/magnitude.py | 124 +++++++++++++++++++++++++++++++ 1 file changed, 124 insertions(+) create mode 100644 pylot/core/analysis/magnitude.py diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py new file mode 100644 index 00000000..c9db0ef2 --- /dev/null +++ b/pylot/core/analysis/magnitude.py @@ -0,0 +1,124 @@ +# -*- coding: utf-8 -*- +""" +Created August/September 2015. + +:author: Ludger Küperkoch / MAGS2 EP3 working group +""" + +import matplotlib.pyplot as plt +import numpy as np +from obspy.core import Stream +from pylot.core.pick.utils import getsignalwin + +class Magnitude(object): + ''' + Superclass for calculating Wood-Anderson peak-to-peak + amplitudes, local magnitudes and moment magnitudes. + ''' + + def __init__(self, wfstream, To, pwin, iplot): + ''' + :param: wfstream + :type: `~obspy.core.stream.Stream + + :param: To, onset time, P- or S phase + :type: float + + :param: pwin, pick window [To To+pwin] to get maximum + peak-to-peak amplitude + :type: float + + :param: iplot, no. of figure window for plotting interims results + :type: integer + + ''' + + assert isinstance(wfstream, Stream), "%s is not a stream object" % str(wfstream) + + self.setwfstream(wfstream) + self.setTo(To) + self.setpwin(pwin) + self.setiplot(iplot) + self.calcwapp() + + + def getwfstream(self): + return self.wfstream + + def setwfstream(self, wfstream): + self.wfstream = wfstream + + def getTo(self): + return self.To + + def setTo(self, To): + self.To = To + + def getpwin(self): + return self.pwin + + def setpwin(self, pwin): + self.pwin = pwin + + def getiplot(self): + return self.iplot + + def setiplot(self, iplot): + self.iplot = iplot + + def getwapp(self): + return self.wapp + + def calcwapp(self): + self.wapp = None + +class WApp(Magnitude): + ''' + Method to derive peak-to-peak amplitude as seen on a Wood-Anderson- + seismograph. Has to be derived from corrected traces! + ''' + + def calcwapp(self): + print "Getting Wood-Anderson peak-to-peak amplitude ..." + print "Simulating Wood-Anderson seismograph ..." + + self.wapp = None + stream = self.getwfstream() + + # 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 + 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) + # get maximum peak within pick window + iwin = getsignalwin(th, self.getTo(), self.getpwin()) + self.wapp = np.max(sqH[iwin]) + print "Determined Wood-Anderson peak-to-peak amplitude: %f mm" % self.wapp + if self.getiplot() > 1: + stream.plot() + 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.title('Station %s, RMS Horizontal Traces, WA-peak-to-peak=%4.1f mm' \ + % (stream[0].stats.station, self.wapp)) + plt.xlabel('Time [s]') + plt.ylabel('Displacement [mm]') + plt.show() + raw_input() + plt.close(f) + + + From 533ccc7b5c383d4fd6139f2bc9b12e38d532444d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 28 Aug 2015 11:29:00 +0200 Subject: [PATCH 2/7] Implemented new class wapp to calculate Wood-Anderson amplitudes for local magnitude calculation. Before calculating Wood-Anderson amplitude the certain traces are instrument corrected. --- pylot/core/pick/autopick.py | 36 ++++++++++++++++++++++++++++++++++++ 1 file changed, 36 insertions(+) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 3c074a7d..7113d90a 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -15,6 +15,8 @@ from pylot.core.pick.Picker import AICPicker, PragPicker 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.read.data import Data +from pylot.core.analysis.magnitude import WApp def autopickevent(data, param): stations = [] @@ -109,6 +111,8 @@ def autopickstation(wfstream, pickparam): nfacsl = pickparam.getParam('noisefactor') # parameter to check for spuriously picked S onset zfac = pickparam.getParam('zfac') + # path to inventory-, dataless- or resp-files + invdir = pickparam.getParam('invdir') # initialize output Pweight = 4 # weight for P onset @@ -132,6 +136,7 @@ def autopickstation(wfstream, pickparam): Pflag = 0 Sflag = 0 Pmarker = [] + Ao = None # split components zdat = wfstream.select(component="Z") @@ -501,6 +506,20 @@ def autopickstation(wfstream, pickparam): print 'autopickstation: S-weight: %d, SNR: %f, SNR[dB]: %f' % ( Sweight, SNRS, SNRSdB) + ################################################################## + # get Wood-Anderson peak-to-peak amplitude + print "################################################" + # initialize Data object + data = Data() + # re-create stream object including both horizontal components + hdat = edat.copy() + hdat += ndat + h_copy = hdat.copy() + cordat = data.restituteWFData(invdir, h_copy) + # calculate WA-peak-to-peak amplitude + # using subclass WApp of superclass Magnitude + wapp = WApp(cordat, mpickS, 10, iplot) + Ao = wapp.getwapp() else: print 'Bad initial (AIC) S-pick, skipping this onset!' @@ -509,6 +528,21 @@ def autopickstation(wfstream, pickparam): print '(min. AIC-SNR=', minAICSSNR, ', min. AIC-Slope=', \ minAICSslope, 'counts/s)' + ############################################################ + # get Wood-Anderson peak-to-peak amplitude + print "################################################" + # initialize Data object + data = Data() + # re-create stream object including both horizontal components + hdat = edat.copy() + hdat += ndat + h_copy = hdat.copy() + cordat = data.restituteWFData(invdir, h_copy) + # calculate WA-peak-to-peak amplitude + # using subclass WApp of superclass Magnitude + wapp = WApp(cordat, mpickP, 20, iplot) + Ao = wapp.getwapp() + else: print 'autopickstation: No horizontal component data available or ' \ 'bad P onset, skipping S picking!' @@ -693,5 +727,7 @@ def autopickstation(wfstream, pickparam): phasepick = {'lpp': lpickS, 'epp': epickS, 'mpp': mpickS, 'spe': Serror, \ 'snr': SNRS, 'snrdb': SNRSdB, 'weight': Sweight, 'fm': None} picks[phase] = phasepick + # add Wood-Anderson amplitude + picks[phase]['Ao'] = Ao return picks From 827a07a7ef652c9ba8417a3d89460070625c8b66 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 28 Aug 2015 11:31:19 +0200 Subject: [PATCH 3/7] Skipped restitution of entire traces at the beginning of processing as this is too time consuming. Instead, only traces are corrected within autopickstation, where at least automatically a P pick has been set. --- autoPyLoT.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/autoPyLoT.py b/autoPyLoT.py index c4e941ea..3084f1ad 100755 --- a/autoPyLoT.py +++ b/autoPyLoT.py @@ -70,8 +70,6 @@ def autoPyLoT(inputfile): # get path to inventory or dataless-seed file with station meta data invdir = parameter.getParam('invdir') - # get corner frequencies for pre-filtering traces - prefilt = parameter.getParam('prefilt') # multiple event processing # read each event in database @@ -83,8 +81,6 @@ def autoPyLoT(inputfile): print data wfdat = data.getWFData() # all available streams - # restitute waveform data getting responses from inventory-file - wfdat = data.restituteWFData(invdir, prefilt) ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) @@ -100,8 +96,6 @@ def autoPyLoT(inputfile): print data wfdat = data.getWFData() # all available streams - # restitute waveform data getting responses from inventory-file - wfdat = data.restituteWFData(invdir, prefilt) ########################################################## # !automated picking starts here! picks = autopickevent(wfdat, parameter) From eb592a3426d970da0c5d6524a43a58ecc77d732f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Fri, 28 Aug 2015 11:39:39 +0200 Subject: [PATCH 4/7] Claculation of Wood-Anderson amplitude only, if S-weight < 4. --- pylot/core/pick/autopick.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 7113d90a..9c84d38f 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -507,19 +507,20 @@ def autopickstation(wfstream, pickparam): print 'autopickstation: S-weight: %d, SNR: %f, SNR[dB]: %f' % ( Sweight, SNRS, SNRSdB) ################################################################## - # get Wood-Anderson peak-to-peak amplitude - print "################################################" - # initialize Data object - data = Data() - # re-create stream object including both horizontal components - hdat = edat.copy() - hdat += ndat - h_copy = hdat.copy() - cordat = data.restituteWFData(invdir, h_copy) - # calculate WA-peak-to-peak amplitude - # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickS, 10, iplot) - Ao = wapp.getwapp() + if Sweight < 4: + # get Wood-Anderson peak-to-peak amplitude + print "################################################" + # initialize Data object + data = Data() + # re-create stream object including both horizontal components + hdat = edat.copy() + hdat += ndat + h_copy = hdat.copy() + cordat = data.restituteWFData(invdir, h_copy) + # calculate WA-peak-to-peak amplitude + # using subclass WApp of superclass Magnitude + wapp = WApp(cordat, mpickS, 10, iplot) + Ao = wapp.getwapp() else: print 'Bad initial (AIC) S-pick, skipping this onset!' From 1f7049691c5ed78145836165668bf2d5dc6c6f5f Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Fri, 28 Aug 2015 16:01:42 +0200 Subject: [PATCH 5/7] [addresses #167] started fixing the multiple phase saving issue --- QtPyLoT.py | 13 +++++++++---- pylot/RELEASE-VERSION | 2 +- pylot/core/read/data.py | 13 +++++++++++-- 3 files changed, 21 insertions(+), 7 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 254e88f7..d14c51a0 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -405,13 +405,16 @@ class MainWindow(QMainWindow): except OverwriteError: msgBox = QMessageBox() msgBox.setText("Picks have been modified!") - msgBox.setInformativeText("Do you want to overwrite the picks and save?") - msgBox.setStandardButtons(QMessageBox.Save | QMessageBox.Discard | - QMessageBox.Cancel) + msgBox.setInformativeText("Do you want to save the changes and overwrite the picks?") + msgBox.setDetailedText(self.getData().getPicksStr()) + msgBox.setStandardButtons(QMessageBox.Save | QMessageBox.Cancel) msgBox.setDefaultButton(QMessageBox.Save) ret = msgBox.exec_() if ret == QMessageBox.Save: - print('Overwrite and Save') + self.getData().resetPicks() + self.saveData() + elif ret == QMessageBox.Cancel: + return False try: self.getData().exportEvent(self.fname, exform) except FormatError: @@ -425,6 +428,8 @@ class MainWindow(QMainWindow): fbasename, exform = os.path.splitext(fname[0]) if not fbasename: return False + elif not exform: + exform = fname[1].split('*')[1][:-1] self.getData().exportEvent(fbasename, exform) return True diff --git a/pylot/RELEASE-VERSION b/pylot/RELEASE-VERSION index 52c00f77..541bd94d 100644 --- a/pylot/RELEASE-VERSION +++ b/pylot/RELEASE-VERSION @@ -1 +1 @@ -1abc-dirty +497c-dirty diff --git a/pylot/core/read/data.py b/pylot/core/read/data.py index f3e07d30..3c83e9aa 100644 --- a/pylot/core/read/data.py +++ b/pylot/core/read/data.py @@ -52,6 +52,13 @@ class Data(object): def __str__(self): return str(self.wfdata) + def getPicksStr(self): + picks_str = '' + for pick in self.getEvtData().picks: + picks_str += str(pick) + '\n' + return picks_str + + def getParent(self): """ @@ -360,9 +367,11 @@ class Data(object): def applyPicks(picks): """ - + Creates ObsPy pick objects and append it to the picks list from the + PyLoT dictionary contain all picks. :param picks: - :raise OverwriteError: + :raise OverwriteError: raises an OverwriteError if the picks list is + not empty. The GUI will then ask for a decision. """ firstonset = None if self.getEvtData().picks: From d756f5d2e1c45b717f899ae2434bc533a4b0bb6f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Mon, 31 Aug 2015 10:10:42 +0200 Subject: [PATCH 6/7] Replaced hard coded window length for getting Wood-Anderson peak-to-peak amplitude with formerly set window length for calculating CF for S-picking. --- pylot/core/pick/autopick.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index 9c84d38f..b8532638 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -519,7 +519,7 @@ def autopickstation(wfstream, pickparam): cordat = data.restituteWFData(invdir, h_copy) # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickS, 10, iplot) + wapp = WApp(cordat, mpickS, mpickP + sstop, iplot) Ao = wapp.getwapp() else: @@ -541,7 +541,7 @@ def autopickstation(wfstream, pickparam): cordat = data.restituteWFData(invdir, h_copy) # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude - wapp = WApp(cordat, mpickP, 20, iplot) + wapp = WApp(cordat, mpickP, mpickP + sstop, iplot) Ao = wapp.getwapp() else: From fb3b599f50dcd7660fc03c1e56f7e0e8f97d63f4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Mon, 31 Aug 2015 10:24:17 +0200 Subject: [PATCH 7/7] restituteWFData: If input streams is None, a copy of streams derived by self.getWFData() is used for further processing. --- pylot/core/read/data.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/pylot/core/read/data.py b/pylot/core/read/data.py index 3c83e9aa..7e12f793 100644 --- a/pylot/core/read/data.py +++ b/pylot/core/read/data.py @@ -229,7 +229,8 @@ class Data(object): :return: """ if streams is None: - st = self.getWFData() + st_raw = self.getWFData() + st = st_raw.copy() else: st = streams