Merge branch 'feature/obspy_dmt_interface' into develop

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 eventpath:
-            if not directory:
                 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
+            basepath = eventlist[0].split(os.path.basename(eventlist[0]))[0]
-            eventlist = [item for item in eventlist if item.split('/')[-1].startswith('e')
+            if check_obspydmt_structure(basepath):
-                         and len(item.split('/')[-1].split('.')) == 3
+                print('Recognized obspyDMT structure in selected files.')
-                         and len(item.split('/')[-1]) == 12]
+                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)
+        nth_sample = int(settings.value("nth_sample")) if settings.value("nth_sample") else 1
-        npts_max = 1e6
+        npts_max = 1e7
         npts = self.get_npts_to_plot()
         npts2plot = npts/nth_sample
         if npts2plot < npts_max:
-            return
+            settings.setValue('large_dataset', False)
-        nth_sample_new = int(np.ceil(npts/npts_max))
+        else:
-        message = "You are about to plot a huge dataset with {npts} datapoints. With a current setting of " \
+            settings.setValue('large_dataset', True)
-                  "nth_sample = {nth_sample} a total of {npts2plot} points will be plotted which is more " \
+            self.update_status('Dataset is very large. Using fast plotting method (MIN/MAX)', 10000)
-                  "than the maximum setting of {npts_max}. " \
+        settings.sync()
-                  "PyLoT recommends to raise nth_sample from {nth_sample} to {nth_sample_new}. Continue?"
+        # 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 " \
-        ans = QMessageBox.question(self, self.tr("Optimize plot performance..."),
+        #           "nth_sample = {nth_sample} a total of {npts2plot} points will be plotted which is more " \
-                                   self.tr(message.format(npts=npts,
+        #           "than the maximum setting of {npts_max}. " \
-                                                          nth_sample=nth_sample,
+        #           "PyLoT recommends to raise nth_sample from {nth_sample} to {nth_sample_new}. Do you want "\
-                                                          npts_max=npts_max,
+        #           "to change nth_sample to {nth_sample_new} now?"
-                                                          nth_sample_new=nth_sample_new,
+        #
-                                                          npts2plot=npts2plot)),
+        # ans = QMessageBox.question(self, self.tr("Optimize plot performance..."),
-                                   QMessageBox.Yes | QMessageBox.No,
+        #                            self.tr(message.format(npts=npts,
-                                   QMessageBox.Yes)
+        #                                                   nth_sample=nth_sample,
-        if ans == QMessageBox.Yes:
+        #                                                   npts_max=npts_max,
-            settings.setValue("nth_sample", nth_sample_new)
+        #                                                   nth_sample_new=nth_sample_new,
-            settings.sync()
+        #                                                   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,
+        if real_Bool(settings.value('large_dataset')) == True:
-                                      nth_sample=int(nth_sample))
+            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)

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.wfsyn = Stream()
-            self.appendWFData(fnames)
+        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, synthetic=False):
-    def appendWFData(self, fnames):
         """
         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

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
 
 

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')

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}

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
 

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]
+                    if st_syn:
-                data = [datum + n for index, datum in enumerate(trace.data) if not index % nth_sample]
+                        trace_syn.detrend('constant')
-                plots.append((times, data))
+                        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