diff --git a/PyLoT.py b/PyLoT.py index cd166817..51b13f70 100755 --- a/PyLoT.py +++ b/PyLoT.py @@ -46,6 +46,9 @@ from obspy import UTCDateTime from obspy.core.event import Magnitude, Origin from obspy.core.util import AttribDict +from pylot.core.util.obspyDMT_interface import check_obspydmt_structure + + try: import pyqtgraph as pg except Exception as e: @@ -75,7 +78,8 @@ from pylot.core.util.dataprocessing import read_metadata, restitute_data from pylot.core.util.utils import fnConstructor, getLogin, \ full_range, readFilterInformation, trim_station_components, check4gaps, make_pen, pick_color_plt, \ pick_linestyle_plt, remove_underscores, check4doubled, identifyPhaseID, excludeQualityClasses, has_spe, \ - check4rotated, transform_colors_mpl, transform_colors_mpl_str, getAutoFilteroptions + check4rotated, transform_colors_mpl, transform_colors_mpl_str, getAutoFilteroptions, check_all_obspy, \ + check_all_pylot, real_Bool from pylot.core.util.event import Event from pylot.core.io.location import create_creation_info, create_event from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \ @@ -143,6 +147,7 @@ class MainWindow(QMainWindow): # default factor for dataplot e.g. enabling/disabling scrollarea self.height_factor = 12 + self.plot_method = 'normal' # UI has to be set up before(!) children widgets are about to show up self.createAction = createAction @@ -987,11 +992,13 @@ class MainWindow(QMainWindow): ''' Return waveform filenames from event in eventbox. ''' + # TODO: add dataStructure class for obspyDMT here, this is just a workaround! + eventpath = self.get_current_event_path(eventbox) + basepath = eventpath.split(os.path.basename(eventpath))[0] if self.dataStructure: - directory = self.get_current_event_path(eventbox) - if not directory: + if not eventpath: return - fnames = [os.path.join(directory, f) for f in os.listdir(directory)] + fnames = [os.path.join(eventpath, f) for f in os.listdir(eventpath)] else: raise DatastructureError('not specified') return fnames @@ -1036,13 +1043,15 @@ class MainWindow(QMainWindow): ed = getExistingDirectories(self, 'Select event directories...') if ed.exec_(): eventlist = ed.selectedFiles() - # select only folders that start with 'e', containin two dots and have length 12 - eventlist = [item for item in eventlist if item.split('/')[-1].startswith('e') - and len(item.split('/')[-1].split('.')) == 3 - and len(item.split('/')[-1]) == 12] + basepath = eventlist[0].split(os.path.basename(eventlist[0]))[0] + if check_obspydmt_structure(basepath): + print('Recognized obspyDMT structure in selected files.') + eventlist = check_all_obspy(eventlist) + else: + eventlist = check_all_pylot(eventlist) if not eventlist: print('No events found! Expected structure for event folders: [eEVID.DOY.YR],\n' - ' e.g. eventID=1, doy=2, yr=2016: e0001.002.16') + ' e.g. eventID=1, doy=2, yr=2016: e0001.002.16 or obspyDMT database') return else: return @@ -1649,9 +1658,13 @@ class MainWindow(QMainWindow): # else: # ans = False self.fnames = self.getWFFnames_from_eventbox() + eventpath = self.get_current_event_path() + basepath = eventpath.split(os.path.basename(eventpath))[0] + obspy_dmt = check_obspydmt_structure(basepath) self.data.setWFData(self.fnames, checkRotated=True, - metadata=self.metadata) + metadata=self.metadata, + obspy_dmt=obspy_dmt) def check_plot_quantity(self): """ @@ -1659,29 +1672,34 @@ class MainWindow(QMainWindow): :rtype: None """ settings = QSettings() - nth_sample = int(settings.value("nth_sample") if settings.value("nth_sample") else 1) - npts_max = 1e6 + nth_sample = int(settings.value("nth_sample")) if settings.value("nth_sample") else 1 + npts_max = 1e7 npts = self.get_npts_to_plot() npts2plot = npts/nth_sample if npts2plot < npts_max: - return - nth_sample_new = int(np.ceil(npts/npts_max)) - message = "You are about to plot a huge dataset with {npts} datapoints. With a current setting of " \ - "nth_sample = {nth_sample} a total of {npts2plot} points will be plotted which is more " \ - "than the maximum setting of {npts_max}. " \ - "PyLoT recommends to raise nth_sample from {nth_sample} to {nth_sample_new}. Continue?" - - ans = QMessageBox.question(self, self.tr("Optimize plot performance..."), - self.tr(message.format(npts=npts, - nth_sample=nth_sample, - npts_max=npts_max, - nth_sample_new=nth_sample_new, - npts2plot=npts2plot)), - QMessageBox.Yes | QMessageBox.No, - QMessageBox.Yes) - if ans == QMessageBox.Yes: - settings.setValue("nth_sample", nth_sample_new) - settings.sync() + settings.setValue('large_dataset', False) + else: + settings.setValue('large_dataset', True) + self.update_status('Dataset is very large. Using fast plotting method (MIN/MAX)', 10000) + settings.sync() + # nth_sample_new = int(np.ceil(npts/npts_max)) + # message = "You are about to plot a huge dataset with {npts} datapoints. With a current setting of " \ + # "nth_sample = {nth_sample} a total of {npts2plot} points will be plotted which is more " \ + # "than the maximum setting of {npts_max}. " \ + # "PyLoT recommends to raise nth_sample from {nth_sample} to {nth_sample_new}. Do you want "\ + # "to change nth_sample to {nth_sample_new} now?" + # + # ans = QMessageBox.question(self, self.tr("Optimize plot performance..."), + # self.tr(message.format(npts=npts, + # nth_sample=nth_sample, + # npts_max=npts_max, + # nth_sample_new=nth_sample_new, + # npts2plot=npts2plot)), + # QMessageBox.Yes | QMessageBox.No, + # QMessageBox.Yes) + # if ans == QMessageBox.Yes: + # settings.setValue("nth_sample", nth_sample_new) + # settings.sync() def get_npts_to_plot(self): return sum(trace.stats.npts for trace in self.data.getWFData()) @@ -1740,9 +1758,12 @@ class MainWindow(QMainWindow): def finish_pg_plot(self): self.getPlotWidget().updateWidget() plots = self.wfp_thread.data - for times, data in plots: + for times, data, times_syn, data_syn in plots: self.dataPlot.plotWidget.getPlotItem().plot(times, data, pen=self.dataPlot.pen_linecolor) + if len(data_syn) > 0: + self.dataPlot.plotWidget.getPlotItem().plot(times_syn, data_syn, + pen=self.dataPlot.pen_linecolor_syn) self.dataPlot.reinitMoveProxy() self.dataPlot.plotWidget.showAxis('left') self.dataPlot.plotWidget.showAxis('bottom') @@ -1870,6 +1891,7 @@ class MainWindow(QMainWindow): comp = self.getComponent() title = 'section: {0} components'.format(zne_text[comp]) wfst = self.get_data().getWFData() + wfsyn = self.get_data().getSynWFData() if self.filterActionP.isChecked() and filter: self.filterWaveformData(plot=False, phase='P') elif self.filterActionS.isChecked() and filter: @@ -1878,8 +1900,12 @@ class MainWindow(QMainWindow): # wfst += self.get_data().getWFData().select(component=alter_comp) plotWidget = self.getPlotWidget() self.adjustPlotHeight() - plots = plotWidget.plotWFData(wfdata=wfst, title=title, mapping=False, component=comp, - nth_sample=int(nth_sample)) + if real_Bool(settings.value('large_dataset')) == True: + self.plot_method = 'fast' + else: + self.plot_method = 'normal' + plots = plotWidget.plotWFData(wfdata=wfst, wfsyn=wfsyn, title=title, mapping=False, component=comp, + nth_sample=int(nth_sample), method=self.plot_method) return plots def adjustPlotHeight(self): @@ -3153,6 +3179,10 @@ class MainWindow(QMainWindow): def draw(self): self.fill_eventbox() self.getPlotWidget().draw() + if self.plot_method == 'fast': + self.dataPlot.setPermText('MIN/MAX plot', color='red') + else: + self.dataPlot.setPermText() def _setDirty(self): self.setDirty(True) diff --git a/pylot/core/io/data.py b/pylot/core/io/data.py index 60d017d1..a7cb85d5 100644 --- a/pylot/core/io/data.py +++ b/pylot/core/io/data.py @@ -368,7 +368,7 @@ class Data(object): data.filter(**kwargs) self.dirty = True - def setWFData(self, fnames, checkRotated=False, metadata=None): + def setWFData(self, fnames, checkRotated=False, metadata=None, obspy_dmt=False): """ Clear current waveform data and set given waveform data :param fnames: waveform data names to append @@ -376,8 +376,22 @@ class Data(object): """ self.wfdata = Stream() self.wforiginal = None - if fnames is not None: - self.appendWFData(fnames) + self.wfsyn = Stream() + wffnames = None + wffnames_syn = None + wfdir = 'processed' if 'processed' in [fname.split('/')[-1] for fname in fnames] else 'raw' + if obspy_dmt: + for fpath in fnames: + if fpath.endswith(wfdir): + wffnames = [os.path.join(fpath, fname) for fname in os.listdir(fpath)] + if 'syngine' in fpath.split('/')[-1]: + wffnames_syn = [os.path.join(fpath, fname) for fname in os.listdir(fpath)] + else: + wffnames = fnames + if wffnames is not None: + self.appendWFData(wffnames) + if wffnames_syn is not None: + self.appendWFData(wffnames_syn, synthetic=True) else: return False @@ -399,8 +413,7 @@ class Data(object): return True - - def appendWFData(self, fnames): + def appendWFData(self, fnames, synthetic=False): """ Read waveform data from fnames and append it to current wf data :param fnames: waveform data to append @@ -413,13 +426,16 @@ class Data(object): if self.dirty: self.resetWFData() + real_or_syn_data = {True: self.wfsyn, + False: self.wfdata} + warnmsg = '' for fname in fnames: try: - self.wfdata += read(fname) + real_or_syn_data[synthetic] += read(fname) except TypeError: try: - self.wfdata += read(fname, format='GSE2') + real_or_syn_data[synthetic] += read(fname, format='GSE2') except Exception as e: warnmsg += '{0}\n{1}\n'.format(fname, e) except SacIOError as se: @@ -434,6 +450,9 @@ class Data(object): def getOriginalWFData(self): return self.wforiginal + def getSynWFData(self): + return self.wfsyn + def resetWFData(self): """ Set waveform data to original waveform data diff --git a/pylot/core/util/dataprocessing.py b/pylot/core/util/dataprocessing.py index 1d13e952..6b952164 100644 --- a/pylot/core/util/dataprocessing.py +++ b/pylot/core/util/dataprocessing.py @@ -192,7 +192,7 @@ def read_metadata(path_to_inventory): robj = inv[invtype] else: print("Reading metadata information from inventory-xml file ...") - robj = inv[invtype] + robj = read_inventory(inv[invtype]) return invtype, robj diff --git a/pylot/core/util/obspyDMT_interface.py b/pylot/core/util/obspyDMT_interface.py new file mode 100644 index 00000000..661e4f5b --- /dev/null +++ b/pylot/core/util/obspyDMT_interface.py @@ -0,0 +1,28 @@ +#!/usr/bin/env python +# -*- coding: utf-8 -*- + +import os +from obspy import UTCDateTime + +def check_obspydmt_structure(path): + ''' + Check path for obspyDMT event structure. + :param path: + :return: + ''' + ev_info = os.path.join(path, 'EVENTS-INFO') + if os.path.isdir(ev_info): + if os.path.isfile(os.path.join(ev_info, 'logger_command.txt')): + return True + return False + +def check_obspydmt_eventfolder(folder): + try: + time = folder.split('.')[0] + time = time.replace('_', 'T') + time = UTCDateTime(time) + return True, time + except Exception as e: + return False, e + +check_obspydmt_eventfolder('20110311_054623.a') \ No newline at end of file diff --git a/pylot/core/util/structure.py b/pylot/core/util/structure.py index 68a16552..f556eb7c 100644 --- a/pylot/core/util/structure.py +++ b/pylot/core/util/structure.py @@ -9,4 +9,4 @@ Created on Wed Jan 26 17:47:25 2015 from pylot.core.io.data import SeiscompDataStructure, PilotDataStructure DATASTRUCTURE = {'PILOT': PilotDataStructure, 'SeisComP': SeiscompDataStructure, - None: None} + 'obspyDMT': None, None: None} diff --git a/pylot/core/util/utils.py b/pylot/core/util/utils.py index 2c93286e..79872e31 100644 --- a/pylot/core/util/utils.py +++ b/pylot/core/util/utils.py @@ -6,6 +6,7 @@ import os import platform import re import subprocess +import warnings import numpy as np from obspy import UTCDateTime, read @@ -13,6 +14,8 @@ from obspy.core import AttribDict from obspy.signal.rotate import rotate2zne from obspy.io.xseed.utils import SEEDParserException +from pylot.core.util.obspyDMT_interface import check_obspydmt_eventfolder + from pylot.core.io.inputs import PylotParameter, FilterOptions from pylot.styles import style_settings @@ -589,7 +592,7 @@ def prepTimeAxis(stime, trace, verbosity=0): srate = trace.stats.sampling_rate tincr = trace.stats.delta etime = stime + nsamp / srate - time_ax = np.arange(stime, etime, tincr) + time_ax = np.linspace(stime, etime, nsamp) if len(time_ax) < nsamp: if verbosity: print('elongate time axes by one datum') @@ -1216,6 +1219,46 @@ def has_spe(pick): return pick['spe'] +def check_all_obspy(eventlist): + ev_type = 'obspydmt' + return check_event_folders(eventlist, ev_type) + + +def check_all_pylot(eventlist): + ev_type = 'pylot' + return check_event_folders(eventlist, ev_type) + + +def check_event_folders(eventlist, ev_type): + checklist = [] + clean_eventlist = [] + for path in eventlist: + folder_check = check_event_folder(path) + if not folder_check: + warnings.warn('Unrecognized event folder: {}'.format(path)) + continue + checklist.append(folder_check == ev_type) + clean_eventlist.append(path) + if all(checklist) or len(checklist) == 0: + return clean_eventlist + else: + warnings.warn('Not all selected folders of type {}'.format(ev_type)) + return [] + + +def check_event_folder(path): + ev_type = None + folder = path.split('/')[-1] + # for pylot: select only folders that start with 'e', containin two dots and have length 12 + if (folder.startswith('e') + and len(folder.split('.')) == 3 + and len(folder) == 12): + ev_type = 'pylot' + elif check_obspydmt_eventfolder(folder)[0]: + ev_type = 'obspydmt' + return ev_type + + if __name__ == "__main__": import doctest diff --git a/pylot/core/util/widgets.py b/pylot/core/util/widgets.py index 04092146..56dc16dd 100644 --- a/pylot/core/util/widgets.py +++ b/pylot/core/util/widgets.py @@ -25,6 +25,7 @@ except ImportError: from matplotlib.backends.backend_qt4agg import NavigationToolbar2QT from matplotlib.widgets import MultiCursor from matplotlib.tight_layout import get_renderer, get_subplotspec_list, get_tight_layout_figure +from scipy.signal import argrelmin, argrelmax from PySide import QtCore, QtGui from PySide.QtGui import QAction, QApplication, QCheckBox, QComboBox, \ @@ -34,7 +35,7 @@ from PySide.QtGui import QAction, QApplication, QCheckBox, QComboBox, \ QPushButton, QFileDialog, QInputDialog, QKeySequence from PySide.QtCore import QSettings, Qt, QUrl, Signal, Slot from PySide.QtWebKit import QWebView -from obspy import Stream, UTCDateTime +from obspy import Stream, Trace, UTCDateTime from obspy.core.util import AttribDict from obspy.taup import TauPyModel from obspy.taup.utils import get_phase_names @@ -448,20 +449,25 @@ class WaveformWidgetPG(QtGui.QWidget): self.plotdict = dict() # create plot self.main_layout = QtGui.QVBoxLayout() - self.label = QtGui.QLabel() + self.label_layout = QtGui.QHBoxLayout() + self.status_label = QtGui.QLabel() + self.perm_label = QtGui.QLabel() self.setLayout(self.main_layout) self.plotWidget = self.pg.PlotWidget(self.parent(), title=title) self.main_layout.addWidget(self.plotWidget) - self.main_layout.addWidget(self.label) + self.main_layout.addLayout(self.label_layout) + self.label_layout.addWidget(self.status_label) + self.label_layout.addWidget(self.perm_label) self.plotWidget.showGrid(x=False, y=True, alpha=0.3) self.plotWidget.hideAxis('bottom') self.plotWidget.hideAxis('left') self.wfstart, self.wfend = 0, 0 self.pen_multicursor = self.pg.mkPen(self.parent()._style['multicursor']['rgba']) self.pen_linecolor = self.pg.mkPen(self.parent()._style['linecolor']['rgba']) + self.pen_linecolor_syn = self.pg.mkPen((100, 0, 255, 255)) self.reinitMoveProxy() self._proxy = self.pg.SignalProxy(self.plotWidget.scene().sigMouseMoved, rateLimit=60, slot=self.mouseMoved) - self.plotWidget.getPlotItem().setDownsampling(auto=True) + #self.plotWidget.getPlotItem().setDownsampling(auto=True) def reinitMoveProxy(self): self.vLine = self.pg.InfiniteLine(angle=90, movable=False, pen=self.pen_multicursor) @@ -479,22 +485,27 @@ class WaveformWidgetPG(QtGui.QWidget): station = self.orig_parent.getStationName(wfID) abstime = self.wfstart + x if self.orig_parent.get_current_event(): - self.label.setText("station = {}, T = {}, t = {} [s]".format(station, abstime, x)) + self.status_label.setText("station = {}, T = {}, t = {} [s]".format(station, abstime, x)) self.vLine.setPos(mousePoint.x()) self.hLine.setPos(mousePoint.y()) def getPlotDict(self): return self.plotdict + def setPermText(self, text=None, color='black'): + self.perm_label.setText(text) + self.perm_label.setStyleSheet('color: {}'.format(color)) + def setPlotDict(self, key, value): self.plotdict[key] = value def clearPlotDict(self): self.plotdict = dict() - def plotWFData(self, wfdata, title=None, zoomx=None, zoomy=None, + def plotWFData(self, wfdata, wfsyn=None, title=None, zoomx=None, zoomy=None, noiselevel=None, scaleddata=False, mapping=True, - component='*', nth_sample=1, iniPick=None, verbosity=0): + component='*', nth_sample=1, iniPick=None, verbosity=0, + method='normal'): if not wfdata: print('Nothing to plot.') return @@ -537,7 +548,10 @@ class WaveformWidgetPG(QtGui.QWidget): for n, (network, station, channel) in enumerate(nsc): n+=1 st = st_select.select(network=network, station=station, channel=channel) - trace = st[0] + trace = st[0].copy() + st_syn = wfsyn.select(network=network, station=station, channel=channel) + if st_syn: + trace_syn = st_syn[0].copy() if mapping: comp = channel[-1] n = compclass.getPlotPosition(str(comp)) @@ -549,13 +563,28 @@ class WaveformWidgetPG(QtGui.QWidget): print(msg) stime = trace.stats.starttime - self.wfstart time_ax = prepTimeAxis(stime, trace) - if time_ax is not None: + if st_syn: + stime_syn = trace_syn.stats.starttime - self.wfstart + time_ax_syn = prepTimeAxis(stime_syn, trace_syn) + + if method == 'fast': + trace.data, time_ax = self.minMax(trace, time_ax) + + if time_ax not in [None, []]: if not scaleddata: trace.detrend('constant') trace.normalize(np.max(np.abs(trace.data)) * 2) - times = [time for index, time in enumerate(time_ax) if not index % nth_sample] - data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample] - plots.append((times, data)) + if st_syn: + trace_syn.detrend('constant') + trace_syn.normalize(np.max(np.abs(trace_syn.data)) * 2) + # TODO: change this to numpy operations instead of lists? + times = np.array([time for index, time in enumerate(time_ax) if not index % nth_sample]) + times_syn = np.array([time for index, time in enumerate(time_ax_syn) if not index % nth_sample] if st_syn else []) + trace.data = np.array([datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]) + trace.data_syn = np.array([datum + n for index, datum in enumerate(trace.data_syn) + if not index % nth_sample] if st_syn else []) + plots.append((times, trace.data, + times_syn, trace.data_syn)) self.setPlotDict(n, (station, channel, network)) self.xlabel = 'seconds since {0}'.format(self.wfstart) self.ylabel = '' @@ -563,6 +592,41 @@ class WaveformWidgetPG(QtGui.QWidget): self.setYLims([0.5, nmax + 0.5]) return plots + def minMax(self, trace, time_ax): + ''' + create min/max array for fast plotting (approach based on obspy __plot_min_max function) + :returns data, time_ax + ''' + npixel = self.width() + ndata = len(trace.data) + pts_per_pixel = ndata/npixel + if pts_per_pixel < 2: + return trace.data, time_ax + remaining_samples = ndata%pts_per_pixel + npixel = ndata//pts_per_pixel + if remaining_samples: + data = trace.data[:-remaining_samples] + else: + data = trace.data + data = data.reshape(npixel, pts_per_pixel) + min_ = data.min(axis=1) + max_ = data.max(axis=1) + if remaining_samples: + extreme_values = np.empty((npixel + 1, 2), dtype=np.float) + extreme_values[:-1, 0] = min_ + extreme_values[:-1, 1] = max_ + extreme_values[-1, 0] = \ + trace.data[-remaining_samples:].min() + extreme_values[-1, 1] = \ + trace.data[-remaining_samples:].max() + else: + extreme_values = np.empty((npixel, 2), dtype=np.float) + extreme_values[:, 0] = min_ + extreme_values[:, 1] = max_ + data = extreme_values.flatten() + time_ax = np.linspace(time_ax[0], time_ax[-1], num=len(data)) + return data, time_ax + # def getAxes(self): # return self.axes