[update] improved array map to identify and display other phase names than P / S

[minor] some tweaks and optimisations

PyLoT.py

@@ -76,7 +76,7 @@ from pylot.core.util.utils import fnConstructor, getLogin, \
     full_range, readFilterInformation, pick_color_plt, \
     pick_linestyle_plt, identifyPhaseID, excludeQualityClasses, \
     transform_colors_mpl, transform_colors_mpl_str, getAutoFilteroptions, check_all_obspy, \
-    check_all_pylot, get_bool, get_None
+    check_all_pylot, get_bool, get_None, get_pylot_eventfile_with_extension
 from pylot.core.util.gui import make_pen
 from pylot.core.util.event import Event
 from pylot.core.io.location import create_creation_info, create_event
@@ -84,7 +84,7 @@ from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \
     PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \
     ComparisonWidget, TuneAutopicker, PylotParaBox, AutoPickDlg, CanvasWidget, AutoPickWidget, \
     CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget, PickQualitiesFromXml, \
-    SourceSpecWindow, ChooseWaveFormWindow, SpectrogramTab
+    SourceSpecWindow, ChooseWaveFormWindow, SpectrogramTab, SearchFileByExtensionDialog
 from pylot.core.util.array_map import Array_map
 from pylot.core.util.structure import DATASTRUCTURE
 from pylot.core.util.thread import Thread, Worker
@@ -1006,16 +1006,15 @@ class MainWindow(QMainWindow):
             return
         refresh = False
         events = self.project.eventlist
-        sld = SingleTextLineDialog(label='Specify file extension: ', default_text='.xml')
+        sld = SearchFileByExtensionDialog(label='Specify file extension: ', default_text='.xml',
+                                          events=events)
         if not sld.exec_():
             return
         fext = sld.lineEdit.text()
         # fext = '.xml'
         for event in events:
-            path = event.path
-            eventname = path.split('/')[-1]  # or event.pylot_id
-            filename = os.path.join(path, 'PyLoT_' + eventname + fext)
-            if os.path.isfile(filename):
+            filename = get_pylot_eventfile_with_extension(event, fext)
+            if filename:
                 self.load_data(filename, draw=False, event=event, overwrite=True)
                 refresh = True
         if not refresh:
@@ -1887,6 +1886,7 @@ class MainWindow(QMainWindow):
         # which will read in data input twice. Therefore current tab is changed to 0
         # in loadProject before calling this function.
         self.fill_eventbox()
+        #print(f'{self.get_current_event()=}')
         plotted = False
         if self.tabs.currentIndex() == 2:
             self.init_event_table()
@@ -1921,7 +1921,6 @@ class MainWindow(QMainWindow):
             self.spectro_layout.addWidget(newSpectroWidget)
             self.spectroWidget = newSpectroWidget
 
-
     def newWF(self, event=None, plot=True):
         '''
         Load new data and plot if necessary.
@@ -1972,6 +1971,8 @@ class MainWindow(QMainWindow):
         #     ans = False
 
         settings = QSettings()
+        # process application events to wait for event items to appear in event box
+        QApplication.processEvents()
         curr_event = self.get_current_event()
         if not curr_event:
             print('Could not find current event. Try reload?')
@@ -2127,7 +2128,7 @@ class MainWindow(QMainWindow):
 
     def finish_pg_plot(self):
         self.getPlotWidget().updateWidget()
-        plots, gaps = self.wfp_thread.data
+        plots = self.wfp_thread.data
         # do not show plot if no data are given
         self.wf_scroll_area.setVisible(len(plots) > 0)
         self.no_data_label.setVisible(not len(plots) > 0)
@@ -2301,8 +2302,8 @@ class MainWindow(QMainWindow):
             self.plot_method = 'normal'
         rval = plotWidget.plotWFData(wfdata=wfst, wfsyn=wfsyn, title=title, mapping=False, component=comp,
                                      nth_sample=int(nth_sample), method=self.plot_method, gain=self.gain)
-        plots, gaps = rval if rval else ([], [])
-        return plots, gaps
+        plots = rval if rval else []
+        return plots
 
     def adjustPlotHeight(self):
         if self.pg:
@@ -3633,7 +3634,7 @@ class MainWindow(QMainWindow):
             return True
         return False
 
-    def update_status(self, message, duration=5000):
+    def update_status(self, message, duration=10000):
         self.statusBar().showMessage(message, duration)
         if self.get_data() is not None:
             if not self.get_current_event() or not self.project.location:
@@ -3807,7 +3808,8 @@ class MainWindow(QMainWindow):
 
     def closeEvent(self, event):
         if self.okToContinue():
-            self.logwidget.close()
+            if hasattr(self, 'logwidget'):
+                self.logwidget.close()
             event.accept()
         else:
             event.ignore()

autopylot.sh

@@ -9,4 +9,4 @@
 
 conda activate pylot_38
 
-python ./autoPyLoT.py -i /home/marcel/.pylot/pylot_janis_noisy.in -c $NSLOTS
+python ./autoPyLoT.py -i /home/marcel/.pylot/pylot_adriaarray.in -c 20 -dmt processed

pylot/core/io/data.py

@@ -19,7 +19,7 @@ from pylot.core.util.errors import FormatError, OverwriteError
 from pylot.core.util.event import Event
 from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
 from pylot.core.util.utils import fnConstructor, full_range, check4rotated, \
-    check4gapsAndMerge, trim_station_components
+    check_for_gaps_and_merge, trim_station_components, check_for_nan
 
 
 class Data(object):
@@ -64,7 +64,7 @@ class Data(object):
                     elif 'LOC' in evtdata:
                         raise NotImplementedError('PILOT location information '
                                                   'read support not yet '
-                                                  'implemeted.')
+                                                  'implemented.')
                     elif 'event.pkl' in evtdata:
                         evtdata = qml_from_obspyDMT(evtdata)
                     else:
@@ -457,6 +457,11 @@ class Data(object):
         :param fnames: waveform data names to append
         :type fnames: list
         """
+        def check_fname_exists(filenames: list) -> list:
+            if filenames:
+                filenames = [fn for fn in filenames if os.path.isfile(fn)]
+            return filenames
+
         self.wfdata = Stream()
         self.wforiginal = None
         self.wfsyn = Stream()
@@ -465,6 +470,8 @@ class Data(object):
         self.tstart = tstart
         self.tstop = tstop
 
+        fnames = check_fname_exists(fnames)
+        fnames_syn = check_fname_exists(fnames_syn)
         # if obspy_dmt:
         #     wfdir = 'raw'
         #     self.processed = False
@@ -491,7 +498,9 @@ class Data(object):
         # remove possible underscores in station names
         # self.wfdata = remove_underscores(self.wfdata)
         # check for gaps and merge
-        self.wfdata = check4gapsAndMerge(self.wfdata)
+        self.wfdata, _ = check_for_gaps_and_merge(self.wfdata)
+        # check for nans
+        check_for_nan(self.wfdata)
         # check for stations with rotated components
         if checkRotated and metadata is not None:
             self.wfdata = check4rotated(self.wfdata, metadata, verbosity=0)

pylot/core/pick/autopick.py

@@ -660,7 +660,7 @@ class AutopickStation(object):
             ax1.set_ylim([-1.5, 1.5])
             ax1.set_ylabel('Normalized Counts')
 
-            if self.horizontal_traces_exist() and self.s_data.Sflag == 1:
+            if self.horizontal_traces_exist():# and self.s_data.Sflag == 1:
                 # plot E trace
                 ax2 = fig.add_subplot(3, 1, 2, sharex=ax1)
                 th1data = np.linspace(0, self.etrace.stats.endtime - self.etrace.stats.starttime,

pylot/core/pick/utils.py

@@ -1332,22 +1332,6 @@ def get_quality_class(uncertainty, weight_classes):
     return quality
 
 
-def set_NaNs_to(data, nan_value):
-    """
-    Replace all NaNs in data with nan_value
-    :param data: array holding data
-    :type data: `~numpy.ndarray`
-    :param nan_value: value which all NaNs are set to
-    :type nan_value: float, int
-    :return: data array with all NaNs replaced with nan_value
-    :rtype: `~numpy.ndarray`
-    """
-    nn = np.isnan(data)
-    if np.any(nn):
-        data[nn] = nan_value
-    return data
-
-
 def taper_cf(cf):
     """
     Taper cf data to get rid off of side maximas

pylot/core/util/array_map.py

@@ -13,6 +13,7 @@ import obspy
 from PySide2 import QtWidgets
 from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+from pylot.core.util.utils import identifyPhaseID
 from scipy.interpolate import griddata
 
 from pylot.core.pick.utils import get_quality_class
@@ -279,9 +280,12 @@ class Array_map(QtWidgets.QWidget):
         self.canvas.axes.figure.canvas.draw_idle()
 
     def onpick(self, event):
+        btn_msg = {1: ' in selection. Aborted', 2: ' to delete a pick on. Aborted', 3: ' to display info.'}
         ind = event.ind
         button = event.mouseevent.button
-        if ind == []:
+        msg_reason = None
+        if len(ind) > 1:
+            self._parent.update_status(f'Found more than one station {btn_msg.get(button)}')
             return
         if button == 1:
             self.openPickDlg(ind)
@@ -384,7 +388,14 @@ class Array_map(QtWidgets.QWidget):
                 try:
                     station_name = st_id.split('.')[-1]
                     # current_picks_dict: auto or manual
-                    pick = self.current_picks_dict()[station_name][phase]
+                    station_picks = self.current_picks_dict().get(station_name)
+                    if not station_picks:
+                        continue
+                    for phase_hint, pick in station_picks.items():
+                        if identifyPhaseID(phase_hint) == phase:
+                            break
+                    else:
+                        continue
                     if pick['picker'] == 'auto':
                         if not pick['spe']:
                             continue

pylot/core/util/dataprocessing.py

@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 
 import glob
+import logging
 import os
 import sys
 
@@ -189,7 +190,11 @@ class Metadata(object):
         metadata = self.get_metadata(seed_id, time)
         if not metadata:
             return
-        return metadata['data'].get_coordinates(seed_id, time)
+        try:
+            return metadata['data'].get_coordinates(seed_id, time)
+        # no specific exception defined in obspy inventory
+        except Exception as e:
+            logging.warning(f'Could not get metadata for {seed_id}')
 
     def get_all_coordinates(self):
         def stat_info_from_parser(parser):

pylot/core/util/generate_array_maps.py

@@ -42,7 +42,7 @@ def main(project_file_path, manual=False, auto=True, file_format='png', f_ext=''
         for item in input_list:
             array_map_worker(item)
     else:
-        pool = multiprocessing.Pool(ncores)
+        pool = multiprocessing.Pool(ncores, maxtasksperchild=1000)
         pool.map(array_map_worker, input_list)
         pool.close()
         pool.join()

pylot/core/util/thread.py

@@ -160,7 +160,7 @@ class MultiThread(QThread):
         try:
             if not self.ncores:
                 self.ncores = multiprocessing.cpu_count()
-            pool = multiprocessing.Pool(self.ncores)
+            pool = multiprocessing.Pool(self.ncores, maxtasksperchild=1000)
             self.data = pool.map_async(self.func, self.args, callback=self.emitDone)
             # self.data = pool.apply_async(self.func, self.shotlist, callback=self.emitDone) #emit each time returned
             pool.close()

pylot/core/util/utils.py

@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 
 import hashlib
+import logging
 import os
 import platform
 import re
@@ -119,7 +120,7 @@ def gen_Pool(ncores=0):
 
     print('gen_Pool: Generated multiprocessing Pool with {} cores\n'.format(ncores))
 
-    pool = multiprocessing.Pool(ncores)
+    pool = multiprocessing.Pool(ncores, maxtasksperchild=100)
     return pool
 
 
@@ -818,19 +819,6 @@ def trim_station_components(data, trim_start=True, trim_end=True):
     return data
 
 
-def merge_stream(stream):
-    gaps = stream.get_gaps()
-    if gaps:
-        # list of merged stations (seed_ids)
-        merged = ['{}.{}.{}.{}'.format(*gap[:4]) for gap in gaps]
-        stream.merge(method=1)
-        print('Merged the following stations because of gaps:')
-        for merged_station in merged:
-            print(merged_station)
-
-    return stream, gaps
-
-
 def check4gapsAndRemove(data):
     """
     check for gaps in Stream and remove them
@@ -851,12 +839,12 @@ def check4gapsAndRemove(data):
     return data
 
 
-def check4gapsAndMerge(data):
+def check_for_gaps_and_merge(data):
     """
     check for gaps in Stream and merge if gaps are found
     :param data: stream of seismic data
     :type data: `~obspy.core.stream.Stream`
-    :return: data stream
+    :return: data stream, gaps returned from obspy get_gaps
     :rtype: `~obspy.core.stream.Stream`
     """
     gaps = data.get_gaps()
@@ -867,7 +855,7 @@ def check4gapsAndMerge(data):
         for merged_station in merged:
             print(merged_station)
 
-    return data
+    return data, gaps
 
 
 def check4doubled(data):
@@ -897,6 +885,33 @@ def check4doubled(data):
     return data
 
 
+def check_for_nan(data, nan_value=0.):
+    """
+    Replace all NaNs in data with nan_value (in place)
+    :param data: stream of seismic data
+    :type data: `~obspy.core.stream.Stream`
+    :param nan_value: value which all NaNs are set to
+    :type nan_value: float, int
+    :return: None
+    """
+    if not data:
+        return
+    for trace in data:
+        np.nan_to_num(trace.data, copy=False, nan=nan_value)
+
+
+def get_pylot_eventfile_with_extension(event, fext):
+    if hasattr(event, 'path'):
+        eventpath = event.path
+    else:
+        logging.warning('No attribute path found for event.')
+        return
+    eventname = eventpath.split('/')[-1]  # or event.pylot_id
+    filename = os.path.join(eventpath, 'PyLoT_' + eventname + fext)
+    if os.path.isfile(filename):
+        return filename
+
+
 def get_stations(data):
     """
     Get list of all station names in data stream

pylot/core/util/widgets.py

@@ -36,7 +36,7 @@ from PySide2.QtWidgets import QAction, QApplication, QCheckBox, QComboBox, \
     QGridLayout, QLabel, QLineEdit, QMessageBox, \
     QTabWidget, QToolBar, QVBoxLayout, QHBoxLayout, QWidget, \
     QPushButton, QFileDialog, QInputDialog
-from PySide2.QtCore import QSettings, Qt, QUrl, Signal
+from PySide2.QtCore import QSettings, Qt, QUrl, Signal, QTimer
 from PySide2.QtWebEngineWidgets import QWebEngineView as QWebView
 from obspy import Stream, Trace, UTCDateTime
 from obspy.core.util import AttribDict
@@ -51,10 +51,10 @@ from pylot.core.pick.autopick import fmpicker
 from pylot.core.util.defaults import OUTPUTFORMATS, FILTERDEFAULTS
 from pylot.core.util.utils import prepTimeAxis, full_range, demeanTrace, isSorted, findComboBoxIndex, clims, \
     pick_linestyle_plt, pick_color_plt, \
-    check4rotated, check4doubled, merge_stream, identifyPhase, \
+    check4rotated, check4doubled, check_for_gaps_and_merge, check_for_nan, identifyPhase, \
     loopIdentifyPhase, trim_station_components, transformFilteroptions2String, \
     identifyPhaseID, get_bool, get_None, pick_color, getAutoFilteroptions, SetChannelComponents, \
-    station_id_remove_channel
+    station_id_remove_channel, get_pylot_eventfile_with_extension
 from autoPyLoT import autoPyLoT
 from pylot.core.util.thread import Thread
 from pylot.core.util.dataprocessing import Metadata
@@ -860,13 +860,26 @@ class WaveformWidgetPG(QtWidgets.QWidget):
     def clearPlotDict(self):
         self.plotdict = dict()
 
-    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,
-                   method='normal', gain=1.):
+    def plotWFData(self, wfdata, wfsyn=None, title=None, scaleddata=False, mapping=True,
+                   component='*', nth_sample=1, verbosity=0, method='normal', gain=1., shift_syn=0.2):
+        def station_sort(nslc):
+            """Try to sort after station integer in case of a line array (e.g. active seismics)"""
+            try:
+                rval = sorted(nslc, key=lambda x: int(x.split('.')[1]))
+                return rval
+            except ValueError as e:
+                # this is the standard case for seismological stations
+                pass
+            except Exception as e:
+                print(f'Sorting by station integer failed with unknown exception: {e}')
+
+            # fallback to default sorting
+            return sorted(nslc)
+
         if not wfdata:
             print('Nothing to plot.')
             return
+
         self.title = title
         self.clearPlotDict()
         self.wfstart, self.wfend = full_range(wfdata)
@@ -884,14 +897,14 @@ class WaveformWidgetPG(QtWidgets.QWidget):
         else:
             st_select = wfdata
 
-        st_select, gaps = merge_stream(st_select)
+        # st_select, gaps = check_for_gaps_and_merge(st_select) #MP MP commented because probably done twice
 
         # list containing tuples of network, station, channel (for sorting)
         nslc = []
         for trace in st_select:
             nslc.append(
                 trace.get_id())  # (trace.stats.network, trace.stats.station, trace.stats.location trace.stats.channel))
-        nslc.sort()
+        nslc = station_sort(nslc)
         nslc.reverse()
         plots = []
 
@@ -955,7 +968,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
         self.ylabel = ''
         self.setXLims([0, self.wfend - self.wfstart])
         self.setYLims([0.5, nmax + 0.5])
-        return plots, gaps
+        return plots
 
     def minMax(self, trace, time_ax):
         '''
@@ -977,7 +990,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
         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 = np.empty((npixel + 1, 2), dtype=float)
             extreme_values[:-1, 0] = min_
             extreme_values[:-1, 1] = max_
             extreme_values[-1, 0] = \
@@ -985,7 +998,7 @@ class WaveformWidgetPG(QtWidgets.QWidget):
             extreme_values[-1, 1] = \
                 trace.data[-remaining_samples:].max()
         else:
-            extreme_values = np.empty((npixel, 2), dtype=np.float)
+            extreme_values = np.empty((npixel, 2), dtype=float)
             extreme_values[:, 0] = min_
             extreme_values[:, 1] = max_
         data = extreme_values.flatten()
@@ -1555,6 +1568,96 @@ class PylotCanvas(FigureCanvas):
         self.draw()
 
 
+class SearchFileByExtensionDialog(QtWidgets.QDialog):
+    def __init__(self, parent=None, label='Text: ', default_text='.xml', events=None):
+        super(SearchFileByExtensionDialog, self).__init__(parent)
+        self.events = events
+        self.filepaths = []
+        self.default_text = default_text
+        self.label = label
+        self.setButtons()
+        self.setupUi()
+        self.connectSignals()
+        self.showPaths()
+        self.refresh_timer = QTimer(self)
+        self.refresh_timer.timeout.connect(self.showPaths)
+        self.refresh_timer.start(10000)
+
+        self.resize(800, 450)
+
+
+    def setupUi(self):
+        self.main_layout = QtWidgets.QVBoxLayout()
+        self.header_layout = QtWidgets.QHBoxLayout()
+        #
+        self.setLayout(self.main_layout)
+
+        # widgets inside the dialog
+        self.textLabel = QtWidgets.QLabel(self.label)
+        self.lineEdit = QtWidgets.QLineEdit(self.default_text)
+
+        # optional search button, currently disabled. List refreshed when text changes
+        self.searchButton = QtWidgets.QPushButton('Search')
+        self.searchButton.setVisible(False)
+
+        self.tableWidget = QtWidgets.QTableWidget()
+        tableWidget = self.tableWidget
+        tableWidget.setColumnCount(2)
+        tableWidget.setRowCount(len(self.events))
+        tableWidget.setHorizontalHeaderLabels(('Filename', 'Last modified'))
+        tableWidget.setEditTriggers(tableWidget.NoEditTriggers)
+        tableWidget.setSortingEnabled(True)
+        header = tableWidget.horizontalHeader()
+        header.setSectionResizeMode(QtWidgets.QHeaderView.ResizeToContents)
+        header.setStretchLastSection(True)
+
+        self.statusText = QtWidgets.QLabel()
+
+        self.header_layout.addWidget(self.textLabel)
+        self.header_layout.addWidget(self.lineEdit)
+        self.header_layout.addWidget(self.searchButton)
+
+        self.main_layout.addLayout(self.header_layout)
+        self.main_layout.addWidget(self.tableWidget)
+        self.main_layout.addWidget(self.statusText)
+        self.main_layout.addWidget(self._buttonbox)
+
+    def showPaths(self):
+        self.filepaths = []
+        fext = self.lineEdit.text()
+        self.tableWidget.clearContents()
+        for index, event in enumerate(self.events):
+            filename = get_pylot_eventfile_with_extension(event, fext)
+            if filename:
+                self.filepaths.append(filename)
+                ts = int(os.path.getmtime(filename))
+
+                # create QTableWidgetItems of filepath and last modification time
+                fname_item = QtWidgets.QTableWidgetItem(f'{filename}')
+                ts_item = QtWidgets.QTableWidgetItem(f'{datetime.datetime.fromtimestamp(ts)}')
+                self.tableWidget.setItem(index, 0, fname_item)
+                self.tableWidget.setItem(index, 1, ts_item)
+
+        # TODO: Idea -> only refresh if table contents changed. Use selection to load only a subset of files
+        if len(self.filepaths) > 0:
+            status_text = f'Found {len(self.filepaths)} eventfiles. Do you want to load them?'
+        else:
+            status_text = 'Did not find any files for specified file mask.'
+        self.statusText.setText(status_text)
+
+
+    def setButtons(self):
+        self._buttonbox = QDialogButtonBox(QDialogButtonBox.Ok |
+                                           QDialogButtonBox.Cancel)
+
+    def connectSignals(self):
+        self._buttonbox.accepted.connect(self.accept)
+        self._buttonbox.rejected.connect(self.reject)
+        self.lineEdit.textChanged.connect(self.showPaths)
+        self.searchButton.clicked.connect(self.showPaths)
+
+
+
 class SingleTextLineDialog(QtWidgets.QDialog):
     def __init__(self, parent=None, label='Text: ', default_text='.xml'):
         super(SingleTextLineDialog, self).__init__(parent)
@@ -3568,8 +3671,9 @@ class TuneAutopicker(QWidget):
             # wfdat = remove_underscores(wfdat)
             # rotate misaligned stations to ZNE
             # check for gaps and doubled channels
-            wfdat, gaps = merge_stream(wfdat)
-            # check4gaps(wfdat)
+            wfdat, _ = check_for_gaps_and_merge(wfdat)
+            # check for nans
+            check_for_nan(wfdat)
             check4doubled(wfdat)
             wfdat = check4rotated(wfdat, self.parent().metadata, verbosity=0)
             # trim station components to same start value
@@ -3724,6 +3828,7 @@ class TuneAutopicker(QWidget):
         st = self.data.getWFData()
         tr = st.select(station=self.get_current_station())[0]
         starttime = tr.stats.starttime
+        # create two lists with figure names and subindices (for subplots) to get the correct axes
         p_axes = [
             ('mainFig', 0),
             ('aicFig', 0),