From c8d8525c11f2683480ef40625704434b5ea1d539 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Thu, 15 Sep 2016 14:51:11 +0200 Subject: [PATCH 01/13] [bugs fixed and found] dataprocessing doesn't work as expected, np.bool_ substituted by bool --- pylot/core/analysis/magnitude.py | 2 +- pylot/core/util/dataprocessing.py | 13 ++++++++++--- 2 files changed, 11 insertions(+), 4 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index a22dbeba..41cfd590 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -353,7 +353,7 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp wf_copy = wfstream.copy() [cordat, restflag] = restitute_data(wf_copy, inventory) - if restflag == 1: + if restflag is True: zdat = cordat.select(component="Z") if len(zdat) == 0: zdat = cordat.select(component="3") diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index fb739fab..f59afbf1 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -155,7 +155,7 @@ def evt_head_check(root_dir, out_dir = None): print(nfiles) -def restitute_data(data, path_to_inventory, unit='VEL', force=False): +def restitute_data(data, path_to_inventory, unit='VEL', force=True): """ takes a data stream and a path_to_inventory and returns the corrected waveform data stream @@ -163,7 +163,7 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): :param path_to_inventory: path to inventory folder or file :param unit: unit to correct for (default: 'VEL') :param force: force restitution for already corrected traces (default: - False) + True) :return: corrected data stream """ @@ -193,6 +193,10 @@ def restitute_data(data, path_to_inventory, 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 not force: print("Trace {0} has already been corrected!".format(seed_id)) continue @@ -225,6 +229,9 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): else: finv = invlist[0] inventory = read_inventory(finv, format='STATIONXML') + else: + restflag.append(False) + continue # apply restitution to data if invtype in ['resp', 'dless']: tr.simulate(**kwargs) @@ -237,7 +244,7 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): # better try restitution for smaller subsets of data (e.g. station by # station) if len(restflag) > 0: - restflag = np.all(restflag) + restflag = bool(np.all(restflag)) else: restflag = False return data, restflag From 17585f93813bde0173cb6194c58cda7396d7b02b Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:29:33 +0200 Subject: [PATCH 02/13] [rename] renaming getGlobalTimes for consistency and introduction of similar new function in future commit --- QtPyLoT.py | 4 ++-- pylot/core/active/seismicshot.py | 4 ++-- pylot/core/io/data.py | 4 ++-- pylot/core/io/phases.py | 4 ++-- pylot/core/util/utils.py | 2 +- pylot/core/util/widgets.py | 6 +++--- 6 files changed, 12 insertions(+), 12 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 524643c6..00f4850a 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -55,7 +55,7 @@ from pylot.core.util.errors import FormatError, DatastructureError, \ OverwriteError from pylot.core.util.connection import checkurl from pylot.core.util.utils import fnConstructor, getLogin, \ - getGlobalTimes + full_range from pylot.core.io.location import create_creation_info, create_event from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ WaveformWidget, PropertiesDlg, HelpForm, createAction, PickDlg, \ @@ -643,7 +643,7 @@ class MainWindow(QMainWindow): ans = self.data.setWFData(self.getWFFnames()) else: ans = False - self._stime = getGlobalTimes(self.get_data().getWFData())[0] + self._stime = full_range(self.get_data().getWFData())[0] if ans: self.plotWaveformData() return ans diff --git a/pylot/core/active/seismicshot.py b/pylot/core/active/seismicshot.py index 9f270c76..e74ecaed 100644 --- a/pylot/core/active/seismicshot.py +++ b/pylot/core/active/seismicshot.py @@ -740,11 +740,11 @@ class SeismicShot(object): self._drawCFs(traceID, folm) def _drawStream(self, traceID, refresh=False, ax=None): - from pylot.core.util.utils import getGlobalTimes + from pylot.core.util.utils import full_range from pylot.core.util.utils import prepTimeAxis stream = self.getSingleStream(traceID) - stime = getGlobalTimes(stream)[0] + stime = full_range(stream)[0] timeaxis = prepTimeAxis(stime, stream[0]) timeaxis -= stime diff --git a/pylot/core/io/data.py b/pylot/core/io/data.py index 0b49230f..c3ca8b85 100644 --- a/pylot/core/io/data.py +++ b/pylot/core/io/data.py @@ -10,7 +10,7 @@ from obspy.core.event import Event from pylot.core.io.phases import readPILOTEvent, picks_from_picksdict, \ picksdict_from_pilot, merge_picks from pylot.core.util.errors import FormatError, OverwriteError -from pylot.core.util.utils import fnConstructor, getGlobalTimes +from pylot.core.util.utils import fnConstructor, full_range class Data(object): @@ -133,7 +133,7 @@ class Data(object): """ - self.cuttimes = getGlobalTimes(self.getWFData()) + self.cuttimes = full_range(self.getWFData()) def getEventFileName(self): """ diff --git a/pylot/core/io/phases.py b/pylot/core/io/phases.py index bd255b34..d3064f6f 100644 --- a/pylot/core/io/phases.py +++ b/pylot/core/io/phases.py @@ -12,7 +12,7 @@ from pylot.core.io.inputs import AutoPickParameter from pylot.core.io.location import create_arrival, create_event, \ create_magnitude, create_origin, create_pick from pylot.core.pick.utils import select_for_phase -from pylot.core.util.utils import getOwner, getGlobalTimes, four_digits +from pylot.core.util.utils import getOwner, full_range, four_digits def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs): @@ -325,7 +325,7 @@ def reassess_pilot_event(root_dir, db_dir, event_id, out_dir=None, fn_param=None msg = 'no waveform data found for station {station}'.format(station=station) warnings.warn(msg, RuntimeWarning) continue - stime, etime = getGlobalTimes(sel_st) + stime, etime = full_range(sel_st) rel_pick = mpp - stime epp, lpp, spe = earllatepicker(sel_st, default.get('nfac{0}'.format(phase)), diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 225b276e..db9a3888 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -205,7 +205,7 @@ def four_digits(year): return year -def getGlobalTimes(stream): +def full_range(stream): ''' takes a stream object and returns the latest end and the earliest start time of all contained trace objects diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 0139e457..4d3e3d16 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -32,7 +32,7 @@ from pylot.core.pick.utils import getSNR, earllatepicker, getnoisewin, \ from pylot.core.pick.compare import Comparison from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS, LOCTOOLS, \ COMPPOSITION_MAP -from pylot.core.util.utils import prepTimeAxis, getGlobalTimes, scaleWFData, \ +from pylot.core.util.utils import prepTimeAxis, full_range, scaleWFData, \ demeanTrace, isSorted, findComboBoxIndex, clims @@ -425,7 +425,7 @@ class WaveformWidget(FigureCanvas): noiselevel=None, scaleddata=False, mapping=True): self.getAxes().cla() self.clearPlotDict() - wfstart, wfend = getGlobalTimes(wfdata) + wfstart, wfend = full_range(wfdata) nmax = 0 for n, trace in enumerate(wfdata): channel = trace.stats.channel @@ -539,7 +539,7 @@ class PickDlg(QDialog): else: self.data = data - self.stime, self.etime = getGlobalTimes(self.getWFData()) + self.stime, self.etime = full_range(self.getWFData()) # initialize plotting widget self.multicompfig = WaveformWidget(self) From 8d37e9299c2f35537f4d25661563bf9061d93b4e Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:32:00 +0200 Subject: [PATCH 03/13] [new] added new function to find common time window within a stream --- pylot/core/analysis/magnitude.py | 20 ++------------------ pylot/core/util/utils.py | 18 ++++++++++++++++++ 2 files changed, 20 insertions(+), 18 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 41cfd590..8ec44a92 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -15,6 +15,7 @@ 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 +from pylot.core.util.utils import common_range class Magnitude(object): @@ -360,25 +361,8 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp cordat_copy = cordat.copy() # get equal time stamps and lengths of traces # necessary for rotation of traces - tr0start = cordat_copy[0].stats.starttime - tr0start = tr0start.timestamp - tr0end = cordat_copy[0].stats.endtime - tr0end = tr0end.timestamp - tr1start = cordat_copy[1].stats.starttime - tr1start = tr1start.timestamp - tr1end = cordat_copy[1].stats.endtime - tr1end = tr1end.timestamp - tr2start = cordat_copy[2].stats.starttime - tr2start = tr2start.timestamp - tr2end = cordat_copy[0].stats.endtime - tr2end = tr2end.timestamp - trstart = UTCDateTime(max([tr0start, tr1start, tr2start])) - trend = UTCDateTime(min([tr0end, tr1end, tr2end])) + trstart, trend = common_range(cordat_copy) cordat_copy.trim(trstart, trend) - minlen = min([len(cordat_copy[0]), len(cordat_copy[1]), len(cordat_copy[2])]) - cordat_copy[0].data = cordat_copy[0].data[0:minlen] - cordat_copy[1].data = cordat_copy[1].data[0:minlen] - cordat_copy[2].data = cordat_copy[2].data[0:minlen] try: # rotate into LQT (ray-coordindate-) system using Obspy's rotate # L: P-wave direction diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index db9a3888..8daeaefa 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -205,6 +205,24 @@ def four_digits(year): return year +def common_range(stream): + ''' + takes a stream object and returns the earliest end and the latest start + time of all contained trace objects + :param stream: seismological data stream + :type stream: `~obspy.core.stream.Stream` + :return: maximum start time and minimum end time + ''' + max_start = None + min_end = None + for trace in stream: + if max_start is None or trace.stats.starttime > max_start: + max_start = trace.stats.starttime + if min_end is None or trace.stats.endtime < min_end: + min_end = trace.stats.endtime + return max_start, min_end + + def full_range(stream): ''' takes a stream object and returns the latest end and the earliest start From 5155efc710a843807e74f290e767ad18fb57c2e1 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:33:08 +0200 Subject: [PATCH 04/13] [bugfix] do not try to give a full filepath for searching issue --- QtPyLoT.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index 00f4850a..67ae0ed6 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -520,7 +520,7 @@ class MainWindow(QMainWindow): def getSavePath(e): print('warning: {0}'.format(e)) - directory = os.path.join(self.getRoot(), self.getEventFileName()) + directory = os.path.realpath(self.getRoot()) file_filter = "QuakeML file (*.xml);;VELEST observation file " \ "format (*.cnv);;NonLinLoc observation file (*.obs)" title = 'Save event data ...' From 8ee515e79f32c9888e464f9032b5eb88dd721eeb Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:34:03 +0200 Subject: [PATCH 05/13] [bugfix] do not continue calculation without given data --- QtPyLoT.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/QtPyLoT.py b/QtPyLoT.py index 67ae0ed6..cc4416c2 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -991,6 +991,8 @@ class MainWindow(QMainWindow): 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, fninv, self.inputs.get('vp'), dist, From fa19ae9b9c58e35857bd7f80239e08208edb8531 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:35:59 +0200 Subject: [PATCH 06/13] [bugfix] only try to calculate moment magnitude given w0 and fc --- QtPyLoT.py | 11 +++++++++-- 1 file changed, 9 insertions(+), 2 deletions(-) diff --git a/QtPyLoT.py b/QtPyLoT.py index cc4416c2..1ff21462 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -998,9 +998,16 @@ class MainWindow(QMainWindow): w0, fc = calcsourcespec(wf, onset, fninv, self.inputs.get('vp'), dist, a.azimuth, a.takeoff_angle, self.inputs.get('Qp'), 0) - stat_mags = calcMoMw(wf, w0, 2700., 3000., dist, fninv) - mags[station] = stat_mags + if w0 is None or fc is None: + continue + station_mag = calcMoMw(wf, w0, 2700., 3000., dist, fninv) + 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 ce4ac4fd043d002fcaf3083c608fde3430c33feb Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:36:51 +0200 Subject: [PATCH 07/13] [pep8] use naming and style conventions --- pylot/core/analysis/magnitude.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 8ec44a92..a7a2befa 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -380,10 +380,10 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp # integrate to displacement # unrotated vertical component (for copmarison) - inttrz = signal.detrend(integrate.cumtrapz(zdat[0].data, None, \ + 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 = signal.detrend(integrate.cumtrapz(ldat[0].data, None, ldat[0].stats.delta)) # get window after P pulse for @@ -455,7 +455,8 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp # use of implicit scipy otimization function fit = synthsourcespec(F, w0in, Fcin) - [optspecfit, pcov] = curve_fit(synthsourcespec, F, YYcor, [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" From f34262d9319433f0083ae1d5baa46dfa4ffeca36 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Mon, 19 Sep 2016 11:39:15 +0200 Subject: [PATCH 08/13] [enhancement] catch possible exceptions during restitution process without losing code's verbosity --- pylot/core/util/dataprocessing.py | 33 ++++++++++++++++++++----------- 1 file changed, 22 insertions(+), 11 deletions(-) diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index f59afbf1..aff245c8 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -155,7 +155,7 @@ def evt_head_check(root_dir, out_dir = None): print(nfiles) -def restitute_data(data, path_to_inventory, unit='VEL', force=True): +def restitute_data(data, path_to_inventory, unit='VEL', force=False): """ takes a data stream and a path_to_inventory and returns the corrected waveform data stream @@ -188,6 +188,9 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=True): print("Neither dataless-SEED file,inventory-xml file nor RESP-file " "found!") return data, False + elif invtype == 'dless': # prevent multiple read of large dlsv + if len(inv[invtype]) == 1: + fname = Parser(inv[invtype][0]) # loop over traces for tr in data: @@ -197,11 +200,12 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=True): # 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 not force: + if 'processing' in tr.stats.keys() \ + and np.all(['remove' in p for p in tr.stats.processing]) \ + and not force: print("Trace {0} has already been corrected!".format(seed_id)) continue stime = tr.stats.starttime - seedresp = None prefilt = get_prefilt(tr) if invtype == 'resp': fresp = find_in_list(inv[invtype], seed_id) @@ -216,8 +220,6 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=True): elif invtype == 'dless': if len(inv[invtype]) > 1: fname = Parser(find_in_list(inv[invtype], seed_id)) - else: - fname = Parser(inv[invtype][0]) seedresp = dict(filename=fname, date=stime, units=unit) @@ -233,12 +235,21 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=True): restflag.append(False) continue # apply restitution to data - if invtype in ['resp', 'dless']: - tr.simulate(**kwargs) - else: - tr.attach_response(inventory) - tr.remove_response(output=unit, - pre_filt=prefilt) + try: + if invtype in ['resp', 'dless']: + tr.simulate(**kwargs) + else: + tr.attach_response(inventory) + tr.remove_response(output=unit, + pre_filt=prefilt) + except ValueError as e: + msg0 = 'Response for {0} not found in Parser'.format(seed_id) + msg1 = 'evalresp failed to calculate response' + if msg0 not in e.message or msg1 not in e.message: + raise + else: + restflag.append(False) + continue restflag.append(True) # check if ALL traces could be restituted, take care of large datasets # better try restitution for smaller subsets of data (e.g. station by From c73435dec3c655577b892fac6446df50843d4428 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Tue, 20 Sep 2016 09:54:14 +0200 Subject: [PATCH 09/13] [fix] do not calculate moment magnitude for S phases --- QtPyLoT.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/QtPyLoT.py b/QtPyLoT.py index 1ff21462..d946de14 100755 --- a/QtPyLoT.py +++ b/QtPyLoT.py @@ -988,6 +988,8 @@ class MainWindow(QMainWindow): settings.setValue("inventoryFile", fninv) settings.sync() 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) From 4a6b653a7257b32189c989f591e29db05046f298 Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Tue, 20 Sep 2016 09:55:54 +0200 Subject: [PATCH 10/13] [new] added new function to read metadata from disk this new function prevents multiple reading of large dataless seed volume to enhance overall performance --- pylot/core/util/dataprocessing.py | 65 ++++++++++++++++++++----------- 1 file changed, 42 insertions(+), 23 deletions(-) diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index aff245c8..497f4762 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -155,22 +155,16 @@ def evt_head_check(root_dir, out_dir = None): print(nfiles) -def restitute_data(data, path_to_inventory, unit='VEL', force=False): +def read_metadata(path_to_inventory): """ - takes a data stream and a path_to_inventory and returns the corrected - waveform data stream - :param data: seismic data stream - :param path_to_inventory: path to inventory folder or file - :param unit: unit to correct for (default: 'VEL') - :param force: force restitution for already corrected traces (default: - True) - :return: corrected data stream + take path_to_inventory and return either the corresponding list of files + found or the Parser object for a network dataless seed volume to prevent + read overhead for large dataless seed volumes + :param path_to_inventory: + :return: tuple containing a either list of files or `obspy.io.xseed.Parser` + object and the inventory type found + :rtype: tuple """ - - restflag = list() - - data = remove_underscores(data) - dlfile = list() invfile = list() respfile = list() @@ -185,12 +179,35 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): invtype = key_for_set_value(inv) if invtype is None: - print("Neither dataless-SEED file,inventory-xml file nor RESP-file " - "found!") - return data, False - elif invtype == 'dless': # prevent multiple read of large dlsv + raise IOError("Neither dataless-SEED file, inventory-xml file nor " + "RESP-file found!") + elif invtype == 'dless': # prevent multiple read of large dlsv if len(inv[invtype]) == 1: - fname = Parser(inv[invtype][0]) + robj = Parser(inv[invtype][0]) + else: + robj = inv[invtype] + else: + robj = inv[invtype] + return invtype, robj + + +def restitute_data(data, invtype, inobj, unit='VEL', force=False): + """ + takes a data stream and a path_to_inventory and returns the corrected + waveform data stream + :param data: seismic data stream + :param invtype: type of found metadata + :param inobj: either list of metadata files or `obspy.io.xseed.Parser` + object + :param unit: unit to correct for (default: 'VEL') + :param force: force restitution for already corrected traces (default: + False) + :return: corrected data stream + """ + + restflag = list() + + data = remove_underscores(data) # loop over traces for tr in data: @@ -208,7 +225,7 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): stime = tr.stats.starttime prefilt = get_prefilt(tr) if invtype == 'resp': - fresp = find_in_list(inv[invtype], seed_id) + fresp = find_in_list(inobj, seed_id) if not fresp: raise IOError('no response file found ' 'for trace {0}'.format(seed_id)) @@ -218,14 +235,16 @@ def restitute_data(data, path_to_inventory, unit='VEL', force=False): units=unit) kwargs = dict(paz_remove=None, pre_filt=prefilt, seedresp=seedresp) elif invtype == 'dless': - if len(inv[invtype]) > 1: - fname = Parser(find_in_list(inv[invtype], seed_id)) + if type(inobj) is list: + fname = Parser(find_in_list(inobj, seed_id)) + else: + fname = inobj seedresp = dict(filename=fname, date=stime, units=unit) kwargs = dict(pre_filt=prefilt, seedresp=seedresp) elif invtype == 'xml': - invlist = inv[invtype] + invlist = inobj if len(invlist) > 1: finv = find_in_list(invlist, seed_id) else: From 7e76bf75773912f6838023e4db57bf33bc8dadeb Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Tue, 20 Sep 2016 09:57:09 +0200 Subject: [PATCH 11/13] [change] make use of new metadata reading utility function to improve performance --- pylot/core/pick/autopick.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/pylot/core/pick/autopick.py b/pylot/core/pick/autopick.py index da496bdf..3b42eb1a 100755 --- a/pylot/core/pick/autopick.py +++ b/pylot/core/pick/autopick.py @@ -17,7 +17,7 @@ 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 +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 @@ -122,6 +122,7 @@ def autopickstation(wfstream, pickparam, verbose=False): 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,7 +566,7 @@ def autopickstation(wfstream, pickparam, verbose=False): hdat = edat.copy() hdat += ndat h_copy = hdat.copy() - [cordat, restflag] = restitute_data(h_copy, invdir) + [cordat, restflag] = restitute_data(h_copy, invtype, inventory) # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude if restflag: @@ -602,7 +603,7 @@ def autopickstation(wfstream, pickparam, verbose=False): hdat = edat.copy() hdat += ndat h_copy = hdat.copy() - [cordat, restflag] = restitute_data(h_copy, invdir) + [cordat, restflag] = restitute_data(h_copy, invtype, inventory) if restflag == 1: # calculate WA-peak-to-peak amplitude # using subclass WApp of superclass Magnitude From df002ce9acbd1543c21c920b0a22706cd0926c1e Mon Sep 17 00:00:00 2001 From: Sebastianw Wehling-Benatelli Date: Tue, 20 Sep 2016 09:58:33 +0200 Subject: [PATCH 12/13] [change] use read in metadata information instead of reading metadata each time invoked --- pylot/core/analysis/magnitude.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index a7a2befa..6c2ab843 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -306,7 +306,8 @@ def calcMoMw(wfstream, w0, rho, vp, delta, inv): return Mo, Mw -def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp, iplot): +def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, + qp, iplot): ''' 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 @@ -353,7 +354,9 @@ def calcsourcespec(wfstream, onset, inventory, vp, delta, azimuth, incidence, qp data = Data() wf_copy = wfstream.copy() - [cordat, restflag] = restitute_data(wf_copy, inventory) + invtype, inventory = metadata + + [cordat, restflag] = restitute_data(wf_copy, invtype, inventory) if restflag is True: zdat = cordat.select(component="Z") if len(zdat) == 0: From 6a2bbe3f91d9ff77ce4a4188fed886e1da484807 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Ludger=20K=C3=BCperkoch?= Date: Tue, 20 Sep 2016 10:40:21 +0200 Subject: [PATCH 13/13] Stabilized zero-crosings determination for source spectrum estimation from P pulse. --- pylot/core/analysis/magnitude.py | 14 ++------------ 1 file changed, 2 insertions(+), 12 deletions(-) diff --git a/pylot/core/analysis/magnitude.py b/pylot/core/analysis/magnitude.py index 6c2ab843..ad54b138 100644 --- a/pylot/core/analysis/magnitude.py +++ b/pylot/core/analysis/magnitude.py @@ -390,21 +390,11 @@ def calcsourcespec(wfstream, onset, metadata, vp, delta, azimuth, incidence, ldat[0].stats.delta)) # get window after P pulse for - # calculating source spectrum - if zdat[0].stats.sampling_rate <= 100: - winzc = zdat[0].stats.sampling_rate - elif zdat[0].stats.sampling_rate > 100 and \ - zdat[0].stats.sampling_rate <= 200: - winzc = 0.5 * zdat[0].stats.sampling_rate - elif zdat[0].stats.sampling_rate > 200 and \ - zdat[0].stats.sampling_rate <= 400: - winzc = 0.2 * zdat[0].stats.sampling_rate - elif zdat[0].stats.sampling_rate > 400: - winzc = zdat[0].stats.sampling_rate + # calculating source spectrum tstart = UTCDateTime(zdat[0].stats.starttime) tonset = onset.timestamp - tstart.timestamp impickP = tonset * zdat[0].stats.sampling_rate - wfzc = Ldat[impickP: impickP + winzc] + wfzc = Ldat[impickP: len(Ldat) - 1] # get time array t = np.arange(0, len(inttrz) * zdat[0].stats.delta, \ zdat[0].stats.delta)