[initial] first implementation of residual plotting (WIP)

# Conflicts:
#	PyLoT.py
#	README.md
#	pylot/core/util/widgets.py

.gitignore

@@ -2,3 +2,4 @@
 *~
 .idea
 pylot/RELEASE-VERSION
+/tests/test_autopicker/dmt_database_test/

.mailmap

@@ -0,0 +1,39 @@
+Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <Darius_A@web.de>
+Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <darius.arnold@rub.de>
+Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <darius.arnold@ruhr-uni-bochum.de>
+Darius Arnold <Darius.Arnold@ruhr-uni-bochum.de> <mail@dariusarnold.de>
+
+Dennis Wlecklik <dennisw@minos02.geophysik.ruhr-uni-bochum.de>
+
+Jeldrik Gaal <jeldrikgaal@gmail.com>
+
+Kaan Coekerim <kaan.coekerim@ruhr-uni-bochum.de>
+Kaan Coekerim <kaan.coekerim@ruhr-uni-bochum.de> <kaan.coekerim@rub.de>
+
+Ludger Kueperkoch <kueperkoch@igem-energie.de> <kueperkoch@bestec-for-nature.com>
+Ludger Kueperkoch <kueperkoch@igem-energie.de> <ludger@quake2.(none)>
+Ludger Kueperkoch <kueperkoch@igem-energie.de> <ludger@sauron.bestec-for-nature>
+
+Marc S. Boxberg <marc.boxberg@rub.de>
+
+Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel.paffrath@rub.de> 
+Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos01.geophysik.ruhr-uni-bochum.de>
+Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos02.geophysik.ruhr-uni-bochum.de>
+Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@minos25.geophysik.ruhr-uni-bochum.de>
+Marcel Paffrath <marcel.paffrath@ruhr-uni-bochum.de> <marcel@email.com>
+
+Sally Zimmermann <sally.zimmermann@ruhr-uni-bochum.de>
+
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos01.geophysik.ruhr-uni-bochum.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos02.geophysik.ruhr-uni-bochum.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastianw@minos22.geophysik.ruhr-uni-bochum.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling-benatelli@scisys.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling@rub.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <sebastian.wehling@rub.de>
+Sebastian Wehling-Benatelli <sebastian.wehling-benatelli@cgi.com> <DarkBeQst@users.noreply.github.com>
+
+Thomas Moeller <thomas.moeller@rub.de>
+
+Ann-Christin Koch <ann-christin.koch@ruhr-uni-bochum.de> <Ann-Christin.Koch@ruhr-uni-bochum.de>
+
+Sebastian Priebe <sebastian.priebe@rub.de>

PyLoT.py

@@ -25,6 +25,7 @@ https://www.iconfinder.com/iconsets/flavour
 
 import argparse
 import json
+import logging
 import os
 import platform
 import shutil
@@ -72,24 +73,25 @@ from pylot.core.util.errors import DatastructureError, \
     OverwriteError
 from pylot.core.util.connection import checkurl
 from pylot.core.util.dataprocessing import Metadata, restitute_data
-from pylot.core.util.utils import fnConstructor, getLogin, \
+from pylot.core.util.utils import fnConstructor, get_login, \
     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
 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
 from pylot.core.util.widgets import FilterOptionsDialog, NewEventDlg, \
     PylotCanvas, WaveformWidgetPG, PropertiesDlg, HelpForm, createAction, PickDlg, \
-    ComparisonWidget, TuneAutopicker, PylotParaBox, AutoPickDlg, CanvasWidget, AutoPickWidget, \
+    ComparisonWidget, TuneAutopicker, PylotParameterWidget, AutoPickDlg, CanvasWidget, AutoPickWidget, \
-    CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget
+    CompareEventsWidget, ProgressBarWidget, AddMetadataWidget, SingleTextLineDialog, LogWidget, PickQualitiesFromXml, \
+    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
 from pylot.core.util.version import get_git_version as _getVersionString
 from pylot.core.io.getEventListFromXML import geteventlistfromxml
-from pylot.core.io.getQualitiesfromxml import getQualitiesfromxml
+from pylot.core.io.phases import getQualitiesfromxml
 
 from pylot.styles import style_settings
 
@@ -112,21 +114,19 @@ class MainWindow(QMainWindow):
     def __init__(self, parent=None, infile=None, reset_qsettings=False):
         super(MainWindow, self).__init__(parent)
 
-        # check for default pylot.in-file
+        if infile and os.path.isfile(infile) is False:
-        if not infile:
+            infile = QFileDialog().getOpenFileName(caption='Choose PyLoT-input file')[0]
-            infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
-            print('Using default input file {}'.format(infile))
-        if os.path.isfile(infile) == False:
-            infile = QFileDialog().getOpenFileName(caption='Choose PyLoT-input file')
 
-            if not os.path.exists(infile[0]):
+            if not os.path.exists(infile):
                 QMessageBox.warning(self, "PyLoT Warning",
-                                    "No PyLoT-input file declared!")
+                                    "No PyLoT-input file declared! Using default parameters!")
-                sys.exit(0)
+                infile = None
-            self.infile = infile[0]
+
-        else:
-            self.infile = infile
         self._inputs = PylotParameter(infile)
+        if not infile:
+            self._inputs.reset_defaults()
+
+        self.infile = infile
         self._props = None
 
         self.gain = 1.
@@ -136,7 +136,7 @@ class MainWindow(QMainWindow):
         self.project.parameter = self._inputs
         self.tap = None
         self.apw = None
-        self.paraBox = None
+        self.parameterWidget = None
         self.array_map = None
         self._metadata = Metadata(verbosity=0)
         self._eventChanged = [False, False]
@@ -178,6 +178,7 @@ class MainWindow(QMainWindow):
         self.autodata = Data(self)
 
         self.fnames = None
+        self.fnames_comp = None
         self._stime = None
 
         # track deleted picks for logging
@@ -187,13 +188,12 @@ class MainWindow(QMainWindow):
         self.table_headers = ['', 'Event', 'Time', 'Lat', 'Lon', 'Depth', 'Ml', 'Mw', '[N] MP', '[N] AP', 'Tuning Set',
                               'Test Set', 'Notes']
 
-        # TODO: refactor rootpath to datapath
         while True:
             try:
                 if settings.value("user/FullName", None) is None:
                     fulluser = QInputDialog.getText(self, "Enter Name:", "Full name")
                     settings.setValue("user/FullName", fulluser)
-                    settings.setValue("user/Login", getLogin())
+                    settings.setValue("user/Login", get_login())
                 if settings.value("agency_id", None) is None:
                     agency = QInputDialog.getText(self,
                                                   "Enter authority/institution name:",
@@ -250,11 +250,11 @@ class MainWindow(QMainWindow):
             self._inputs.reset_defaults()
             # check for default pylot.in-file
             infile = os.path.join(pylot_config_dir, '.pylot.in')
-            print('Using default input file {}'.format(infile))
+            logging.warning('Using default input file {}'.format(infile))
             self._inputs.export2File(infile)
         self.infile = infile
 
-    def setupUi(self):
+    def setupUi(self, use_logwidget=False):
         try:
             self.startTime = min(
                 [tr.stats.starttime for tr in self.data.wfdata])
@@ -480,7 +480,7 @@ class MainWindow(QMainWindow):
                                                               "automatic pick "
                                                               "data.", False)
         self.compare_action.setEnabled(False)
-        self.qualities_action = self.createAction(parent=self, text='Show pick qualitites...',
+        self.qualities_action = self.createAction(parent=self, text='Show pick qualities...',
                                                   slot=self.pickQualities, shortcut='Alt+Q',
                                                   icon=qualities_icon, tip='Histogram of pick qualities')
         self.qualities_action.setEnabled(False)
@@ -493,7 +493,6 @@ class MainWindow(QMainWindow):
                                                       icon=eventlist_xml_icon,
                                                       tip='Create an Eventlist from a XML File')
         self.eventlist_xml_action.setEnabled(False)
-
         printAction = self.createAction(self, "&Print event ...",
                                         self.show_event_information, QKeySequence.Print,
                                         print_icon,
@@ -506,6 +505,8 @@ class MainWindow(QMainWindow):
         logAction = self.createAction(self, "&Show Log", self.showLogWidget,
                                       tip="""Display Log""")
 
+        logAction.setEnabled(use_logwidget)
+
         # create button group for component selection
 
         componentGroup = QActionGroup(self)
@@ -684,10 +685,9 @@ class MainWindow(QMainWindow):
         # add scroll area used in case number of traces gets too high
         self.wf_scroll_area = QtWidgets.QScrollArea(self)
         self.wf_scroll_area.setVisible(False)
-        self.no_data_label = QLabel('No Data')
+        self.no_data_label = QLabel('No Data. If data were already loaded, try to select the event again in the eventbox.')
         self.no_data_label.setStyleSheet('color: red')
         self.no_data_label.setAlignment(Qt.AlignCenter)
-
         # create central matplotlib figure canvas widget
         self.init_wfWidget()
 
@@ -695,14 +695,17 @@ class MainWindow(QMainWindow):
         wf_tab = QtWidgets.QWidget(self)
         array_tab = QtWidgets.QWidget(self)
         events_tab = QtWidgets.QWidget(self)
+        spectro_tab = QtWidgets.QWidget(self)
 
         # init main widgets layouts
         self.wf_layout = QtWidgets.QVBoxLayout()
         self.array_layout = QtWidgets.QVBoxLayout()
         self.events_layout = QtWidgets.QVBoxLayout()
+        self.spectro_layout = QtWidgets.QVBoxLayout()
         wf_tab.setLayout(self.wf_layout)
         array_tab.setLayout(self.array_layout)
         events_tab.setLayout(self.events_layout)
+        spectro_tab.setLayout(self.spectro_layout)
 
         # tighten up layouts inside tabs
         for layout in [self.wf_layout, self.array_layout, self.events_layout]:
@@ -713,12 +716,14 @@ class MainWindow(QMainWindow):
         self.tabs.addTab(wf_tab, 'Waveform Plot')
         self.tabs.addTab(array_tab, 'Array Map')
         self.tabs.addTab(events_tab, 'Eventlist')
+        #self.tabs.addTab(spectro_tab, 'Spectro')
 
         self.wf_layout.addWidget(self.no_data_label)
         self.wf_layout.addWidget(self.wf_scroll_area)
         self.wf_scroll_area.setWidgetResizable(True)
         self.init_array_tab()
         self.init_event_table()
+        #self.init_spectro_tab()
         self.tabs.setCurrentIndex(0)
 
         self.eventLabel = QLabel()
@@ -731,13 +736,19 @@ class MainWindow(QMainWindow):
         _widget.setLayout(self._main_layout)
         _widget.showFullScreen()
 
-        self.logwidget = LogWidget(parent=None)
+
-        self.logwidget.show()
+        if use_logwidget:
-        sys.stdout = self.logwidget.stdout
+            self.logwidget = LogWidget(parent=None)
-        sys.stderr = self.logwidget.stderr
+            self.logwidget.show()
+            self.stdout = self.logwidget.stdout
+            self.stderr = self.logwidget.stderr
 
         self.setCentralWidget(_widget)
 
+        # Need to store PickQualities Window somewhere so it doesnt disappear
+        self.pickQualitiesWindow = None
+
+
     def init_wfWidget(self):
         xlab = self.startTime.strftime('seconds since %Y/%m/%d %H:%M:%S (%Z)')
         plottitle = None  # "Overview: {0} components ".format(self.getComponent())
@@ -991,18 +1002,17 @@ 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'
+        filenames = sld.getChecked()
         for event in events:
-            path = event.path
+            for filename in filenames:
-            eventname = path.split('/')[-1]  # or event.pylot_id
+                if os.path.isfile(filename) and event.pylot_id in filename:
-            filename = os.path.join(path, 'PyLoT_' + eventname + fext)
+                    self.load_data(filename, draw=False, event=event, ask_user=True, merge_strategy=sld.merge_strategy)
-            if os.path.isfile(filename):
+                    refresh = True
-                self.load_data(filename, draw=False, event=event, overwrite=True)
-                refresh = True
         if not refresh:
             return
         if self.get_current_event().pylot_picks:
@@ -1010,8 +1020,8 @@ class MainWindow(QMainWindow):
         self.fill_eventbox()
         self.setDirty(True)
 
-    def load_data(self, fname=None, loc=False, draw=True, event=None, overwrite=False):
+    def load_data(self, fname=None, loc=False, draw=True, event=None, ask_user=False, merge_strategy='Overwrite'):
-        if not overwrite:
+        if not ask_user:
             if not self.okToContinue():
                 return
         if fname is None:
@@ -1020,14 +1030,33 @@ class MainWindow(QMainWindow):
                 fname = self.filename_from_action(action)
                 if not fname:
                     return
+
         if not event:
             event = self.get_current_event()
+
+        if event.picks:
+            qmb = QMessageBox(self, icon=QMessageBox.Question,
+                              text='Do you want to overwrite the data?',)
+            overwrite_button = qmb.addButton('Overwrite', QMessageBox.YesRole)
+            merge_button = qmb.addButton('Merge', QMessageBox.NoRole)
+            qmb.exec_()
+
+            if qmb.clickedButton() == overwrite_button:
+                merge_strategy = 'Overwrite'
+            elif qmb.clickedButton() == merge_button:
+                merge_strategy = 'Merge'
+            else:
+                return
+
         data = Data(self, event)
         try:
             data_new = Data(self, evtdata=str(fname))
-            # MP MP commented because adding several picks might cause inconsistencies
+            if merge_strategy == 'Overwrite':
-            data = data_new
+                data = data_new
-            # data += data_new
+            elif merge_strategy == 'Merge':
+                data += data_new
+            else:
+                raise NotImplementedError(f'Unknown merge strategy: {merge_strategy}')
         except ValueError:
             qmb = QMessageBox(self, icon=QMessageBox.Question,
                               text='Warning: Missmatch in event identifiers {} and {}. Continue?'.format(
@@ -1115,16 +1144,19 @@ class MainWindow(QMainWindow):
             else:
                 return
 
-    def getWFFnames_from_eventbox(self, eventbox=None):
+    def getWFFnames_from_eventbox(self, eventbox: str = None, subpath: str = None) -> list:
         '''
         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 subpath:
+            eventpath = os.path.join(eventpath, subpath)
+        if not os.path.isdir(eventpath):
+            return []
         if self.dataStructure:
             if not eventpath:
-                return
+                return []
             fnames = [os.path.join(eventpath, f) for f in os.listdir(eventpath)]
         else:
             raise DatastructureError('not specified')
@@ -1168,7 +1200,7 @@ class MainWindow(QMainWindow):
         '''
         if not self.project:
             self.createNewProject()
-        ed = getExistingDirectories(self, 'Select event directories...')
+        ed = GetExistingDirectories(self, 'Select event directories...')
         if ed.exec_():
             eventlist = [event for event in ed.selectedFiles() if not event.endswith('EVENTS-INFO')]
             basepath = eventlist[0].split(os.path.basename(eventlist[0]))[0]
@@ -1178,7 +1210,7 @@ class MainWindow(QMainWindow):
                 with open(eventlist_file, 'r') as infile:
                     eventlist_subset = [os.path.join(basepath, filename.split('\n')[0]) for filename in
                                         infile.readlines()]
-                    msg = 'Found file "eventlist.txt" in database path. WILL ONLY USE SELECTED EVENTS out of {} events ' \
+                    msg = 'Found file "eventlist.txt" in datapath. WILL ONLY USE SELECTED EVENTS out of {} events ' \
                           'contained in this subset'
                     print(msg.format(len(eventlist_subset)))
                     eventlist = [eventname for eventname in eventlist if eventname in eventlist_subset]
@@ -1203,49 +1235,34 @@ class MainWindow(QMainWindow):
         # get path from first event in list and split them
         path = eventlist[0]
         try:
-            system_name = platform.system()
+            datapath = os.path.split(path)[0]
-            if system_name in ["Linux", "Darwin"]:
+            dirs = {
-                dirs = {
+                'datapath': datapath,
-                    'database': path.split('/')[-2],
+            }
-                    'datapath': os.path.split(path)[0],  # path.split('/')[-3],
-                    'rootpath': '/' + os.path.join(*path.split('/')[:-3])
-                }
-            elif system_name == "Windows":
-                rootpath = path.split('/')[:-3]
-                rootpath[0] += '/'
-                dirs = {
-                    # TODO: Arrange path to meet Win standards
-                    'database': path.split('/')[-2],
-                    'datapath': path.split('/')[-3],
-                    'rootpath': os.path.join(*rootpath)
-                }
         except Exception as e:
             dirs = {
-                'database': '',
                 'datapath': '',
-                'rootpath': ''
             }
             print('Warning: Could not automatically init folder structure. ({})'.format(e))
 
         settings = QSettings()
-        settings.setValue("data/dataRoot", dirs['datapath'])  # d irs['rootpath'])
+        settings.setValue("data/dataRoot", dirs['datapath'])
         settings.sync()
 
         if not self.project.eventlist:
             # init parameter object
             self.setParameter(show=False)
             # hide all parameter (show all needed parameter later)
-            self.paraBox.hide_parameter()
+            self.parameterWidget.hide_parameter()
             for directory in dirs.keys():
                 # set parameter
-                box = self.paraBox.boxes[directory]
+                box = self.parameterWidget.boxes[directory]
-                self.paraBox.setValue(box, dirs[directory])
+                self.parameterWidget.setValue(box, dirs[directory])
                 # show needed parameter in box
-                self.paraBox.show_parameter(directory)
+                self.parameterWidget.show_parameter(directory)
-            dirs_box = self.paraBox.get_groupbox_dialog('Directories')
+            dirs_box = self.parameterWidget.get_groupbox_dialog('Directories')
             if not dirs_box.exec_():
                 return
-            self.project.rootpath = dirs['rootpath']
             self.project.datapath = dirs['datapath']
         else:
             if hasattr(self.project, 'datapath'):
@@ -1254,7 +1271,6 @@ class MainWindow(QMainWindow):
                                         'Datapath missmatch to current project!')
                     return
             else:
-                self.project.rootpath = dirs['rootpath']
                 self.project.datapath = dirs['datapath']
 
         self.project.add_eventlist(eventlist)
@@ -1342,11 +1358,10 @@ class MainWindow(QMainWindow):
             return True
 
     def modify_project_path(self, new_rootpath):
-        # TODO: change root to datapath
+        self.project.datapath = new_rootpath
-        self.project.rootpath = new_rootpath
         for event in self.project.eventlist:
-            event.rootpath = new_rootpath
+            event.datapath = new_rootpath
-            event.path = os.path.join(event.rootpath, event.datapath, event.database, event.pylot_id)
+            event.path = os.path.join(event.datapath, event.pylot_id)
             event.path = event.path.replace('\\', '/')
             event.path = event.path.replace('//', '/')
 
@@ -1369,7 +1384,7 @@ class MainWindow(QMainWindow):
         index = eventBox.currentIndex()
         tv = QtWidgets.QTableView()
         header = tv.horizontalHeader()
-        header.setResizeMode(QtWidgets.QHeaderView.ResizeToContents)
+        header.setSectionResizeMode(QtWidgets.QHeaderView.ResizeToContents)
         header.setStretchLastSection(True)
         header.hide()
         tv.verticalHeader().hide()
@@ -1389,25 +1404,28 @@ class MainWindow(QMainWindow):
 
         for id, event in enumerate(self.project.eventlist):
             event_path = event.path
-            phaseErrors = {'P': self._inputs['timeerrorsP'],
+            #phaseErrors = {'P': self._inputs['timeerrorsP'],
-                           'S': self._inputs['timeerrorsS']}
+            #               'S': self._inputs['timeerrorsS']}
 
-            ma_props = {'manual': event.pylot_picks,
+            man_au_picks = {'manual': event.pylot_picks,
-                        'auto': event.pylot_autopicks}
+                            'auto': event.pylot_autopicks}
-            ma_count = {'manual': 0,
+            npicks = {'manual': {'P': 0, 'S': 0},
-                        'auto': 0}
+                      'auto': {'P': 0, 'S': 0}}
-            ma_count_total = {'manual': 0,
+            npicks_total = {'manual': {'P': 0, 'S': 0},
-                              'auto': 0}
+                            'auto': {'P': 0, 'S': 0}}
 
-            for ma in ma_props.keys():
+            for ma in man_au_picks.keys():
-                if ma_props[ma]:
+                if man_au_picks[ma]:
-                    for picks in ma_props[ma].values():
+                    for picks in man_au_picks[ma].values():
                         for phasename, pick in picks.items():
                             if not type(pick) in [dict, AttribDict]:
                                 continue
+                            phase_ID = identifyPhaseID(phasename)
+                            if not phase_ID in npicks[ma].keys():
+                                continue
                             if pick.get('spe'):
-                                ma_count[ma] += 1
+                                npicks[ma][phase_ID] += 1
-                            ma_count_total[ma] += 1
+                            npicks_total[ma][phase_ID] += 1
 
             event_ref = event.isRefEvent()
             event_test = event.isTestEvent()
@@ -1442,16 +1460,23 @@ class MainWindow(QMainWindow):
             if event.dirty:
                 event_str += '*'
             item_path = QStandardItem(event_str)
-            item_time = QStandardItem('{}'.format(time))
+            item_time = QStandardItem('{}'.format(time.strftime("%Y-%m-%d %H:%M:%S") if time else ''))
             item_lat = QStandardItem('{}'.format(lat))
             item_lon = QStandardItem('{}'.format(lon))
             item_depth = QStandardItem('{}'.format(depth))
             item_localmag = QStandardItem('{}'.format(localmag))
             item_momentmag = QStandardItem('{}'.format(momentmag))
-            item_nmp = QStandardItem('{}({})'.format(ma_count['manual'], ma_count_total['manual']))
+
+            item_nmp = QStandardItem()
+            item_nap = QStandardItem()
             item_nmp.setIcon(self.manupicksicon_small)
-            item_nap = QStandardItem('{}({})'.format(ma_count['auto'], ma_count_total['auto']))
             item_nap.setIcon(self.autopicksicon_small)
+
+            for picktype, item_np in [('manual', item_nmp), ('auto', item_nap)]:
+                npicks_str = f"{npicks[picktype]['P']}|{npicks[picktype]['S']}"
+                #npicks_str += f"({npicks_total[picktype]['P']}/{npicks_total[picktype]['S']})"
+                item_np.setText(npicks_str)
+
             item_ref = QStandardItem()  # str(event_ref))
             item_test = QStandardItem()  # str(event_test))
             if event_ref:
@@ -1512,8 +1537,8 @@ class MainWindow(QMainWindow):
             return True
         return not bool(os.listdir(wf_path))
 
-    def filename_from_action(self, action):
+    def filename_from_action(self, action=None):
-        if action.data() is None:
+        if not action or action.data() is None:
             filt = "Supported file formats" \
                    " (*.mat *.qml *.xml *.kor *.evt)"
             caption = "Open an event file"
@@ -1530,7 +1555,7 @@ class MainWindow(QMainWindow):
             self.set_fname(self.get_data().getEventFileName(), type)
         return self.get_fnames(type)
 
-    def saveData(self, event=None, directory=None, outformats=['.xml', '.cnv', '.obs', '_focmec.in', '.pha']):
+    def saveData(self, event=None, directory=None, outformats=None):
         '''
         Save event data to directory with specified output formats.
         :param event: PyLoT Event, if not set current event will be used
@@ -1538,6 +1563,8 @@ class MainWindow(QMainWindow):
         :param outformats: str/list of output formats
         :return:
         '''
+        if outformats is None:
+            outformats = ['.xml', '.cnv', '.obs', '_focmec.in', '.pha']
         if not event:
             event = self.get_current_event()
         if not type(outformats) == list:
@@ -1665,16 +1692,55 @@ class MainWindow(QMainWindow):
         self.cmpw.show()
 
     def pickQualities(self):
-        path = self._inputs['rootpath'] + '/' + self._inputs['datapath'] + '/' + self._inputs['database']
+        path = self.get_current_event_path()
-        getQualitiesfromxml(path)
+        (_, plot) = getQualitiesfromxml(path, self._inputs.get('timeerrorsP'), self._inputs.get('timeerrorsS'),plotflag=1)
+        self.pickQualitiesWindow = PickQualitiesFromXml(figure=plot, path=self.get_current_event_path(),inputVar=self._inputs)
+        self.pickQualitiesWindow.showUI()
         return 
 
+    # WIP JG
     def eventlistXml(self):
-        path = self._inputs['rootpath'] + '/' + self._inputs['datapath'] + '/' + self._inputs['database']
+        path = self._inputs['datapath']
         outpath = self.project.location[:self.project.location.rfind('/')]
         geteventlistfromxml(path, outpath)
         return
 
+        # WIP JG
+    def spectogramView(self):
+        global test
+        stations = []
+        names = []
+        traces = {}
+        for tr in self.get_data().wfdata.traces:
+            if not tr.stats.station in stations:
+                stations.append(tr.stats.station)
+                names.append(tr.stats.network + '.' + tr.stats.station)
+        for station in stations:
+            traces[station] = {}
+        for ch in ['Z', 'N', 'E']:
+            for tr in self.get_data().wfdata.select(component=ch).traces:
+                traces[tr.stats.station][ch] = tr
+
+        names.sort()
+        a = self.get_current_event()
+
+        print (self.get_data().wfdata.traces[0])
+
+        test = SpectrogramTab(traces, self.get_data().wfdata)
+        height = self.tabs.widget(0).height()
+        width = self.tabs.widget(0).width()
+        self.tabs.setCurrentIndex(3)
+        figCanvas = test.makeSpecFig(direction=self.dispComponent, height = height, width = width, parent = self.tabs.widget)
+        return figCanvas
+        #self.spectro_layout.addWidget()
+        # self.get_data().wfdata.spectrogram()
+        # self.tabs.addTab(figCanvas, 'Spectrogram')
+        # self.tabs[3] = figCanvas
+        # self.refreshTabs()
+        # test.show()
+
+
+
     def compareMulti(self):
         if not self.compareoptions:
             return
@@ -1833,6 +1899,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()
@@ -1860,6 +1927,12 @@ class MainWindow(QMainWindow):
                     self.newWF(plot=False)
                 self.update_obspy_dmt()
                 self.refresh_array_map()
+        if self.tabs.currentIndex() == 3:
+            if self.spectroWidget != None:
+                self.spectro_layout.removeWidget(self.spectroWidget)
+            newSpectroWidget = self.spectogramView()
+            self.spectro_layout.addWidget(newSpectroWidget)
+            self.spectroWidget = newSpectroWidget
 
     def newWF(self, event=None, plot=True):
         '''
@@ -1890,7 +1963,13 @@ class MainWindow(QMainWindow):
 
     def prepareLoadWaveformData(self):
         self.fnames = self.getWFFnames_from_eventbox()
-        self.fnames_syn = []
+        self.fnames_comp = []
+        fnames_comp = self.getWFFnames_from_eventbox(subpath='compare')
+        self.dataPlot.activateCompareOptions(bool(fnames_comp))
+        if fnames_comp:
+            if self.dataPlot.comp_checkbox.isChecked():
+                self.fnames_comp = fnames_comp
+
         eventpath = self.get_current_event_path()
         basepath = eventpath.split(os.path.basename(eventpath))[0]
         self.obspy_dmt = check_obspydmt_structure(basepath)
@@ -1911,6 +1990,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?')
@@ -1918,8 +1999,8 @@ class MainWindow(QMainWindow):
 
         if len(curr_event.origins) > 0:
             origin_time = curr_event.origins[0].time
-            tstart = settings.value('tstart') if get_None(settings.value('tstart')) else 0
+            tstart = settings.value('tstart') if get_none(settings.value('tstart')) else 0
-            tstop = settings.value('tstop') if get_None(settings.value('tstop')) else 0
+            tstop = settings.value('tstop') if get_none(settings.value('tstop')) else 0
             tstart = origin_time + float(tstart)
             tstop = origin_time + float(tstop)
         else:
@@ -1927,7 +2008,7 @@ class MainWindow(QMainWindow):
             tstop = None
 
         self.data.setWFData(self.fnames,
-                            self.fnames_syn,
+                            self.fnames_comp,
                             checkRotated=True,
                             metadata=self.metadata,
                             tstart=tstart,
@@ -1935,7 +2016,7 @@ class MainWindow(QMainWindow):
 
     def prepareObspyDMT_data(self, eventpath):
         qcbox_processed = self.dataPlot.qcombo_processed
-        qcheckb_syn = self.dataPlot.syn_checkbox
+        qcheckb_syn = self.dataPlot.comp_checkbox
         qcbox_processed.setEnabled(False)
         qcheckb_syn.setEnabled(False)
         for fpath in os.listdir(eventpath):
@@ -1943,8 +2024,8 @@ class MainWindow(QMainWindow):
             if 'syngine' in fpath:
                 eventpath_syn = os.path.join(eventpath, fpath)
                 qcheckb_syn.setEnabled(True)
-                if self.dataPlot.syn_checkbox.isChecked():
+                if self.dataPlot.comp_checkbox.isChecked():
-                    self.fnames_syn = [os.path.join(eventpath_syn, filename) for filename in os.listdir(eventpath_syn)]
+                    self.fnames_comp = [os.path.join(eventpath_syn, filename) for filename in os.listdir(eventpath_syn)]
             if 'processed' in fpath:
                 qcbox_processed.setEnabled(True)
         if qcbox_processed.isEnabled():
@@ -2066,15 +2147,16 @@ 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)
         for times, data, times_syn, data_syn in plots:
-            self.dataPlot.plotWidget.getPlotItem().plot(times, data,
+            self.dataPlot.plotWidget.getPlotItem().plot(np.array(times), np.array(data),
-                                                        pen=self.dataPlot.pen_linecolor)
+                                                        pen=self.dataPlot.pen_linecolor,
+                                                        skipFiniteCheck=True)
             if len(data_syn) > 0:
-                self.dataPlot.plotWidget.getPlotItem().plot(times_syn, data_syn,
+                self.dataPlot.plotWidget.getPlotItem().plot(np.array(times_syn), np.array(data_syn),
                                                             pen=self.dataPlot.pen_linecolor_syn)
         self.dataPlot.reinitMoveProxy()
         self.highlight_stations()
@@ -2114,9 +2196,11 @@ class MainWindow(QMainWindow):
         if event.pylot_autopicks:
             self.drawPicks(picktype='auto')
         if event.pylot_picks or event.pylot_autopicks:
-            self.locateEventAction.setEnabled(True)
+            if not self._inputs.get('extent') == 'global':
+                self.locateEventAction.setEnabled(True)
             self.qualities_action.setEnabled(True)
             self.eventlist_xml_action.setEnabled(True)
+
         if True in self.comparable.values():
             self.compare_action.setEnabled(True)
         self.draw()
@@ -2224,7 +2308,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()
+        wfsyn = self.get_data().getAltWFdata()
         if self.filterActionP.isChecked() and filter:
             self.filterWaveformData(plot=False, phase='P')
         elif self.filterActionS.isChecked() and filter:
@@ -2233,14 +2317,14 @@ class MainWindow(QMainWindow):
         # wfst += self.get_data().getWFData().select(component=alter_comp)
         plotWidget = self.getPlotWidget()
         self.adjustPlotHeight()
-        if get_Bool(settings.value('large_dataset')) == True:
+        if get_bool(settings.value('large_dataset')) == True:
             self.plot_method = 'fast'
         else:
             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 ([], [])
+        plots = rval if rval else []
-        return plots, gaps
+        return plots
 
     def adjustPlotHeight(self):
         if self.pg:
@@ -2343,7 +2427,7 @@ class MainWindow(QMainWindow):
         filterS = filteroptions['S']
         minP, maxP = filterP.getFreq()
         minS, maxS = filterS.getFreq()
-        self.paraBox.params_to_gui()
+        self.parameterWidget.params_to_gui()
 
     def getFilterOptions(self):
         return self.filteroptions
@@ -2536,17 +2620,20 @@ class MainWindow(QMainWindow):
             print("Warning! No network, station, and location info available!")
             return
         self.update_status('picking on station {0}'.format(station))
-        data = self.get_data().getOriginalWFData().copy()
+        wfdata = self.get_data().getOriginalWFData().copy()
+        wfdata_comp = self.get_data().getAltWFdata().copy()
         event = self.get_current_event()
         wftype = self.dataPlot.qcombo_processed.currentText() if self.obspy_dmt else None
         pickDlg = PickDlg(self, parameter=self._inputs,
-                          data=data.select(station=station),
+                          data=wfdata.select(station=station),
+                          data_compare=wfdata_comp.select(station=station),
                           station=station, network=network,
                           location=location,
                           picks=self.getPicksOnStation(station, 'manual'),
                           autopicks=self.getPicksOnStation(station, 'auto'),
                           metadata=self.metadata, event=event,
-                          filteroptions=self.filteroptions, wftype=wftype)
+                          filteroptions=self.filteroptions, wftype=wftype,
+                          show_comp_data=self.dataPlot.comp_checkbox.isChecked())
         if self.filterActionP.isChecked():
             pickDlg.currentPhase = "P"
             pickDlg.filterWFData()
@@ -2659,7 +2746,7 @@ class MainWindow(QMainWindow):
         self.init_fig_dict()
         # if not self.tap:
         # init TuneAutopicker object
-        self.tap = TuneAutopicker(self)
+        self.tap = TuneAutopicker(self, self.obspy_dmt)
         # first call of update to init tabs with empty canvas
         self.update_autopicker()
         # connect update signal of TuneAutopicker with update function
@@ -2858,7 +2945,9 @@ class MainWindow(QMainWindow):
         self.log_deleted_picks([deleted_pick])
 
     def log_deleted_picks(self, deleted_picks, event_path=None):
-        ''' Log deleted picks to list self.deleted_picks '''
+        '''
+        Log deleted picks to list self.deleted_picks
+        '''
         if not event_path:
             event_path = self.get_current_event_path()
         for deleted_pick in deleted_picks:
@@ -2872,7 +2961,9 @@ class MainWindow(QMainWindow):
             self.deleted_picks[event_path].append(deleted_pick)
 
     def dump_deleted_picks(self, event_path):
-        ''' Save deleted picks to json file for event in event_path. Load old file before and merge'''
+        '''
+        Save deleted picks to json file for event in event_path. Load old file before and merge
+        '''
         try:
             deleted_picks_from_file = self.load_deleted_picks(event_path)
         except Exception as e:
@@ -2921,10 +3012,16 @@ class MainWindow(QMainWindow):
             event = self.get_current_event()
         event.pylot_picks = {}
         event.pylot_autopicks = {}
-        picksdict = picksdict_from_picks(evt=self.get_data().get_evt_data())
+        picksdict = picksdict_from_picks(evt=self.get_data().get_evt_data(), parameter=self.getParameter())
         event.addPicks(picksdict['manual'])
         event.addAutopicks(picksdict['auto'])
 
+    def getParameter(self):
+        if hasattr(self.project, 'parameter') and isinstance(self.project.parameter, PylotParameter):
+            return self.project.parameter
+        else:
+            return self._inputs
+
     def drawPicks(self, station=None, picktype=None, stime=None):
         # if picktype not specified, draw both
         if not stime:
@@ -2999,7 +3096,7 @@ class MainWindow(QMainWindow):
 
                 if self.pg:
                     if spe:
-                        if picks['epp'] and picks['lpp']:
+                        if not self.plot_method == 'fast' and picks['epp'] and picks['lpp']:
                             pen = make_pen(picktype, phaseID, 'epp', quality)
                             self.drawnPicks[picktype][station].append(pw.plot([epp, epp], ylims,
                                                                               alpha=.25, pen=pen, name='EPP'))
@@ -3068,22 +3165,24 @@ class MainWindow(QMainWindow):
         lt = locateTool[loctool]
         # get working directory
         locroot = parameter['nllocroot']
+        #locroot = 'E:/NLL/src/Insheim'
         if locroot is None:
             self.PyLoTprefs()
             self.locate_event()
 
         ctrfile = os.path.join(locroot, 'run', parameter['ctrfile'])
+        #ctrfile = 'E:/NLL/src/Insheim/run/Insheim_min1d032016.in'
         ttt = parameter['ttpatter']
         outfile = parameter['outpatter']
         eventname = self.get_current_event_name()
-        obsdir = os.path.join(self._inputs['rootpath'], self._inputs['datapath'], self._inputs['database'], eventname)
+        obsdir = os.path.join(self._inputs['datapath'], eventname)
-        self.saveData(event=self.get_current_event(), directory=obsdir, outformats='.obs')
+        self.saveData(event=self.get_current_event(), directory=obsdir, outformats=['.obs'])
         filename = 'PyLoT_' + eventname
         locpath = os.path.join(locroot, 'loc', filename)
         phasefile = os.path.join(obsdir, filename + '.obs')
         lt.modify_inputs(ctrfile, locroot, filename, phasefile, ttt)
         try:
-            lt.locate(ctrfile)
+            lt.locate(ctrfile, self._inputs)
         except RuntimeError as e:
             print(e.message)
         # finally:
@@ -3137,7 +3236,7 @@ class MainWindow(QMainWindow):
         '''
         if checked: pass  # dummy argument for QAction trigger signal
         self.tabs.setCurrentIndex(1)
-        # if there is no metadata (invetories is an empty list), just initialize the default empty tab
+        # if there is no metadata (inventories is an empty list), just initialize the default empty tab
         if not self.metadata.inventories:
             self.init_array_tab()
             return
@@ -3159,6 +3258,15 @@ class MainWindow(QMainWindow):
         self.tabs.setCurrentIndex(index)
         self.refresh_array_map()
 
+    def init_spectro_tab(self):
+        '''
+        Init spectrogram tab with currently selected event.
+        '''
+        self.spectroWidget = None
+        #self.spectro_layout.addWidget( self.spectogramView() )
+        pass
+
+
     def array_map_thread(self):
         '''
         Start modal thread to init the array_map object.
@@ -3381,7 +3489,7 @@ class MainWindow(QMainWindow):
                     self.event_table.setCellWidget(r_index, c_index, item)
 
         header = self.event_table.horizontalHeader()
-        header.setResizeMode(QtWidgets.QHeaderView.ResizeToContents)
+        header.setSectionResizeMode(QtWidgets.QHeaderView.ResizeToContents)
         header.setStretchLastSection(True)
         self.event_table.cellChanged[int, int].connect(cell_changed)
         self.event_table.cellClicked[int, int].connect(cell_clicked)
@@ -3488,7 +3596,7 @@ class MainWindow(QMainWindow):
 
         wf_select = Stream()
         # restitute only picked traces
-        for station in np.unique(self.getPicks('manual').keys() + self.getPicks('auto').keys()):
+        for station in np.unique(list(self.getPicks('manual').keys()) + list(self.getPicks('auto').keys())):
             wf_select += wf_copy.select(station=station)
 
         corr_wf = restitute_data(wf_select, self.metadata)
@@ -3555,7 +3663,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:
@@ -3608,10 +3716,13 @@ class MainWindow(QMainWindow):
         if not self.okToContinue():
             return
         if not fnm:
-            dlg = QFileDialog(parent=self)
+            settings = QSettings()
+            dir = settings.value('current_project_path')
+            dlg = QFileDialog(parent=self, directory=dir)
             fnm = dlg.getOpenFileName(self, 'Open project file...', filter='Pylot project (*.plp)')[0]
             if not fnm:
                 return
+            settings.setValue('current_project_path', os.path.split(fnm)[0])
         if not os.path.exists(fnm):
             QMessageBox.warning(self, 'Could not open file',
                                 'Could not open project file {}. File does not exist.'.format(fnm))
@@ -3622,6 +3733,7 @@ class MainWindow(QMainWindow):
                 if self.project.parameter:
                     # do this step to update default parameter on older PyLoT projects
                     self.project.parameter.reinit_default_parameters()
+                    PylotParameter.check_deprecated_parameters(self.project.parameter)
 
                     self._inputs = self.project.parameter
                     self.updateFilteroptions()
@@ -3667,7 +3779,7 @@ class MainWindow(QMainWindow):
             filename = fnm[0] + '.plp'
         self.project.parameter = self._inputs
         settings = QSettings()
-        autosaveXML = get_Bool(settings.value('autosaveXML', True))
+        autosaveXML = get_bool(settings.value('autosaveXML', True))
         if autosaveXML:
             self.exportEvents()
         if not self.project.save(filename): return False
@@ -3691,7 +3803,7 @@ class MainWindow(QMainWindow):
                 self.metadata.clear_inventory()
                 self.project.parameter = self._inputs
                 settings = QSettings()
-                autosaveXML = get_Bool(settings.value('autosaveXML', True))
+                autosaveXML = get_bool(settings.value('autosaveXML', True))
                 if autosaveXML:
                     self.exportEvents()
                 if not self.project.save(): return False
@@ -3729,7 +3841,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()
@@ -3738,13 +3851,13 @@ class MainWindow(QMainWindow):
 
     def setParameter(self, checked=0, show=True):
         if checked: pass  # dummy argument to receive trigger signal (checked) if called by QAction
-        if not self.paraBox:
+        if not self.parameterWidget:
-            self.paraBox = PylotParaBox(self._inputs, parent=self, windowflag=Qt.Window)
+            self.parameterWidget = PylotParameterWidget(self._inputs, parent=self, windowflag=Qt.Window)
-            self.paraBox.accepted.connect(self._setDirty)
+            self.parameterWidget.accepted.connect(self._setDirty)
-            self.paraBox.accepted.connect(self.filterOptionsFromParameter)
+            self.parameterWidget.accepted.connect(self.filterOptionsFromParameter)
         if show:
-            self.paraBox.params_to_gui()
+            self.parameterWidget.params_to_gui()
-            self.paraBox.show()
+            self.parameterWidget.show()
 
     def deleteAllAutopicks(self):
         qmb = QMessageBox(self, icon=QMessageBox.Question,
@@ -3780,7 +3893,7 @@ class MainWindow(QMainWindow):
     def helpHelp(self):
         if checkurl():
             form = HelpForm(self,
-                            'https://ariadne.geophysik.ruhr-uni-bochum.de/trac/PyLoT/wiki')
+                            'https://github.com/seismology-RUB/PyLoT')
         else:
             form = HelpForm(self, ':/help.html')
         form.show()
@@ -3791,16 +3904,18 @@ class Project(object):
     Pickable class containing information of a PyLoT project, like event lists and file locations.
     '''
 
-    # TODO: remove rootpath
     def __init__(self):
         self.eventlist = []
         self.location = None
-        self.rootpath = None
         self.datapath = None
         self.dirty = False
         self.parameter = None
         self._table = None
 
+    @property
+    def rootpath(self):
+        return self.datapath
+
     def add_eventlist(self, eventlist):
         '''
         Add events from an eventlist containing paths to event directories.
@@ -3810,8 +3925,6 @@ class Project(object):
             return
         for item in eventlist:
             event = Event(item)
-            event.rootpath = self.parameter['rootpath']
-            event.database = self.parameter['database']
             event.datapath = self.parameter['datapath']
             if not event.path in self.getPaths():
                 self.eventlist.append(event)
@@ -3960,13 +4073,13 @@ class Project(object):
         return project
 
 
-class getExistingDirectories(QFileDialog):
+class GetExistingDirectories(QFileDialog):
     '''
     File dialog with possibility to select multiple folders.
     '''
 
     def __init__(self, *args):
-        super(getExistingDirectories, self).__init__(*args)
+        super(GetExistingDirectories, self).__init__(*args)
         self.setOption(self.DontUseNativeDialog, True)
         self.setOption(self.ReadOnly, True)
         self.setFileMode(self.Directory)
@@ -3992,16 +4105,7 @@ def create_window():
     return app, app_created
 
 
-def main(args=None):
+def main(project_filename=None, pylot_infile=None, reset_qsettings=False):
-    project_filename = None
-    # args.project_filename = 'C:/Shared/AlpArray/alparray_data/project_alparray_test.plp'
-    pylot_infile = None
-    if args:
-        if args.project_filename:
-            project_filename = args.project_filename
-        if args.input_filename:
-            pylot_infile = args.input_filename
-        reset_qsettings = args.reset_qsettings
 
     # create the Qt application
     pylot_app, app_created = create_window()
@@ -4050,4 +4154,5 @@ if __name__ == "__main__":
     parser.add_argument('--reset_qsettings', default=False, action='store_true',
                         help='reset qsettings (debug option)')
     args = parser.parse_args()
-    sys.exit(main(args))
+    sys.exit(main(project_filename=args.project_filename, pylot_infile=args.input_filename,
+                  reset_qsettings=args.reset_qsettings))

README.md

@@ -11,7 +11,7 @@ PILOT has originally been developed in Mathworks' MatLab. In order to distribute
 problems, it has been decided to redevelop the software package in Python. The great work of the ObsPy group allows easy
 handling of a bunch of seismic data and PyLoT will benefit a lot compared to the former MatLab version.
 
-The development of PyLoT is part of the joint research project MAGS2 and AlpArray.
+The development of PyLoT is part of the joint research project MAGS2, AlpArray and AdriaArray.
 
 ## Installation
 
@@ -24,61 +24,44 @@ It is highly recommended to use Anaconda for a simple creation of a Python insta
 
 Afterwards run (from the PyLoT main directory where the files *requirements.txt* and *pylot.yml* are located)
 
-    conda create --name pylot_38 --file requirements.txt
+    conda env create -f pylot.yml
 or
     
-    conda env create -f pylot.yml
+    conda create -c conda-forge --name pylot_311 python=3.11 --file requirements.txt
 
-to create a new Anaconda environment called "pylot_38".
+to create a new Anaconda environment called *pylot_311*.
 
 Afterwards activate the environment by typing
 
-    conda activate pylot_38
+    conda activate pylot_311
 
 #### Prerequisites:
 
 In order to run PyLoT you need to install:
 
 - Python 3
-- obspy
-- pyside2
-- pyqtgraph
 - cartopy
+- joblib
+- obspy
+- pyaml
+- pyqtgraph
+- pyside2
 
 (the following are already dependencies of the above packages):
 - scipy
 - numpy
-- matplotlib <= 3.3.x
+- matplotlib
 
 #### Some handwork:
 
-PyLoT needs a properties folder on your system to work. It should be situated in your home directory
+Some extra information on error estimates (just needed for reading old PILOT data) and the Richter magnitude scaling
-(on Windows usually C:/Users/*username*):
-
-    mkdir ~/.pylot
-
-In the next step you have to copy some files to this directory:
-
-*for local distance seismicity*
-
-    cp path-to-pylot/inputs/pylot_local.in ~/.pylot/pylot.in
-
-*for regional distance seismicity*
-
-    cp path-to-pylot/inputs/pylot_regional.in ~/.pylot/pylot.in
-
-*for global distance seismicity*
-
-    cp path-to-pylot/inputs/pylot_global.in ~/.pylot/pylot.in
-
-and some extra information on error estimates (just needed for reading old PILOT data) and the Richter magnitude scaling
 relation
 
     cp path-to-pylot/inputs/PILOT_TimeErrors.in path-to-pylot/inputs/richter_scaling.data ~/.pylot/
 
 You may need to do some modifications to these files. Especially folder names should be reviewed.
 
-PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10.
+PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10/11.
 
 ## Release notes
 
@@ -87,6 +70,7 @@ PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10.
 - event organisation in project files and waveform visualisation
 - consistent manual phase picking through predefined SNR dependant zoom level
 - consistent automatic phase picking routines using Higher Order Statistics, AIC and Autoregression
+- pick correlation correction for teleseismic waveforms
 - interactive tuning of auto-pick parameters
 - uniform uncertainty estimation from waveform's properties for automatic and manual picks
 - pdf representation and comparison of picks taking the uncertainty intrinsically into account
@@ -95,17 +79,17 @@ PyLoT has been tested on Mac OSX (10.11), Debian Linux 8 and on Windows 10.
 
 #### Known issues:
 
-We hope to solve these with the next release.
+Current release is still in development progress and has several issues. We are currently lacking manpower, but hope to assess many of the issues in the near future.
 
 ## Staff
 
-Original author(s): L. Kueperkoch, S. Wehling-Benatelli, M. Bischoff (PILOT)
+Original author(s): M. Rische, S. Wehling-Benatelli, L. Kueperkoch, M. Bischoff (PILOT)
 
-Developer(s): S. Wehling-Benatelli, L. Kueperkoch, K. Olbert, M. Bischoff, C. Wollin, M. Rische, M. Paffrath
+Developer(s): S. Wehling-Benatelli, M. Paffrath, L. Kueperkoch, K. Olbert, M. Bischoff, C. Wollin, M. Rische, D. Arnold, K. Cökerim, S. Zimmermann
 
 Others: A. Bruestle, T. Meier, W. Friederich
 
 
 [ObsPy]: http://github.com/obspy/obspy/wiki
 
-September 2017
+September 2024

autoPyLoT.py

@@ -28,7 +28,7 @@ from pylot.core.util.dataprocessing import restitute_data, Metadata
 from pylot.core.util.defaults import SEPARATOR
 from pylot.core.util.event import Event
 from pylot.core.util.structure import DATASTRUCTURE
-from pylot.core.util.utils import get_None, trim_station_components, check4gapsAndRemove, check4doubled, \
+from pylot.core.util.utils import get_none, trim_station_components, check4gapsAndRemove, check4doubled, \
     check4rotated
 from pylot.core.util.version import get_git_version as _getVersionString
 
@@ -91,9 +91,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
                                                               sp=sp_info)
     print(splash)
 
-    parameter = get_None(parameter)
+    parameter = get_none(parameter)
-    inputfile = get_None(inputfile)
+    inputfile = get_none(inputfile)
-    eventid = get_None(eventid)
+    eventid = get_none(eventid)
 
     fig_dict = None
     fig_dict_wadatijack = None
@@ -119,13 +119,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
             obspyDMT_wfpath = input_dict['obspyDMT_wfpath']
 
     if not parameter:
-        if inputfile:
+        if not inputfile:
-            parameter = PylotParameter(inputfile)
+            print('Using default input parameter')
-            # iplot = parameter['iplot']
+        parameter = PylotParameter(inputfile)
-        else:
-            infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
-            print('Using default input file {}'.format(infile))
-            parameter = PylotParameter(infile)
     else:
         if not type(parameter) == PylotParameter:
             print('Wrong input type for parameter: {}'.format(type(parameter)))
@@ -140,11 +136,9 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
     if parameter.hasParam('datastructure'):
         # getting information on data structure
         datastructure = DATASTRUCTURE[parameter.get('datastructure')]()
-        dsfields = {'root': parameter.get('rootpath'),
+        dsfields = {'dpath': parameter.get('datapath'),}
-                    'dpath': parameter.get('datapath'),
-                    'dbase': parameter.get('database')}
 
-        exf = ['root', 'dpath', 'dbase']
+        exf = ['dpath']
 
         if parameter['eventID'] != '*' and fnames == 'None':
             dsfields['eventID'] = parameter['eventID']
@@ -154,7 +148,7 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
         datastructure.setExpandFields(exf)
 
         # check if default location routine NLLoc is available and all stations are used
-        if get_None(parameter['nllocbin']) and station == 'all':
+        if get_none(parameter['nllocbin']) and station == 'all':
             locflag = 1
             # get NLLoc-root path
             nllocroot = parameter.get('nllocroot')
@@ -190,15 +184,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
         if not input_dict:
             # started in production mode
             datapath = datastructure.expandDataPath()
-            if fnames == 'None' and parameter['eventID'] == '*':
+            if fnames in [None, 'None'] and parameter['eventID'] == '*':
                 # multiple event processing
                 # read each event in database
                 events = [event for event in glob.glob(os.path.join(datapath, '*')) if
                           (os.path.isdir(event) and not event.endswith('EVENTS-INFO'))]
-            elif fnames == 'None' and parameter['eventID'] != '*' and not type(parameter['eventID']) == list:
+            elif fnames in [None, 'None'] and parameter['eventID'] != '*' and not type(parameter['eventID']) == list:
                 # single event processing
                 events = glob.glob(os.path.join(datapath, parameter['eventID']))
-            elif fnames == 'None' and type(parameter['eventID']) == list:
+            elif fnames in [None, 'None'] and type(parameter['eventID']) == list:
                 # multiple event processing
                 events = []
                 for eventID in parameter['eventID']:
@@ -210,12 +204,10 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
                 locflag = 2
         else:
             # started in tune or interactive mode
-            datapath = os.path.join(parameter['rootpath'],
+            datapath = parameter['datapath']
-                                    parameter['datapath'])
             events = []
             for eventID in eventid:
                 events.append(os.path.join(datapath,
-                                           parameter['database'],
                                            eventID))
 
         if not events:
@@ -242,12 +234,15 @@ def autoPyLoT(input_dict=None, parameter=None, inputfile=None, fnames=None, even
                 data.get_evt_data().path = eventpath
                 print('Reading event data from filename {}...'.format(filename))
             except Exception as e:
-                print('Could not read event from file {}: {}'.format(filename, e))
+                if type(e) == FileNotFoundError:
+                    print('Creating new event file.')
+                else:
+                    print('Could not read event from file {}: {}'.format(filename, e))
                 data = Data()
                 pylot_event = Event(eventpath)  # event should be path to event directory
                 data.setEvtData(pylot_event)
-            if fnames == 'None':
+            if fnames in [None, 'None']:
-                data.setWFData(glob.glob(os.path.join(datapath, event_datapath, '*')))
+                data.setWFData(glob.glob(os.path.join(event_datapath, '*')))
                 # the following is necessary because within
                 # multiple event processing no event ID is provided
                 # in autopylot.in

autopylot.sh

@@ -3,8 +3,10 @@
 #$ -l low
 #$ -cwd
 #$ -pe smp 40
-#$ -l mem=2G
+##$ -l mem=3G
-#$ -l h_vmem=2G
+#$ -l h_vmem=6G
 #$ -l os=*stretch
 
-python ./autoPyLoT.py -i /home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in -dmt processed -c $NSLOTS
+conda activate pylot_311
+
+python ./autoPyLoT.py -i /home/marcel/.pylot/pylot_adriaarray.in -c 20 -dmt processed

docs/gui.md

@@ -203,8 +203,6 @@ The meaning of the header entries is:
 PyLoT GUI starts with an empty project. To add events, use the add event data button. Select one or multiple folders
 containing events.
 
-[//]: <> (TODO: explain _Directories: Root path, Data path, Database path_)
-
 ### Saving projects
 
 Save the current project from the menu with File->Save project or File->Save project as. PyLoT uses ``.plp`` files to

inputs/pylot_global.in

@@ -4,10 +4,8 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-                                                   #rootpath#       %project path                                                                                        
                                                    #datapath#       %data path
-                                                   #database#       %name of data base                                                                                   
+                                                   #eventID#        %event ID for single event processing (* for all events found in datapath)
-                                                   #eventID#        %event ID for single event processing (* for all events found in database)                           
                                                    #invdir#         %full path to  inventory or dataless-seed file                                                       
 PILOT                                              #datastructure#  %choose data structure                                                                               
 True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
@@ -43,6 +41,7 @@ global                                             #extent#         %extent of a
 1150.0                                             #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
 True                                               #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
 iasp91                                             #taup_model#     %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
+P,Pdiff                                            #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
 0.05 0.5                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
 0.001 0.5                                          #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
 0.05 0.5                                           #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     

inputs/pylot_local.in

@@ -4,10 +4,8 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-/DATA/Insheim                                      #rootpath#       %project path                                                                                        
+/DATA/Insheim/EVENT_DATA/LOCAL/2018.02_Insheim     #datapath#       %data path
-EVENT_DATA/LOCAL                                   #datapath#       %data path                                                                                           
+e0006.038.18                                       #eventID#        %event ID for single event processing (* for all events found in datapath)
-2018.02_Insheim                                    #database#       %name of data base                                                                                   
-e0006.038.18                                       #eventID#        %event ID for single event processing (* for all events found in database)                           
 /DATA/Insheim/STAT_INFO                            #invdir#         %full path to  inventory or dataless-seed file                                                       
 PILOT                                              #datastructure#  %choose data structure                                                                               
 True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
@@ -43,6 +41,7 @@ local                                              #extent#         %extent of a
 10.0                                               #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
 False                                              #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
 iasp91                                             #taup_model#     %define TauPy model for traveltime estimation
+P                                                  #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
 2.0 20.0                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
 2.0 30.0                                           #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
 2.0 10.0                                           #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     

inputs/pylot_regional.in

@@ -4,10 +4,8 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-                                                   #rootpath#       %project path                                                                                        
                                                    #datapath#       %data path
-                                                   #database#       %name of data base                                                                                   
+                                                   #eventID#        %event ID for single event processing (* for all events found in datapath)
-                                                   #eventID#        %event ID for single event processing (* for all events found in database)                           
                                                    #invdir#         %full path to  inventory or dataless-seed file                                                       
 PILOT                                              #datastructure#  %choose data structure                                                                               
 True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
@@ -43,6 +41,7 @@ local                                              #extent#         %extent of a
 10.0                                               #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
 True                                               #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
 iasp91                                             #taup_model#     %define TauPy model for traveltime estimation
+P                                                  #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
 2.0 10.0                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
 2.0 12.0                                           #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
 2.0 8.0                                            #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     

pylot.yml

@@ -1,14 +1,12 @@
-name: pylot_38
+name: pylot_311
 channels:
   - conda-forge
   - defaults
 dependencies:
-  - cartopy=0.20.2
+  - cartopy=0.23.0=py311hcf9f919_1
-  - matplotlib-base=3.3.4
+  - joblib=1.4.2=pyhd8ed1ab_0
-  - numpy=1.22.3
+  - obspy=1.4.1=py311he736701_3
-  - obspy=1.3.0
+  - pyaml=24.7.0=pyhd8ed1ab_0
-  - pyqtgraph=0.12.4
+  - pyqtgraph=0.13.7=pyhd8ed1ab_0
-  - pyside2=5.13.2
+  - pyside2=5.15.8=py311h3d699ce_4
-  - python=3.8.12
+  - pytest=8.3.2=pyhd8ed1ab_0
-  - qt=5.12.9
-  - scipy=1.8.0

pylot/__init__.py

@@ -9,7 +9,7 @@ PyLoT - the Python picking and Localization Tool
 
 This python library contains a graphical user interfaces for picking
 seismic phases. This software needs ObsPy (http://github.com/obspy/obspy/wiki)
-and the Qt4 libraries to be installed first.
+and the Qt libraries to be installed first.
 
 PILOT has been developed in Mathworks' MatLab. In order to distribute
 PILOT without facing portability problems, it has been decided to re-

pylot/core/analysis/magnitude.py

@@ -144,6 +144,10 @@ class Magnitude(object):
                     azimuthal_gap=self.origin_id.get_referred_object().quality.azimuthal_gap)
             else:
                 # no scaling necessary
+                # Temporary fix needs rework
+                if (len(self.magnitudes.keys()) == 0):
+                    print("Error in local magnitude calculation ")
+                    return None
                 mag = ope.Magnitude(
                     mag=np.median([M.mag for M in self.magnitudes.values()]),
                     magnitude_type=self.type,
@@ -414,6 +418,10 @@ class MomentMagnitude(Magnitude):
             distance = degrees2kilometers(a.distance)
             azimuth = a.azimuth
             incidence = a.takeoff_angle
+            if not 0. <= incidence <= 360.:
+                if self.verbose:
+                    print(f'WARNING: Incidence angle outside bounds - {incidence}')
+                return
             w0, fc = calcsourcespec(scopy, onset, self.p_velocity, distance,
                                     azimuth, incidence, self.p_attenuation,
                                     self.plot_flag, self.verbose)
@@ -432,6 +440,40 @@ class MomentMagnitude(Magnitude):
             self.event.station_magnitudes.append(magnitude)
             self.magnitudes = (station, magnitude)
 
+    # WIP JG
+    def getSourceSpec(self):
+        for a in self.arrivals:
+            if a.phase not in 'pP':
+                continue
+            # make sure calculating Mo only from reliable onsets
+            # NLLoc: time_weight = 0 => do not use onset!
+            if a.time_weight == 0:
+                continue
+            pick = a.pick_id.get_referred_object()
+            station = pick.waveform_id.station_code
+            if len(self.stream) <= 2:
+                print("Station:" '{0}'.format(station))
+                print("WARNING: No instrument corrected data available,"
+                      " no magnitude calculation possible! Go on.")
+                continue
+            wf = self.stream.select(station=station)
+            if not wf:
+                continue
+            try:
+                scopy = wf.copy()
+            except AssertionError:
+                print("WARNING: Something's wrong with the data,"
+                      "station {},"
+                      "no calculation of moment magnitude possible! Go on.".format(station))
+                continue
+            onset = pick.time
+            distance = degrees2kilometers(a.distance)
+            azimuth = a.azimuth
+            incidence = a.takeoff_angle
+            w0, fc, plt = calcsourcespec(scopy, onset, self.p_velocity, distance,
+                                    azimuth, incidence, self.p_attenuation,
+                                    3, self.verbose)
+            return w0, fc, plt
 
 def calcMoMw(wfstream, w0, rho, vp, delta, verbosity=False):
     '''
@@ -602,7 +644,7 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
         N = min(int(np.power(m, 2)), 16384)
         # N = int(np.power(m, 2))
         y = dt * np.fft.fft(xdat, N)
-        Y = abs(y[: N / 2])
+        Y = abs(y[: int(N / 2)])
         L = (N - 1) / freq
         f = np.arange(0, fny, 1 / L)
 
@@ -679,6 +721,8 @@ def calcsourcespec(wfstream, onset, vp, delta, azimuth, incidence,
             plt.xlabel('Frequency [Hz]')
             plt.ylabel('Amplitude [m/Hz]')
             plt.grid()
+            if iplot == 3:
+                return w0, Fc, plt
             plt.show()
             try:
                 input()

pylot/core/io/data.py

@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 
 import copy
+import logging
 import os
 
 from PySide2.QtWidgets import QMessageBox
@@ -19,7 +20,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 +65,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:
@@ -260,7 +261,6 @@ class Data(object):
         can be a str or a list of strings of ['manual', 'auto', 'origin', 'magnitude']
         """
         from pylot.core.util.defaults import OUTPUTFORMATS
-
         if not type(fcheck) == list:
             fcheck = [fcheck]
 
@@ -321,35 +321,61 @@ class Data(object):
             if lendiff != 0:
                 print("Manual as well as automatic picks available. Prefered the {} manual ones!".format(lendiff))
 
+
+            no_uncertainties_p = []
+            no_uncertainties_s = []
             if upperErrors:
                 # check for pick uncertainties exceeding adjusted upper errors
                 # Picks with larger uncertainties will not be saved in output file!
                 for j in range(len(picks)):
                     for i in range(len(picks_copy)):
                         if picks_copy[i].phase_hint[0] == 'P':
-                            if (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[0]) or \
+                            # Skipping pick if no upper_uncertainty is found and warning user
-                                    (picks_copy[i].time_errors['uncertainty'] is None):
+                            if picks_copy[i].time_errors['upper_uncertainty'] is None:
+                                #print("{1} P-Pick of station {0} does not have upper_uncertainty and cant be checked".format(
+                                 #   picks_copy[i].waveform_id.station_code,
+                                  #  picks_copy[i].method_id))
+                                if not picks_copy[i].waveform_id.station_code in no_uncertainties_p:
+                                    no_uncertainties_p.append(picks_copy[i].waveform_id.station_code)
+                                continue
+
+                            #print ("checking for upper_uncertainty")
+                            if  (picks_copy[i].time_errors['uncertainty'] is None) or \
+                                    (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[0]):
                                 print("Uncertainty exceeds or equal adjusted upper time error!")
                                 print("Adjusted uncertainty: {}".format(upperErrors[0]))
                                 print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
                                 print("{1} P-Pick of station {0} will not be saved in outputfile".format(
                                     picks_copy[i].waveform_id.station_code,
                                     picks_copy[i].method_id))
-                                print("#")
                                 del picks_copy[i]
                                 break
                         if picks_copy[i].phase_hint[0] == 'S':
-                            if (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[1]) or \
+
-                                    (picks_copy[i].time_errors['uncertainty'] is None):
+                            # Skipping pick if no upper_uncertainty is found and warning user
+                            if picks_copy[i].time_errors['upper_uncertainty'] is None:
+                                #print("{1} S-Pick of station {0} does not have upper_uncertainty and cant be checked".format(
+                                    #picks_copy[i].waveform_id.station_code,
+                                    #picks_copy[i].method_id))
+                                if not picks_copy[i].waveform_id.station_code in no_uncertainties_s:
+                                    no_uncertainties_s.append(picks_copy[i].waveform_id.station_code)
+                                continue
+
+
+                            if (picks_copy[i].time_errors['uncertainty'] is None) or \
+                                (picks_copy[i].time_errors['upper_uncertainty'] >= upperErrors[1]):
                                 print("Uncertainty exceeds or equal adjusted upper time error!")
                                 print("Adjusted uncertainty: {}".format(upperErrors[1]))
                                 print("Pick uncertainty: {}".format(picks_copy[i].time_errors['uncertainty']))
                                 print("{1} S-Pick of station {0} will not be saved in outputfile".format(
                                     picks_copy[i].waveform_id.station_code,
                                     picks_copy[i].method_id))
-                                print("#")
                                 del picks_copy[i]
                                 break
+            for s in no_uncertainties_p:
+                print("P-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s))
+            for s in no_uncertainties_s:
+                print("S-Pick of station {0} does not have upper_uncertainty and cant be checked".format(s))
 
             if fnext == '.obs':
                 try:
@@ -383,18 +409,16 @@ class Data(object):
                                      not implemented: {1}'''.format(evtformat, e))
             if fnext == '_focmec.in':
                 try:
-                    infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
+                    parameter = PylotParameter()
-                    print('Using default input file {}'.format(infile))
+                    logging.warning('Using default input parameter')
-                    parameter = PylotParameter(infile)
                     focmec.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
                 except KeyError as e:
                     raise KeyError('''{0} export format
                                      not implemented: {1}'''.format(evtformat, e))
             if fnext == '.pha':
                 try:
-                    infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
+                    parameter = PylotParameter()
-                    print('Using default input file {}'.format(infile))
+                    logging.warning('Using default input parameter')
-                    parameter = PylotParameter(infile)
                     hypodd.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
                 except KeyError as e:
                     raise KeyError('''{0} export format
@@ -426,20 +450,30 @@ class Data(object):
         data.filter(**kwargs)
         self.dirty = True
 
-    def setWFData(self, fnames, fnames_syn=None, checkRotated=False, metadata=None, tstart=0, tstop=0):
+    def setWFData(self, fnames, fnames_alt=None, checkRotated=False, metadata=None, tstart=0, tstop=0):
         """
         Clear current waveform data and set given waveform data
         :param fnames: waveform data names to append
+        :param fnames_alt: alternative data to show (e.g. synthetic/processed)
         :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()
+        self.wf_alt = Stream()
         if tstart == tstop:
             tstart = tstop = None
         self.tstart = tstart
         self.tstop = tstop
 
+        # remove directories
+        fnames = check_fname_exists(fnames)
+        fnames_alt = check_fname_exists(fnames_alt)
+
         # if obspy_dmt:
         #     wfdir = 'raw'
         #     self.processed = False
@@ -457,8 +491,8 @@ class Data(object):
         #     wffnames = fnames
         if fnames is not None:
             self.appendWFData(fnames)
-            if fnames_syn is not None:
+            if fnames_alt is not None:
-                self.appendWFData(fnames_syn, synthetic=True)
+                self.appendWFData(fnames_alt, alternative=True)
         else:
             return False
 
@@ -466,7 +500,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)
@@ -478,7 +514,7 @@ class Data(object):
         self.dirty = False
         return True
 
-    def appendWFData(self, fnames, synthetic=False):
+    def appendWFData(self, fnames, alternative=False):
         """
         Read waveform data from fnames and append it to current wf data
         :param fnames: waveform data to append
@@ -491,20 +527,20 @@ class Data(object):
         if self.dirty:
             self.resetWFData()
 
-        real_or_syn_data = {True: self.wfsyn,
+        orig_or_alternative_data = {True: self.wf_alt,
-                            False: self.wfdata}
+                                    False: self.wfdata}
 
         warnmsg = ''
         for fname in set(fnames):
             try:
-                real_or_syn_data[synthetic] += read(fname, starttime=self.tstart, endtime=self.tstop)
+                orig_or_alternative_data[alternative] += read(fname, starttime=self.tstart, endtime=self.tstop)
             except TypeError:
                 try:
-                    real_or_syn_data[synthetic] += read(fname, format='GSE2', starttime=self.tstart, endtime=self.tstop)
+                    orig_or_alternative_data[alternative] += read(fname, format='GSE2', starttime=self.tstart, endtime=self.tstop)
                 except Exception as e:
                     try:
-                        real_or_syn_data[synthetic] += read(fname, format='SEGY', starttime=self.tstart,
+                        orig_or_alternative_data[alternative] += read(fname, format='SEGY', starttime=self.tstart,
-                                                            endtime=self.tstop)
+                                                                      endtime=self.tstop)
                     except Exception as e:
                         warnmsg += '{0}\n{1}\n'.format(fname, e)
             except SacIOError as se:
@@ -519,8 +555,8 @@ class Data(object):
     def getOriginalWFData(self):
         return self.wforiginal
 
-    def getSynWFData(self):
+    def getAltWFdata(self):
-        return self.wfsyn
+        return self.wf_alt
 
     def resetWFData(self):
         """

pylot/core/io/default_parameters.py

@@ -6,24 +6,14 @@ import numpy as np
 Default parameters used for picking
 """
 
-defaults = {'rootpath': {'type': str,
+defaults = {'datapath': {'type': str,
-                         'tooltip': 'project path',
+                         'tooltip': 'path to eventfolders',
-                         'value': '',
-                         'namestring': 'Root path'},
-
-            'datapath': {'type': str,
-                         'tooltip': 'data path',
                          'value': '',
                          'namestring': 'Data path'},
 
-            'database': {'type': str,
-                         'tooltip': 'name of data base',
-                         'value': '',
-                         'namestring': 'Database path'},
-
             'eventID': {'type': str,
-                        'tooltip': 'event ID for single event processing (* for all events found in database)',
+                        'tooltip': 'event ID for single event processing (* for all events found in datapath)',
-                        'value': '',
+                        'value': '*',
                         'namestring': 'Event ID'},
 
             'extent': {'type': str,
@@ -511,7 +501,7 @@ defaults = {'rootpath': {'type': str,
 
             'taup_model': {'type': str,
                            'tooltip': 'Define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6',
-                           'value': None,
+                           'value': 'iasp91',
                            'namestring': 'TauPy model'},
 
             'taup_phases': {'type': str,
@@ -522,9 +512,7 @@ defaults = {'rootpath': {'type': str,
 
 settings_main = {
     'dirs': [
-        'rootpath',
         'datapath',
-        'database',
         'eventID',
         'invdir',
         'datastructure',

pylot/core/io/getEventListFromXML.py

@@ -8,6 +8,11 @@
    Edited for use in PyLoT
    JG, igem, 01/2022
 """
+
+import os
+import argparse
+import numpy as np
+import matplotlib.pyplot as plt
 import glob
 
 from obspy.core.event import read_events

pylot/core/io/getQualitiesfromxml.py

@@ -1,138 +0,0 @@
-#!/usr/bin/python
-# -*- coding: utf-8 -*-
-
-"""
-   Script to get onset uncertainties from Quakeml.xml files created by PyLoT.
-   Uncertainties are tranformed into quality classes and visualized via histogram if desired.
-   Ludger Küperkoch, BESTEC GmbH, 07/2017
-   rev.: Ludger Küperkoch, igem, 10/2020
-   Edited for usage in PyLoT: Jeldrik Gaal, igem, 01/2022
-"""
-
-import glob
-
-import matplotlib.pyplot as plt
-import numpy as np
-from obspy.core.event import read_events
-
-
-def getQualitiesfromxml(path):
-    # uncertainties
-    ErrorsP = [0.02, 0.04, 0.08, 0.16]
-    ErrorsS = [0.04, 0.08, 0.16, 0.32]
-
-    Pw0 = []
-    Pw1 = []
-    Pw2 = []
-    Pw3 = []
-    Pw4 = []
-    Sw0 = []
-    Sw1 = []
-    Sw2 = []
-    Sw3 = []
-    Sw4 = []
-
-    # data path
-    dp = path + '/e*/*.xml'
-    # list of all available xml-files
-    xmlnames = glob.glob(dp)
-
-    # read all onset weights
-    for names in xmlnames:
-        print("Getting onset weights from {}".format(names))
-        cat = read_events(names)
-        arrivals = cat.events[0].picks
-        for Pick in arrivals:
-            if Pick.phase_hint[0] == 'P':
-                if Pick.time_errors.uncertainty <= ErrorsP[0]:
-                    Pw0.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsP[0] and \
-                        Pick.time_errors.uncertainty <= ErrorsP[1]:
-                    Pw1.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsP[1] and \
-                        Pick.time_errors.uncertainty <= ErrorsP[2]:
-                    Pw2.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsP[2] and \
-                        Pick.time_errors.uncertainty <= ErrorsP[3]:
-                    Pw3.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsP[3]:
-                    Pw4.append(Pick.time_errors.uncertainty)
-                else:
-                    pass
-            elif Pick.phase_hint[0] == 'S':
-                if Pick.time_errors.uncertainty <= ErrorsS[0]:
-                    Sw0.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsS[0] and \
-                        Pick.time_errors.uncertainty <= ErrorsS[1]:
-                    Sw1.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsS[1] and \
-                        Pick.time_errors.uncertainty <= ErrorsS[2]:
-                    Sw2.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsS[2] and \
-                        Pick.time_errors.uncertainty <= ErrorsS[3]:
-                    Sw3.append(Pick.time_errors.uncertainty)
-                elif Pick.time_errors.uncertainty > ErrorsS[3]:
-                    Sw4.append(Pick.time_errors.uncertainty)
-                else:
-                    pass
-            else:
-                print("Phase hint not defined for picking!")
-                pass
-    # get percentage of weights
-    numPweights = np.sum([len(Pw0), len(Pw1), len(Pw2), len(Pw3), len(Pw4)])
-    numSweights = np.sum([len(Sw0), len(Sw1), len(Sw2), len(Sw3), len(Sw4)])
-    try:
-        P0perc = 100.0 / numPweights * len(Pw0)
-    except:
-        P0perc = 0
-    try:
-        P1perc = 100.0 / numPweights * len(Pw1)
-    except:
-        P1perc = 0
-    try:
-        P2perc = 100.0 / numPweights * len(Pw2)
-    except:
-        P2perc = 0
-    try:
-        P3perc = 100.0 / numPweights * len(Pw3)
-    except:
-        P3perc = 0
-    try:
-        P4perc = 100.0 / numPweights * len(Pw4)
-    except:
-        P4perc = 0
-    try:
-        S0perc = 100.0 / numSweights * len(Sw0)
-    except:
-        Soperc = 0
-    try:
-        S1perc = 100.0 / numSweights * len(Sw1)
-    except:
-        S1perc = 0
-    try:
-        S2perc = 100.0 / numSweights * len(Sw2)
-    except:
-        S2perc = 0
-    try:
-        S3perc = 100.0 / numSweights * len(Sw3)
-    except:
-        S3perc = 0
-    try:
-        S4perc = 100.0 / numSweights * len(Sw4)
-    except:
-        S4perc = 0
-
-    weights = ('0', '1', '2', '3', '4')
-    y_pos = np.arange(len(weights))
-    width = 0.34
-    p1, = plt.bar(0 - width, P0perc, width, color='black')
-    p2, = plt.bar(0, S0perc, width, color='red')
-    plt.bar(y_pos - width, [P0perc, P1perc, P2perc, P3perc, P4perc], width, color='black')
-    plt.bar(y_pos, [S0perc, S1perc, S2perc, S3perc, S4perc], width, color='red')
-    plt.ylabel('%')
-    plt.xticks(y_pos, weights)
-    plt.xlim([-0.5, 4.5])
-    plt.xlabel('Qualities')
-    plt.title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
-    plt.legend([p1, p2], ['P-Weights', 'S-Weights'])
-    plt.show()

pylot/core/io/inputs.py

@@ -1,5 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
+import logging
+import os
 
 from pylot.core.io import default_parameters
 from pylot.core.util.errors import ParameterError
@@ -88,10 +90,10 @@ class PylotParameter(object):
         return bool(self.__parameter)
 
     def __getitem__(self, key):
-        try:
+        if key in self.__parameter:
             return self.__parameter[key]
-        except:
+        else:
-            return None
+            logging.warning(f'{key} not found in PylotParameter')
 
     def __setitem__(self, key, value):
         try:
@@ -418,6 +420,28 @@ class PylotParameter(object):
             line = value + name + ttip
             fid.write(line)
 
+    @staticmethod
+    def check_deprecated_parameters(parameters):
+        if parameters.hasParam('database') and parameters.hasParam('rootpath'):
+            parameters['datapath'] = os.path.join(parameters['rootpath'], parameters['datapath'],
+                                                  parameters['database'])
+            logging.warning(
+                f'Parameters database and rootpath are deprecated. '
+                f'Tried to merge them to now path: {parameters["datapath"]}.'
+            )
+
+        remove_keys = []
+        for key in parameters:
+            if not key in default_parameters.defaults.keys():
+                remove_keys.append(key)
+                logging.warning(f'Removing deprecated parameter: {key}')
+
+        for key in remove_keys:
+            del parameters[key]
+
+        parameters._settings_main = default_parameters.settings_main
+        parameters._settings_special_pick = default_parameters.settings_special_pick
+
 
 class FilterOptions(object):
     '''

pylot/core/io/location.py

@@ -1,7 +1,7 @@
 from obspy import UTCDateTime
 from obspy.core import event as ope
 
-from pylot.core.util.utils import getLogin, getHash
+from pylot.core.util.utils import get_login, get_hash
 
 
 def create_amplitude(pickID, amp, unit, category, cinfo):
@@ -61,7 +61,7 @@ def create_creation_info(agency_id=None, creation_time=None, author=None):
     :return:
     '''
     if author is None:
-        author = getLogin()
+        author = get_login()
     if creation_time is None:
         creation_time = UTCDateTime()
     return ope.CreationInfo(agency_id=agency_id, author=author,
@@ -210,7 +210,7 @@ def create_resourceID(timetohash, restype, authority_id=None, hrstr=None):
     '''
     assert isinstance(timetohash, UTCDateTime), "'timetohash' is not an ObsPy" \
                                                 "UTCDateTime object"
-    hid = getHash(timetohash)
+    hid = get_hash(timetohash)
     if hrstr is None:
         resID = ope.ResourceIdentifier(restype + '/' + hid[0:6])
     else:

pylot/core/io/phases.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 import glob
+import logging
 import os
 import warnings
 
@@ -16,8 +17,8 @@ from pylot.core.io.inputs import PylotParameter
 from pylot.core.io.location import create_event, \
     create_magnitude
 from pylot.core.pick.utils import select_for_phase, get_quality_class
-from pylot.core.util.utils import getOwner, full_range, four_digits, transformFilterString4Export, \
+from pylot.core.util.utils import get_owner, full_range, four_digits, transformFilterString4Export, \
-    backtransformFilterString
+    backtransformFilterString, loopIdentifyPhase, identifyPhase
 
 
 def add_amplitudes(event, amplitudes):
@@ -58,7 +59,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
     if phasfn is not None and os.path.isfile(phasfn):
         phases = sio.loadmat(phasfn)
         phasctime = UTCDateTime(os.path.getmtime(phasfn))
-        phasauthor = getOwner(phasfn)
+        phasauthor = get_owner(phasfn)
     else:
         phases = None
         phasctime = None
@@ -66,7 +67,7 @@ def readPILOTEvent(phasfn=None, locfn=None, authority_id='RUB', **kwargs):
     if locfn is not None and os.path.isfile(locfn):
         loc = sio.loadmat(locfn)
         locctime = UTCDateTime(os.path.getmtime(locfn))
-        locauthor = getOwner(locfn)
+        locauthor = get_owner(locfn)
     else:
         loc = None
         locctime = None
@@ -217,12 +218,14 @@ def picksdict_from_obs(fn):
     return picks
 
 
-def picksdict_from_picks(evt):
+def picksdict_from_picks(evt, parameter=None, nwst_id: bool = False):
     """
     Takes an Event object and return the pick dictionary commonly used within
     PyLoT
     :param evt: Event object contain all available information
     :type evt: `~obspy.core.event.Event`
+    :param nwst_id: determines if network and station id are used or only station id
+    :type nwst_id: bool
     :return: pick dictionary (auto and manual)
     """
     picksdict = {
@@ -230,8 +233,12 @@ def picksdict_from_picks(evt):
         'auto': {}
     }
     for pick in evt.picks:
+        errors = None
         phase = {}
-        station = pick.waveform_id.station_code
+        if not nwst_id:
+            station_or_nwst = pick.waveform_id.station_code
+        else:
+            station_or_nwst = f'{pick.waveform_id.station_code}.{pick.waveform_id.network_code}'
         if pick.waveform_id.channel_code is None:
             channel = ''
         else:
@@ -252,7 +259,7 @@ def picksdict_from_picks(evt):
         if pick_method == 'None':
             pick_method = 'manual'
         try:
-            onsets = picksdict[pick_method][station]
+            onsets = picksdict[pick_method][station_or_nwst]
         except KeyError as e:
             # print(e)
             onsets = {}
@@ -273,36 +280,33 @@ def picksdict_from_picks(evt):
         phase['epp'] = epp
         phase['lpp'] = lpp
         phase['spe'] = spe
-        try:
+        weight = phase.get('weight')
-            phase['weight'] = weight
+        if not weight:
-        except:
+            if not parameter:
-            # get onset weight from uncertainty
+                logging.warning('Using ')
-            infile = os.path.join(os.path.expanduser('~'), '.pylot', 'pylot.in')
+                logging.warning('Using default input parameter')
-            print('Using default input file {}'.format(infile))
+                parameter = PylotParameter()
-            parameter = PylotParameter(infile)
+            pick.phase_hint = identifyPhase(pick.phase_hint)
             if pick.phase_hint == 'P':
                 errors = parameter['timeerrorsP']
             elif pick.phase_hint == 'S':
                 errors = parameter['timeerrorsS']
-            weight = get_quality_class(spe, errors)
+            if errors:
-            phase['weight'] = weight
+                weight = get_quality_class(spe, errors)
+                phase['weight'] = weight
         phase['channel'] = channel
         phase['network'] = network
         phase['picker'] = pick_method
-        try:
+        if pick.polarity == 'positive':
-            if pick.polarity == 'positive':
+            phase['fm'] = 'U'
-                phase['fm'] = 'U'
+        elif pick.polarity == 'negative':
-            elif pick.polarity == 'negative':
+            phase['fm'] = 'D'
-                phase['fm'] = 'D'
+        else:
-            else:
-                phase['fm'] = 'N'
-        except:
-            print("No FM info available!")
             phase['fm'] = 'N'
         phase['filter_id'] = filter_id if filter_id is not None else ''
 
         onsets[pick.phase_hint] = phase.copy()
-        picksdict[pick_method][station] = onsets.copy()
+        picksdict[pick_method][station_or_nwst] = onsets.copy()
     return picksdict
 
 
@@ -375,7 +379,6 @@ def picks_from_picksdict(picks, creation_info=None):
 
 
 def reassess_pilot_db(root_dir, db_dir, out_dir=None, fn_param=None, verbosity=0):
-    import glob
     # TODO: change root to datapath
     db_root = os.path.join(root_dir, db_dir)
     evt_list = glob.glob1(db_root, 'e????.???.??')
@@ -514,11 +517,11 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         fid = open("%s" % filename, 'w')
         # write header
         fid.write('# EQEVENT: %s Label: EQ%s  Loc:  X 0.00  Y 0.00  Z 10.00  OT 0.00 \n' %
-                  (parameter.get('database'), parameter.get('eventID')))
+                  (parameter.get('datapath'), parameter.get('eventID')))
         arrivals = chooseArrivals(arrivals)  # MP MP what is chooseArrivals? It is not defined anywhere
         for key in arrivals:
             # P onsets
-            if arrivals[key].has_key('P'):
+            if 'P' in arrivals[key]:
                 try:
                     fm = arrivals[key]['P']['fm']
                 except KeyError as e:
@@ -552,7 +555,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
                                                                                            ss_ms,
                                                                                            pweight))
             # S onsets
-            if arrivals[key].has_key('S') and arrivals[key]['S']['mpp'] is not None:
+            if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None:
                 fm = '?'
                 onset = arrivals[key]['S']['mpp']
                 year = onset.year
@@ -667,11 +670,11 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         print("Writing phases to %s for HYPOSAT" % filename)
         fid = open("%s" % filename, 'w')
         # write header
-        fid.write('%s, event %s \n' % (parameter.get('database'), parameter.get('eventID')))
+        fid.write('%s, event %s \n' % (parameter.get('datapath'), parameter.get('eventID')))
         arrivals = chooseArrivals(arrivals)  # MP MP what is chooseArrivals? It is not defined anywhere
         for key in arrivals:
             # P onsets
-            if arrivals[key].has_key('P') and arrivals[key]['P']['mpp'] is not None:
+            if 'P' in  arrivals[key] and arrivals[key]['P']['mpp'] is not None:
                 if arrivals[key]['P']['weight'] < 4:
                     Ponset = arrivals[key]['P']['mpp']
                     pyear = Ponset.year
@@ -700,7 +703,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
                     fid.write('%-5s P1       %4.0f %02d %02d %02d %02d %05.02f   %5.3f -999.   0.00 -999.  0.00\n'
                               % (key, pyear, pmonth, pday, phh, pmm, Pss, pstd))
             # S onsets
-            if arrivals[key].has_key('S') and arrivals[key]['S']['mpp'] is not None:
+            if 'S' in arrivals[key] and arrivals[key]['S']['mpp'] is not None:
                 if arrivals[key]['S']['weight'] < 4:
                     Sonset = arrivals[key]['S']['mpp']
                     syear = Sonset.year
@@ -769,7 +772,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         usedarrivals = chooseArrivals(arrivals)
         for key in usedarrivals:
             # P onsets
-            if usedarrivals[key].has_key('P'):
+            if 'P' in usedarrivals[key]:
                 if usedarrivals[key]['P']['weight'] < 4:
                     n += 1
                     stat = key
@@ -783,7 +786,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
                     else:
                         fid.write('%-4sP%d%6.2f\n' % (stat, Pweight, Prt))
                         # S onsets
-            if usedarrivals[key].has_key('S'):
+            if 'S' in usedarrivals[key]:
                 if usedarrivals[key]['S']['weight'] < 4:
                     n += 1
                     stat = key
@@ -829,13 +832,13 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         # prefer manual picks
         usedarrivals = chooseArrivals(arrivals)
         for key in usedarrivals:
-            if usedarrivals[key].has_key('P'):
+            if 'P' in usedarrivals[key]:
                 # P onsets
                 if usedarrivals[key]['P']['weight'] < 4:
                     Ponset = usedarrivals[key]['P']['mpp']
                     Prt = Ponset - stime  # onset time relative to source time
                     fid.write('%s    %6.3f  1  P\n' % (key, Prt))
-            if usedarrivals[key].has_key('S'):
+            if 'S' in usedarrivals[key]:
                 # S onsets
                 if usedarrivals[key]['S']['weight'] < 4:
                     Sonset = usedarrivals[key]['S']['mpp']
@@ -880,7 +883,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         # prefer manual picks
         usedarrivals = chooseArrivals(arrivals)
         for key in usedarrivals:
-            if usedarrivals[key].has_key('P'):
+            if 'P' in usedarrivals[key]:
                 if usedarrivals[key]['P']['weight'] < 4 and usedarrivals[key]['P']['fm'] is not None:
                     stat = key
                     for i in range(len(picks)):
@@ -962,7 +965,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         arrivals = chooseArrivals(arrivals)  # MP MP what is chooseArrivals? It is not defined anywhere
         # write phase lines
         for key in arrivals:
-            if arrivals[key].has_key('P'):
+            if 'P' in arrivals[key]:
                 if arrivals[key]['P']['weight'] < 4 and arrivals[key]['P']['fm'] is not None:
                     stat = key
                     ccode = arrivals[key]['P']['channel']
@@ -1015,7 +1018,6 @@ def chooseArrivals(arrivals):
     :return: arrivals but with the manual picks prefered if possible
     """
     # If len of arrivals is greater than 2 it comes from autopicking so only autopicks are available
-    print("=== CHOOSE ===")
     if len(arrivals) > 2:
         return arrivals
     if arrivals['auto'] and arrivals['manual']:
@@ -1056,37 +1058,63 @@ def merge_picks(event, picks):
     return event
 
 
-def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
+def getQualitiesfromxml(path, errorsP, errorsS, plotflag=1, figure=None, verbosity=0):
     """
     Script to get onset uncertainties from Quakeml.xml files created by PyLoT.
     Uncertainties are tranformed into quality classes and visualized via histogram if desired.
     Ludger Küperkoch, BESTEC GmbH, 07/2017
-    :param xmlnames: list of xml obspy event files containing picks
+    :param path: path containing xml files
-    :type xmlnames: list
+    :type path: str
-    :param ErrorsP: time errors of P waves for the four discrete quality classes
+    :param errorsP: time errors of P waves for the four discrete quality classes
-    :type ErrorsP:
+    :type errorsP:
-    :param ErrorsS: time errors of S waves for the four discrete quality classes
+    :param errorsS: time errors of S waves for the four discrete quality classes
-    :type ErrorsS:
+    :type errorsS:
     :param plotflag:
     :type plotflag:
     :return:
     :rtype:
     """
 
-    from pylot.core.pick.utils import get_quality_class
+    def calc_perc(uncertainties, ntotal):
-    from pylot.core.util.utils import loopIdentifyPhase, identifyPhase
+        ''' simple function that calculates percentage of number of uncertainties (list length)'''
+        if len(uncertainties) == 0:
+            return 0
+        else:
+            return 100. / ntotal * len(uncertainties)
+
+    def calc_weight_perc(psweights, weight_ids):
+        ''' calculate percentages of different weights (pick classes!?) of total number of uncertainties of a phase'''
+        # count total number of list items for this phase
+        numWeights = np.sum([len(weight) for weight in psweights.values()])
+
+        # iterate over all available weights to return a list with percentages for plotting
+        plot_list = []
+        for weight_id in weight_ids:
+            plot_list.append(calc_perc(psweights[weight_id], numWeights))
+
+        return plot_list, numWeights
+
+    # get all xmlfiles in path (maybe this should be changed to one xml file for this function, selectable via GUI?)
+    xmlnames = glob.glob(os.path.join(path, '*.xml'))
+    if len(xmlnames) == 0:
+        print(f'No files found in path {path}.')
+        return False
+
+    # first define possible phases here
+    phases = ['P', 'S']
+
+    # define possible weights (0-4)
+    weight_ids = list(range(5))
+
+    # put both error lists in a dictionary with P/S key so that amount of code can be halfed by simply using P/S as key
+    errors = dict(P=errorsP, S=errorsS)
+
+    # create dictionaries for each phase (P/S) with a dictionary of empty list for each weight defined in weights
+    # tuple above
+    weights = {}
+    for phase in phases:
+        weights[phase] = {weight_id: [] for weight_id in weight_ids}
 
-    # read all onset weights
-    Pw0 = []
-    Pw1 = []
-    Pw2 = []
-    Pw3 = []
-    Pw4 = []
-    Sw0 = []
-    Sw1 = []
-    Sw2 = []
-    Sw3 = []
-    Sw4 = []
     for names in xmlnames:
         print("Getting onset weights from {}".format(names))
         cat = read_events(names)
@@ -1094,117 +1122,60 @@ def getQualitiesfromxml(xmlnames, ErrorsP, ErrorsS, plotflag=1):
         arrivals = cat.events[0].picks
         arrivals_copy = cat_copy.events[0].picks
         # Prefere manual picks if qualities are sufficient!
-        for Pick in arrivals:
+        for pick in arrivals:
-            if Pick.method_id.id.split('/')[1] == 'manual':
+            if pick.method_id.id.split('/')[1] == 'manual':
-                mstation = Pick.waveform_id.station_code
+                mstation = pick.waveform_id.station_code
                 mstation_ext = mstation + '_'
                 for mpick in arrivals_copy:
-                    phase = identifyPhase(loopIdentifyPhase(Pick.phase_hint))
+                    phase = identifyPhase(loopIdentifyPhase(pick.phase_hint)) # MP MP catch if this fails?
-                    if phase == 'P':
+                    if ((mpick.waveform_id.station_code == mstation) or
-                        if ((mpick.waveform_id.station_code == mstation) or
+                        (mpick.waveform_id.station_code == mstation_ext)) and \
-                            (mpick.waveform_id.station_code == mstation_ext)) and \
+                            (mpick.method_id.id.split('/')[1] == 'auto') and \
-                                (mpick.method_id.split('/')[1] == 'auto') and \
+                            (mpick.time_errors['uncertainty'] <= errors[phase][3]):
-                                (mpick.time_errors['uncertainty'] <= ErrorsP[3]):
+                        del mpick
-                            del mpick
+                        break
-                            break
-                    elif phase == 'S':
-                        if ((mpick.waveform_id.station_code == mstation) or
-                            (mpick.waveform_id.station_code == mstation_ext)) and \
-                                (mpick.method_id.split('/')[1] == 'auto') and \
-                                (mpick.time_errors['uncertainty'] <= ErrorsS[3]):
-                            del mpick
-                            break
         lendiff = len(arrivals) - len(arrivals_copy)
         if lendiff != 0:
             print("Found manual as well as automatic picks, prefered the {} manual ones!".format(lendiff))
 
-        for Pick in arrivals_copy:
+        for pick in arrivals_copy:
-            phase = identifyPhase(loopIdentifyPhase(Pick.phase_hint))
+            phase = identifyPhase(loopIdentifyPhase(pick.phase_hint))
-            if phase == 'P':
+            uncertainty = pick.time_errors.uncertainty
-                Pqual = get_quality_class(Pick.time_errors.uncertainty, ErrorsP)
+            if not uncertainty:
-                if Pqual == 0:
+                if verbosity > 0:
-                    Pw0.append(Pick.time_errors.uncertainty)
+                    print('No uncertainty, pick {} invalid!'.format(pick.method_id.id))
-                elif Pqual == 1:
+                continue
-                    Pw1.append(Pick.time_errors.uncertainty)
+            # check P/S phase
-                elif Pqual == 2:
+            if phase not in phases:
-                    Pw2.append(Pick.time_errors.uncertainty)
-                elif Pqual == 3:
-                    Pw3.append(Pick.time_errors.uncertainty)
-                elif Pqual == 4:
-                    Pw4.append(Pick.time_errors.uncertainty)
-            elif phase == 'S':
-                Squal = get_quality_class(Pick.time_errors.uncertainty, ErrorsS)
-                if Squal == 0:
-                    Sw0.append(Pick.time_errors.uncertainty)
-                elif Squal == 1:
-                    Sw1.append(Pick.time_errors.uncertainty)
-                elif Squal == 2:
-                    Sw2.append(Pick.time_errors.uncertainty)
-                elif Squal == 3:
-                    Sw3.append(Pick.time_errors.uncertainty)
-                elif Squal == 4:
-                    Sw4.append(Pick.time_errors.uncertainty)
-            else:
                 print("Phase hint not defined for picking!")
-                pass
+                continue
+
+            qual = get_quality_class(uncertainty, errors[phase])
+            weights[phase][qual].append(uncertainty)
 
     if plotflag == 0:
-        Punc = [Pw0, Pw1, Pw2, Pw3, Pw4]
+        p_unc = [weights['P'][weight_id] for weight_id in weight_ids]
-        Sunc = [Sw0, Sw1, Sw2, Sw3, Sw4]
+        s_unc = [weights['S'][weight_id] for weight_id in weight_ids]
-        return Punc, Sunc
+        return p_unc, s_unc
     else:
+        if not figure:
+            fig = plt.figure()
+        ax = fig.add_subplot(111)
         # get percentage of weights
-        numPweights = np.sum([len(Pw0), len(Pw1), len(Pw2), len(Pw3), len(Pw4)])
+        listP, numPweights = calc_weight_perc(weights['P'], weight_ids)
-        numSweights = np.sum([len(Sw0), len(Sw1), len(Sw2), len(Sw3), len(Sw4)])
+        listS, numSweights = calc_weight_perc(weights['S'], weight_ids)
-        if len(Pw0) > 0:
-            P0perc = 100 / numPweights * len(Pw0)
-        else:
-            P0perc = 0
-        if len(Pw1) > 0:
-            P1perc = 100 / numPweights * len(Pw1)
-        else:
-            P1perc = 0
-        if len(Pw2) > 0:
-            P2perc = 100 / numPweights * len(Pw2)
-        else:
-            P2perc = 0
-        if len(Pw3) > 0:
-            P3perc = 100 / numPweights * len(Pw3)
-        else:
-            P3perc = 0
-        if len(Pw4) > 0:
-            P4perc = 100 / numPweights * len(Pw4)
-        else:
-            P4perc = 0
-        if len(Sw0) > 0:
-            S0perc = 100 / numSweights * len(Sw0)
-        else:
-            S0perc = 0
-        if len(Sw1) > 0:
-            S1perc = 100 / numSweights * len(Sw1)
-        else:
-            S1perc = 0
-        if len(Sw2) > 0:
-            S2perc = 100 / numSweights * len(Sw2)
-        else:
-            S2perc = 0
-        if len(Sw3) > 0:
-            S3perc = 100 / numSweights * len(Sw3)
-        else:
-            S3perc = 0
-        if len(Sw4) > 0:
-            S4perc = 100 / numSweights * len(Sw4)
-        else:
-            S4perc = 0
 
-        weights = ('0', '1', '2', '3', '4')
+        y_pos = np.arange(len(weight_ids))
-        y_pos = np.arange(len(weights))
         width = 0.34
-        plt.bar(y_pos - width, [P0perc, P1perc, P2perc, P3perc, P4perc], width, color='black')
+        ax.bar(y_pos - width, listP, width, color='black')
-        plt.bar(y_pos, [S0perc, S1perc, S2perc, S3perc, S4perc], width, color='red')
+        ax.bar(y_pos, listS, width, color='red')
-        plt.ylabel('%')
+        ax.set_ylabel('%')
-        plt.xticks(y_pos, weights)
+        ax.set_xticks(y_pos, weight_ids)
-        plt.xlim([-0.5, 4.5])
+        ax.set_xlim([-0.5, 4.5])
-        plt.xlabel('Qualities')
+        ax.set_xlabel('Qualities')
-        plt.title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
+        ax.set_title('{0} P-Qualities, {1} S-Qualities'.format(numPweights, numSweights))
-        plt.show()
+
+        if not figure:
+            fig.show()
+
+        return listP, listS

pylot/core/loc/nll.py

@@ -82,8 +82,8 @@ def locate(fnin, parameter=None):
     :param fnin: external program name
     :return: None
     """
-    exe_path = which('NLLoc', parameter)
+    exe_path = os.path.join(parameter['nllocbin'], 'NLLoc')
-    if exe_path is None:
+    if not os.path.isfile(exe_path):
         raise NLLocError('NonLinLoc executable not found; check your '
                          'environment variables')
 

pylot/core/pick/autopick.py

@@ -20,9 +20,9 @@ from pylot.core.pick.charfuns import CharacteristicFunction
 from pylot.core.pick.charfuns import HOScf, AICcf, ARZcf, ARHcf, AR3Ccf
 from pylot.core.pick.picker import AICPicker, PragPicker
 from pylot.core.pick.utils import checksignallength, checkZ4S, earllatepicker, \
-    getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class
+    getSNR, fmpicker, checkPonsets, wadaticheck, get_quality_class, PickingFailedException, MissingTraceException
 from pylot.core.util.utils import getPatternLine, gen_Pool, \
-    get_Bool, identifyPhaseID, get_None, correct_iplot
+    get_bool, identifyPhaseID, get_none, correct_iplot
 
 
 def autopickevent(data, param, iplot=0, fig_dict=None, fig_dict_wadatijack=None, ncores=0, metadata=None, origin=None):
@@ -232,20 +232,6 @@ class PickingContainer:
         self.Sflag = 0
 
 
-class MissingTraceException(ValueError):
-    """
-    Used to indicate missing traces in a obspy.core.stream.Stream object
-    """
-    pass
-
-
-class PickingFailedException(Exception):
-    """
-    Raised when picking fails due to missing values etc.
-    """
-    pass
-
-
 class AutopickStation(object):
 
     def __init__(self, wfstream, pickparam, verbose, iplot=0, fig_dict=None, metadata=None, origin=None):
@@ -272,10 +258,14 @@ class AutopickStation(object):
         self.pickparams = copy.deepcopy(pickparam)
         self.verbose = verbose
         self.iplot = correct_iplot(iplot)
-        self.fig_dict = get_None(fig_dict)
+        self.fig_dict = get_none(fig_dict)
         self.metadata = metadata
         self.origin = origin
 
+        # initialize TauPy pick estimates
+        self.estFirstP = None
+        self.estFirstS = None
+
         # initialize picking results
         self.p_results = PickingResults()
         self.s_results = PickingResults()
@@ -429,7 +419,7 @@ class AutopickStation(object):
             if station_coords is None:
                 exit_taupy()
                 raise AttributeError('Warning: Could not find station in metadata')
-            # TODO raise when metadata.get_coordinates returns None
+            # TODO: raise when metadata.get_coordinates returns None
             source_origin = origin[0]
             model = TauPyModel(taup_model)
             taup_phases = self.pickparams['taup_phases']
@@ -457,25 +447,27 @@ class AutopickStation(object):
             for arr in arrivals:
                 phases[identifyPhaseID(arr.phase.name)].append(arr)
             # get first P and S onsets from arrivals list
-            estFirstP = 0
+            arrival_time_p = 0
-            estFirstS = 0
+            arrival_time_s = 0
             if len(phases['P']) > 0:
-                arrP, estFirstP = min([(arr, arr.time) for arr in phases['P']], key=lambda t: t[1])
+                arrP, arrival_time_p = min([(arr, arr.time) for arr in phases['P']], key=lambda t: t[1])
             if len(phases['S']) > 0:
-                arrS, estFirstS = min([(arr, arr.time) for arr in phases['S']], key=lambda t: t[1])
+                arrS, arrival_time_s = min([(arr, arr.time) for arr in phases['S']], key=lambda t: t[1])
             print('autopick: estimated first arrivals for P: {} s, S:{} s after event'
-                  ' origin time using TauPy'.format(estFirstP, estFirstS))
+                  ' origin time using TauPy'.format(arrival_time_p, arrival_time_s))
-            return estFirstP, estFirstS
+            return arrival_time_p, arrival_time_s
 
         def exit_taupy():
             """If taupy failed to calculate theoretical starttimes, picking continues.
             For this a clean exit is required, since the P starttime is no longer relative to the theoretic onset but
             to the vertical trace starttime, eg. it can't be < 0."""
+            # TODO here the pickparams is modified, instead of a copy
             if self.pickparams["pstart"] < 0:
-                # TODO here the pickparams is modified, instead of a copy
                 self.pickparams["pstart"] = 0
+            if self.pickparams["sstart"] < 0:
+                self.pickparams["sstart"] = 0
 
-        if self.pickparams["use_taup"] is False:
+        if get_bool(self.pickparams["use_taup"]) is False:
             # correct user mistake where a relative cuttime is selected (pstart < 0) but use of taupy is disabled/ has
             # not the required parameters
             exit_taupy()
@@ -489,16 +481,34 @@ class AutopickStation(object):
             raise AttributeError('No source origins given!')
 
         arrivals = create_arrivals(self.metadata, self.origin, self.pickparams["taup_model"])
-        estFirstP, estFirstS = first_PS_onsets(arrivals)
+        arrival_P, arrival_S = first_PS_onsets(arrivals)
+
+        self.estFirstP = (self.origin[0].time + arrival_P) - self.ztrace.stats.starttime
+
         # modifiy pstart and pstop relative to estimated first P arrival (relative to station time axis)
-        self.pickparams["pstart"] += (self.origin[0].time + estFirstP) - self.ztrace.stats.starttime
+        self.pickparams["pstart"] += self.estFirstP
-        self.pickparams["pstop"] += (self.origin[0].time + estFirstP) - self.ztrace.stats.starttime
+        self.pickparams["pstop"] += self.estFirstP
         print('autopick: CF calculation times respectively:'
               ' pstart: {} s, pstop: {} s'.format(self.pickparams["pstart"], self.pickparams["pstop"]))
         # make sure pstart and pstop are inside the starttime/endtime of vertical trace
         self.pickparams["pstart"] = max(self.pickparams["pstart"], 0)
         self.pickparams["pstop"] = min(self.pickparams["pstop"], len(self.ztrace) * self.ztrace.stats.delta)
 
+        if self.horizontal_traces_exist():
+            # for the two horizontal components take earliest and latest time to make sure that the s onset is not clipped
+            # if start and endtime of horizontal traces differ, the s windowsize will automatically increase
+            trace_s_start = min([self.etrace.stats.starttime, self.ntrace.stats.starttime])
+            self.estFirstS = (self.origin[0].time + arrival_S) - trace_s_start
+            # modifiy sstart and sstop relative to estimated first S arrival (relative to station time axis)
+            self.pickparams["sstart"] += self.estFirstS
+            self.pickparams["sstop"] += self.estFirstS
+            print('autopick: CF calculation times respectively:'
+                  ' sstart: {} s, sstop: {} s'.format(self.pickparams["sstart"], self.pickparams["sstop"]))
+            # make sure pstart and pstop are inside the starttime/endtime of horizontal traces
+            self.pickparams["sstart"] = max(self.pickparams["sstart"], 0)
+            self.pickparams["sstop"] = min(self.pickparams["sstop"], len(self.ntrace) * self.ntrace.stats.delta,
+                                           len(self.etrace) * self.etrace.stats.delta)
+
     def autopickstation(self):
         """
         Main function of autopickstation, which calculates P and S picks and returns them in a dictionary.
@@ -508,6 +518,17 @@ class AutopickStation(object):
         station's value is the station name on which the picks were calculated.
         :rtype: dict
         """
+
+        if get_bool(self.pickparams['use_taup']) is True and self.origin is not None:
+            try:
+                # modify pstart, pstop, sstart, sstop to be around theoretical onset if taupy should be used,
+                # else do nothing
+                self.modify_starttimes_taupy()
+            except AttributeError as ae:
+                print(ae)
+            except MissingTraceException as mte:
+                print(mte)
+
         try:
             self.pick_p_phase()
         except MissingTraceException as mte:
@@ -515,13 +536,15 @@ class AutopickStation(object):
         except PickingFailedException as pfe:
             print(pfe)
 
-        if self.horizontal_traces_exist() and self.p_results.weight is not None and self.p_results.weight < 4:
+        if self.horizontal_traces_exist():
-            try:
+            if (self.p_results.weight is not None and self.p_results.weight < 4) or \
-                self.pick_s_phase()
+                    get_bool(self.pickparams.get('use_taup')):
-            except MissingTraceException as mte:
+                try:
-                print(mte)
+                    self.pick_s_phase()
-            except PickingFailedException as pfe:
+                except MissingTraceException as mte:
-                print(pfe)
+                    print(mte)
+                except PickingFailedException as pfe:
+                    print(pfe)
 
         self.plot_pick_results()
         self.finish_picking()
@@ -594,6 +617,12 @@ class AutopickStation(object):
             # plot tapered trace filtered with bpz2 filter settings
             ax1.plot(tdata, self.tr_filt_z_bpz2.data / max(self.tr_filt_z_bpz2.data), color=linecolor, linewidth=0.7,
                      label='Data')
+            # plot pickwindows for P
+            pstart, pstop = self.pickparams['pstart'], self.pickparams['pstop']
+            if pstart is not None and pstop is not None:
+                ax1.axvspan(pstart, pstop, color='r', alpha=0.1, zorder=0, label='P window')
+            if self.estFirstP is not None:
+                ax1.axvline(self.estFirstP, ls='dashed', color='r', alpha=0.4, label='TauPy estimate')
             if self.p_results.weight < 4:
                 # plot CF of initial onset (HOScf or ARZcf)
                 ax1.plot(self.cf1.getTimeArray(), self.cf1.getCF() / max(self.cf1.getCF()), 'b', label='CF1')
@@ -631,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,
@@ -661,13 +690,14 @@ class AutopickStation(object):
                         ax2.plot([self.s_results.lpp, self.s_results.lpp], [-1.1, 1.1], 'g--', label='lpp')
                         ax2.plot([self.s_results.epp, self.s_results.epp], [-1.1, 1.1], 'g--', label='epp')
                         title = '{channel}, S weight={sweight}, SNR={snr:7.2}, SNR[dB]={snrdb:7.2}'
-                        ax2.set_title(title.format(channel=self.etrace.stats.channel,
+                        ax2.set_title(title.format(channel=str(self.etrace.stats.channel),
-                                                   sweight=self.s_results.weight,
+                                                   sweight=str(self.s_results.weight),
-                                                   snr=self.s_results.snr,
+                                                   snr=str(self.s_results.snr),
-                                                   snrdb=self.s_results.snrdb))
+                                                   snrdb=str(self.s_results.snrdb)))
                     else:
                         title = '{channel}, S weight={sweight}, SNR=None, SNR[dB]=None'
-                        ax2.set_title(title.format(channel=self.etrace.stats.channel, sweight=self.s_results.weight))
+                        ax2.set_title(title.format(channel=str(self.etrace.stats.channel),
+                                                   sweight=str(self.s_results.weight)))
                 ax2.legend(loc=1)
                 ax2.set_yticks([])
                 ax2.set_ylim([-1.5, 1.5])
@@ -697,6 +727,15 @@ class AutopickStation(object):
                         ax3.plot([refSpick.getpick() - 0.5, refSpick.getpick() + 0.5], [-1.3, -1.3], 'g', linewidth=2)
                         ax3.plot([self.s_results.lpp, self.s_results.lpp], [-1.1, 1.1], 'g--', label='lpp')
                         ax3.plot([self.s_results.epp, self.s_results.epp], [-1.1, 1.1], 'g--', label='epp')
+
+                # plot pickwindows for S
+                sstart, sstop = self.pickparams['sstart'], self.pickparams['sstop']
+                if sstart is not None and sstop is not None:
+                    for axis in [ax2, ax3]:
+                        axis.axvspan(sstart, sstop, color='b', alpha=0.1, zorder=0, label='S window')
+                        if self.estFirstS is not None:
+                            axis.axvline(self.estFirstS, ls='dashed', color='b', alpha=0.4, label='TauPy estimate')
+
                 ax3.legend(loc=1)
                 ax3.set_yticks([])
                 ax3.set_ylim([-1.5, 1.5])
@@ -795,15 +834,6 @@ class AutopickStation(object):
         # save filtered trace in instance for later plotting
         self.tr_filt_z_bpz2 = tr_filt
 
-        if get_Bool(self.pickparams['use_taup']) is True and self.origin is not None:
-            try:
-                # modify pstart, pstop to be around theoretical onset if taupy should be used, else does nothing
-                self.modify_starttimes_taupy()
-            except AttributeError as ae:
-                print(ae)
-            except MissingTraceException as mte:
-                print(mte)
-
         Lc = self.pickparams['pstop'] - self.pickparams['pstart']
 
         Lwf = self.ztrace.stats.endtime - self.ztrace.stats.starttime
@@ -828,14 +858,21 @@ class AutopickStation(object):
             self.cf1 = None
         assert isinstance(self.cf1, CharacteristicFunction), 'cf1 is not set correctly: maybe the algorithm name ({})' \
                                                              ' is corrupted'.format(self.pickparams["algoP"])
+        # get the original waveform stream from first CF class cut to identical length as CF for plotting
+        cut_ogstream = self.cf1.getDataArray(self.cf1.getCut())
+
+        # MP: Rename to cf_stream for further use of z_copy and to prevent chaos when z_copy suddenly becomes a cf
+        # stream and later again a waveform stream
+        cf_stream = z_copy.copy()
+        cf_stream[0].data = self.cf1.getCF()
+
         # calculate AIC cf from first cf (either HOS or ARZ)
-        z_copy[0].data = self.cf1.getCF()
+        aiccf = AICcf(cf_stream, cuttimes)
-        aiccf = AICcf(z_copy, cuttimes)
         # get preliminary onset time from AIC-CF
         self.set_current_figure('aicFig')
         aicpick = AICPicker(aiccf, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
                             Tsmooth=self.pickparams["aictsmooth"], fig=self.current_figure,
-                            linecolor=self.current_linecolor)
+                            linecolor=self.current_linecolor, ogstream=cut_ogstream)
         # save aicpick for plotting later
         self.p_data.aicpick = aicpick
         # add pstart and pstop to aic plot
@@ -848,7 +885,7 @@ class AutopickStation(object):
                           label='P stop')
                 ax.legend(loc=1)
 
-        Pflag = self._pick_p_quality_control(aicpick, z_copy, tr_filt)
+        Pflag = self._pick_p_quality_control(aicpick, cf_stream, tr_filt)
         # go on with processing if AIC onset passes quality control
         slope = aicpick.getSlope()
         if not slope: slope = 0
@@ -887,7 +924,7 @@ class AutopickStation(object):
         refPpick = PragPicker(self.cf2, self.pickparams["tsnrz"], self.pickparams["pickwinP"], self.iplot,
                               self.pickparams["ausP"],
                               self.pickparams["tsmoothP"], aicpick.getpick(), self.current_figure,
-                              self.current_linecolor)
+                              self.current_linecolor, ogstream=cut_ogstream)
         # save PragPicker result for plotting
         self.p_data.refPpick = refPpick
         self.p_results.mpp = refPpick.getpick()
@@ -915,7 +952,7 @@ class AutopickStation(object):
             "minFMSNR"]:
             # if SNR is high enough, try to determine first motion of onset
             self.set_current_figure('fm_picker')
-            self.p_results.fm = fmpicker(self.zstream, z_copy, self.pickparams["fmpickwin"], self.p_results.mpp,
+            self.p_results.fm = fmpicker(self.zstream.copy(), z_copy, self.pickparams["fmpickwin"], self.p_results.mpp,
                                          self.iplot, self.current_figure, self.current_linecolor)
         msg = "autopickstation: P-weight: {}, SNR: {}, SNR[dB]: {}, Polarity: {}"
         msg = msg.format(self.p_results.weight, self.p_results.snr, self.p_results.snrdb, self.p_results.fm)
@@ -1127,9 +1164,11 @@ class AutopickStation(object):
                   ''.format(self.s_results.weight, self.s_results.snr, self.s_results.snrdb))
 
     def pick_s_phase(self):
-
+        if get_bool(self.pickparams.get('use_taup')) is True:
-        # determine time window for calculating CF after P onset
+            cuttimesh = (self.pickparams.get('sstart'), self.pickparams.get('sstop'))
-        cuttimesh = self._calculate_cuttimes(type='S', iteration=1)
+        else:
+            # determine time window for calculating CF after P onset
+            cuttimesh = self._calculate_cuttimes(type='S', iteration=1)
 
         # calculate autoregressive CF
         self.arhcf1 = self._calculate_autoregressive_cf_s_pick(cuttimesh)
@@ -1137,11 +1176,14 @@ class AutopickStation(object):
         # calculate AIC cf
         haiccf = self._calculate_aic_cf_s_pick(cuttimesh)
 
+        # get the original waveform stream cut to identical length as CF for plotting
+        ogstream = haiccf.getDataArray(haiccf.getCut())
+
         # get preliminary onset time from AIC cf
         self.set_current_figure('aicARHfig')
         aicarhpick = AICPicker(haiccf, self.pickparams["tsnrh"], self.pickparams["pickwinS"], self.iplot,
                                Tsmooth=self.pickparams["aictsmoothS"], fig=self.current_figure,
-                               linecolor=self.current_linecolor)
+                               linecolor=self.current_linecolor, ogstream=ogstream)
         # save pick for later plotting
         self.aicarhpick = aicarhpick
 

pylot/core/pick/charfuns.py

@@ -16,10 +16,15 @@ autoregressive prediction: application ot local and regional distances, Geophys.
 
 :author: MAGS2 EP3 working group
 """
-
 import numpy as np
+try:
+    from scipy.signal import tukey
+except ImportError:
+    from scipy.signal.windows import tukey
+
 from obspy.core import Stream
-from scipy import signal
+
+from pylot.core.pick.utils import PickingFailedException
 
 
 class CharacteristicFunction(object):
@@ -55,7 +60,7 @@ class CharacteristicFunction(object):
         self.setOrder(order)
         self.setFnoise(fnoise)
         self.setARdetStep(t2)
-        self.calcCF(self.getDataArray())
+        self.calcCF()
         self.arpara = np.array([])
         self.xpred = np.array([])
 
@@ -207,17 +212,15 @@ class CharacteristicFunction(object):
             data = self.orig_data.copy()
             return data
 
-    def calcCF(self, data=None):
+    def calcCF(self):
-        self.cf = data
+        pass
 
 
 class AICcf(CharacteristicFunction):
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate the Akaike Information Criterion (AIC) after Maeda (1985).
-        :param data: data, time series (whether seismogram or CF)
-        :type data: tuple
         :return: AIC function
         :rtype:
         """
@@ -226,7 +229,7 @@ class AICcf(CharacteristicFunction):
         ind = np.where(~np.isnan(xnp))[0]
         if ind.size:
             xnp[:ind[0]] = xnp[ind[0]]
-        xnp = signal.tukey(len(xnp), alpha=0.05) * xnp
+        xnp = tukey(len(xnp), alpha=0.05) * xnp
         xnp = xnp - np.mean(xnp)
         datlen = len(xnp)
         k = np.arange(1, datlen)
@@ -255,13 +258,11 @@ class HOScf(CharacteristicFunction):
         """
         super(HOScf, self).__init__(data, cut, pickparams["tlta"], pickparams["hosorder"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate skewness (statistics of order 3) or kurtosis
         (statistics of order 4), using one long moving window, as published
-        in Kueperkoch et al. (2010).
+        in Kueperkoch et al. (2010), or order 2, i.e. STA/LTA.
-        :param data: data, time series (whether seismogram or CF)
-        :type data: tuple
         :return: HOS cf
         :rtype:
         """
@@ -316,12 +317,10 @@ class ARZcf(CharacteristicFunction):
         super(ARZcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Parorder"],
                                     fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         function used to calculate the AR prediction error from a single vertical trace. Can be used to pick
         P onsets.
-        :param data:
-        :type data: ~obspy.core.stream.Stream
         :return: ARZ cf
         :rtype:
         """
@@ -452,14 +451,12 @@ class ARHcf(CharacteristicFunction):
         super(ARHcf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
                                     fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate a characteristic function using autoregressive modelling of the waveform of
         both horizontal traces.
         The waveform is predicted in a moving time window using the calculated AR parameters. The difference
         between the predicted and the actual waveform servers as a characteristic function.
-        :param data: wavefor stream
-        :type data: ~obspy.core.stream.Stream
         :return: ARH cf
         :rtype:
         """
@@ -467,6 +464,9 @@ class ARHcf(CharacteristicFunction):
         print('Calculating AR-prediction error from both horizontal traces ...')
 
         xnp = self.getDataArray(self.getCut())
+        if len(xnp[0]) == 0:
+            raise PickingFailedException('calcCF: Found empty data trace for cut times. Return')
+
         n0 = np.isnan(xnp[0].data)
         if len(n0) > 1:
             xnp[0].data[n0] = 0
@@ -605,14 +605,12 @@ class AR3Ccf(CharacteristicFunction):
         super(AR3Ccf, self).__init__(data, cut, t1=t1, t2=t2, order=pickparams["Sarorder"],
                                      fnoise=pickparams["addnoise"])
 
-    def calcCF(self, data):
+    def calcCF(self):
         """
         Function to calculate a characteristic function using autoregressive modelling of the waveform of
         all three traces.
         The waveform is predicted in a moving time window using the calculated AR parameters. The difference
         between the predicted and the actual waveform servers as a characteristic function
-        :param data: stream holding all three traces
-        :type data: ~obspy.core.stream.Stream
         :return: AR3C cf
         :rtype:
         """

pylot/core/pick/picker.py

@@ -37,7 +37,8 @@ class AutoPicker(object):
 
     warnings.simplefilter('ignore')
 
-    def __init__(self, cf, TSNR, PickWindow, iplot=0, aus=None, Tsmooth=None, Pick1=None, fig=None, linecolor='k'):
+    def __init__(self, cf, TSNR, PickWindow, iplot=0, aus=None, Tsmooth=None, Pick1=None,
+                 fig=None, linecolor='k', ogstream=None):
         """
         Create AutoPicker object
         :param cf: characteristic function, on which the picking algorithm is applied
@@ -59,12 +60,15 @@ class AutoPicker(object):
         :type fig: `~matplotlib.figure.Figure`
         :param linecolor: matplotlib line color string
         :type linecolor: str
+        :param ogstream: original stream (waveform), e.g. for plotting purposes
+        :type ogstream: `~obspy.core.stream.Stream`
         """
 
         assert isinstance(cf, CharacteristicFunction), "%s is not a CharacteristicFunction object" % str(cf)
         self._linecolor = linecolor
         self._pickcolor_p = 'b'
         self.cf = cf.getCF()
+        self.ogstream = ogstream
         self.Tcf = cf.getTimeArray()
         self.Data = cf.getXCF()
         self.dt = cf.getIncrement()
@@ -173,12 +177,14 @@ class AICPicker(AutoPicker):
         nn = np.isnan(self.cf)
         if len(nn) > 1:
             self.cf[nn] = 0
-        # taper AIC-CF to get rid off side maxima
+        # taper AIC-CF to get rid of side maxima
         tap = np.hanning(len(self.cf))
         aic = tap * self.cf + max(abs(self.cf))
         # smooth AIC-CF
         ismooth = int(round(self.Tsmooth / self.dt))
-        aicsmooth = np.zeros(len(aic))
+        # MP MP better start with original data than zeros if array shall be smoothed, created artificial value before
+        # when starting with i in range(1...) loop below and subtracting offset afterwards
+        aicsmooth = np.copy(aic)
         if len(aic) < ismooth:
             print('AICPicker: Tsmooth larger than CF!')
             return
@@ -188,7 +194,7 @@ class AICPicker(AutoPicker):
                     ii1 = i - ismooth
                     aicsmooth[i] = aicsmooth[i - 1] + (aic[i] - aic[ii1]) / ismooth
                 else:
-                    aicsmooth[i] = np.mean(aic[1: i])
+                    aicsmooth[i] = np.mean(aic[0: i]) # MP MP created np.nan for i=1
         # remove offset in AIC function
         offset = abs(min(aic) - min(aicsmooth))
         aicsmooth = aicsmooth - offset
@@ -197,7 +203,7 @@ class AICPicker(AutoPicker):
         # minimum in AIC function 
         icfmax = np.argmax(cf)
 
-        # MP MP testing threshold
+        # TODO: If this shall be kept, maybe add thresh_factor to pylot parameters
         thresh_hit = False
         thresh_factor = 0.7
         thresh = thresh_factor * cf[icfmax]
@@ -209,7 +215,6 @@ class AICPicker(AutoPicker):
                 if sample <= cf[index - 1]:
                     icfmax = index - 1
                     break
-        # MP MP ---
 
         # find minimum in AIC-CF front of maximum of HOS/AR-CF
         lpickwindow = int(round(self.PickWindow / self.dt))
@@ -315,16 +320,7 @@ class AICPicker(AutoPicker):
                             plt.close(fig)
                     return
                 iislope = islope[0][0:imax + 1]
-            # MP MP change slope calculation
+            dataslope = self.Data[0].data[iislope]
-            # get all maxima of aicsmooth
-            iaicmaxima = argrelmax(aicsmooth)[0]
-            # get first index of maximum after pickindex (indices saved in iaicmaxima)
-            aicmax = iaicmaxima[np.where(iaicmaxima > pickindex)[0]]
-            if len(aicmax) > 0:
-                iaicmax = aicmax[0]
-            else:
-                iaicmax = -1
-            dataslope = aicsmooth[pickindex: iaicmax]
             # calculate slope as polynomal fit of order 1
             xslope = np.arange(0, len(dataslope), 1)
             try:
@@ -335,8 +331,8 @@ class AICPicker(AutoPicker):
                 else:
                     self.slope = 1 / (len(dataslope) * self.Data[0].stats.delta) * (datafit[-1] - datafit[0])
                     # normalize slope to maximum of cf to make it unit independent
-                    self.slope /= aicsmooth[iaicmax]
+                    self.slope /= self.Data[0].data[icfmax]
-            except ValueError as e:
+            except Exception as e:
                 print("AICPicker: Problems with data fitting! {}".format(e))
 
         else:
@@ -355,6 +351,12 @@ class AICPicker(AutoPicker):
                 self.Tcf = self.Tcf[0:len(self.Tcf) - 1]
             ax1.plot(self.Tcf, cf / max(cf), color=self._linecolor, linewidth=0.7, label='(HOS-/AR-) Data')
             ax1.plot(self.Tcf, aicsmooth / max(aicsmooth), 'r', label='Smoothed AIC-CF')
+            # plot the original waveform also for evaluation of the CF and pick
+            if self.ogstream:
+                data = self.ogstream[0].data
+                if len(data) == len(self.Tcf):
+                    ax1.plot(self.Tcf, 0.5 * data / max(data), 'k', label='Seismogram', alpha=0.3, zorder=0,
+                             lw=0.5)
             if self.Pick is not None:
                 ax1.plot([self.Pick, self.Pick], [-0.1, 0.5], 'b', linewidth=2, label='AIC-Pick')
             ax1.set_xlabel('Time [s] since %s' % self.Data[0].stats.starttime)
@@ -375,7 +377,7 @@ class AICPicker(AutoPicker):
                             label='Signal Window')
                 ax2.axvspan(self.Tcf[iislope[0]], self.Tcf[iislope[-1]], color='g', alpha=0.2, lw=0,
                             label='Slope Window')
-                ax2.plot(self.Tcf[pickindex: iaicmax], datafit, 'g', linewidth=2,
+                ax2.plot(self.Tcf[iislope], datafit, 'g', linewidth=2,
                          label='Slope')  # MP MP changed temporarily!
 
                 if self.slope is not None:

pylot/core/pick/utils.py

@@ -15,7 +15,7 @@ import numpy as np
 from obspy.core import Stream, UTCDateTime
 from scipy.signal import argrelmax
 
-from pylot.core.util.utils import get_Bool, get_None, SetChannelComponents
+from pylot.core.util.utils import get_bool, get_none, SetChannelComponents, common_range
 
 
 def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecolor='k'):
@@ -62,8 +62,8 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
     plt_flag = 0
     try:
         iplot = int(iplot)
-    except:
+    except ValueError:
-        if get_Bool(iplot):
+        if get_bool(iplot):
             iplot = 2
         else:
             iplot = 0
@@ -136,7 +136,7 @@ def earllatepicker(X, nfac, TSNR, Pick1, iplot=0, verbosity=1, fig=None, linecol
     PickError = symmetrize_error(diffti_te, diffti_tl)
 
     if iplot > 1:
-        if get_None(fig) is None:
+        if get_none(fig) is None:
             fig = plt.figure()  # iplot)
             plt_flag = 1
         fig._tight = True
@@ -344,7 +344,7 @@ def fmpicker(Xraw, Xfilt, pickwin, Pick, iplot=0, fig=None, linecolor='k'):
         print("fmpicker: Found polarity %s" % FM)
 
     if iplot > 1:
-        if get_None(fig) is None:
+        if get_none(fig) is None:
             fig = plt.figure()  # iplot)
             plt_flag = 1
         fig._tight = True
@@ -816,7 +816,7 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
     try:
         iplot = int(iplot)
     except:
-        if get_Bool(iplot):
+        if get_bool(iplot):
             iplot = 2
         else:
             iplot = 0
@@ -828,14 +828,22 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
     if len(X) > 1:
         # all three components available
         # make sure, all components have equal lengths
-        ilen = min([len(X[0].data), len(X[1].data), len(X[2].data)])
+        earliest_starttime = min(tr.stats.starttime for tr in X)
-        x1 = X[0][0:ilen]
+        cuttimes = common_range(X)
-        x2 = X[1][0:ilen]
+        X = X.slice(cuttimes[0], cuttimes[1])
-        x3 = X[2][0:ilen]
+        x1, x2, x3 = X[:3]
+
+        if not (len(x1) == len(x2) == len(x3)):
+            raise PickingFailedException('checksignallength: unequal lengths of components!')
+
         # get RMS trace
         rms = np.sqrt((np.power(x1, 2) + np.power(x2, 2) + np.power(x3, 2)) / 3)
+        ilen = len(rms)
+        dt = earliest_starttime - X[0].stats.starttime
+        pick -= dt
     else:
         x1 = X[0].data
+        x2 = x3 = None
         ilen = len(x1)
         rms = abs(x1)
 
@@ -868,12 +876,16 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
             returnflag = 0
 
     if iplot > 1:
-        if get_None(fig) is None:
+        if get_none(fig) is None:
             fig = plt.figure()  # iplot)
             plt_flag = 1
         fig._tight = True
         ax = fig.add_subplot(111)
         ax.plot(t, rms, color=linecolor, linewidth=0.7, label='RMS Data')
+        ax.plot(t, x1, 'k', alpha=0.3, lw=0.3, zorder=0)
+        if x2 is not None and x3 is not None:
+            ax.plot(t, x2, 'r', alpha=0.3, lw=0.3, zorder=0)
+            ax.plot(t, x3, 'g', alpha=0.3, lw=0.3, zorder=0)
         ax.axvspan(t[inoise[0]], t[inoise[-1]], color='y', alpha=0.2, lw=0, label='Noise Window')
         ax.axvspan(t[isignal[0]], t[isignal[-1]], color='b', alpha=0.2, lw=0, label='Signal Window')
         ax.plot([t[isignal[0]], t[isignal[len(isignal) - 1]]],
@@ -883,6 +895,7 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
         ax.set_xlabel('Time [s] since %s' % X[0].stats.starttime)
         ax.set_ylabel('Counts')
         ax.set_title('Check for Signal Length, Station %s' % X[0].stats.station)
+        ax.set_xlim(pickparams["pstart"], pickparams["pstop"])
         ax.set_yticks([])
         if plt_flag == 1:
             fig.show()
@@ -890,6 +903,8 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
                 input()
             except SyntaxError:
                 pass
+            except EOFError:
+                pass
             plt.close(fig)
 
     return returnflag
@@ -1130,7 +1145,7 @@ def checkZ4S(X, pick, pickparams, iplot, fig=None, linecolor='k'):
     try:
         iplot = int(iplot)
     except:
-        if get_Bool(iplot):
+        if get_bool(iplot):
             iplot = 2
         else:
             iplot = 0
@@ -1211,14 +1226,14 @@ def checkZ4S(X, pick, pickparams, iplot, fig=None, linecolor='k'):
             t = np.linspace(diff_dict[key], trace.stats.endtime - trace.stats.starttime + diff_dict[key],
                             trace.stats.npts)
             if i == 0:
-                if get_None(fig) is None:
+                if get_none(fig) is None:
                     fig = plt.figure()  # self.iplot) ### WHY? MP MP
                     plt_flag = 1
                 ax1 = fig.add_subplot(3, 1, i + 1)
                 ax = ax1
                 ax.set_title('CheckZ4S, Station %s' % zdat[0].stats.station)
             else:
-                if get_None(fig) is None:
+                if get_none(fig) is None:
                     fig = plt.figure()  # self.iplot) ### WHY? MP MP
                     plt_flag = 1
                 ax = fig.add_subplot(3, 1, i + 1, sharex=ax1)
@@ -1320,7 +1335,7 @@ def get_quality_class(uncertainty, weight_classes):
     :return: quality of pick (0-4)
     :rtype: int
     """
-    if not uncertainty: return max(weight_classes)
+    if not uncertainty: return len(weight_classes)
     try:
         # create generator expression containing all indices of values in weight classes that are >= than uncertainty.
         # call next on it once to receive first value
@@ -1332,22 +1347,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
@@ -1499,7 +1498,7 @@ def get_pickparams(pickparam):
     first_motion_params = dict(zip(first_motion_names, fm_parameter_values))
     signal_length_params = dict(zip(signal_length_names, sl_parameter_values))
 
-    p_params['use_taup'] = get_Bool(p_params['use_taup'])
+    p_params['use_taup'] = get_bool(p_params['use_taup'])
 
     return p_params, s_params, first_motion_params, signal_length_params
 
@@ -1508,7 +1507,7 @@ def getQualityFromUncertainty(uncertainty, Errors):
     # set initial quality to 4 (worst) and change only if one condition is hit
     quality = 4
 
-    if get_None(uncertainty) is None:
+    if get_none(uncertainty) is None:
         return quality
 
     if uncertainty <= Errors[0]:
@@ -1532,3 +1531,17 @@ if __name__ == '__main__':
     import doctest
 
     doctest.testmod()
+
+
+class PickingFailedException(Exception):
+    """
+    Raised when picking fails due to missing values etc.
+    """
+    pass
+
+
+class MissingTraceException(ValueError):
+    """
+    Used to indicate missing traces in a obspy.core.stream.Stream object
+    """
+    pass

pylot/core/util/array_map.py

@@ -1,18 +1,24 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
+import copy
 import traceback
 
 import cartopy.crs as ccrs
 import cartopy.feature as cf
+from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter
 import matplotlib
 import matplotlib.patheffects as PathEffects
 import matplotlib.pyplot as plt
 import numpy as np
 import obspy
-from PySide2 import QtWidgets
+from PySide2 import QtWidgets, QtGui
 from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
 from mpl_toolkits.axes_grid1.inset_locator import inset_axes
+from obspy import UTCDateTime
+
+from pylot.core.io.data import Data
+from pylot.core.io.phases import picksdict_from_picks
+from pylot.core.util.utils import identifyPhaseID
 from scipy.interpolate import griddata
 
 from pylot.core.pick.utils import get_quality_class
@@ -23,10 +29,10 @@ matplotlib.use('Qt5Agg')
 
 class MplCanvas(FigureCanvas):
 
-    def __init__(self, parent=None, extern_axes=None, width=5, height=4, dpi=100):
+    def __init__(self, extern_axes=None, projection=None, width=15, height=5, dpi=100):
         if extern_axes is None:
             self.fig = plt.figure(figsize=(width, height), dpi=dpi)
-            self.axes = self.fig.add_subplot(111)
+            self.axes = self.fig.add_subplot(111, projection=projection)
         else:
             self.fig = extern_axes.figure
             self.axes = extern_axes
@@ -37,6 +43,7 @@ class MplCanvas(FigureCanvas):
 class Array_map(QtWidgets.QWidget):
     def __init__(self, parent, metadata, parameter=None, axes=None, annotate=True, pointsize=25.,
                  linewidth=1.5, width=5e6, height=2e6):
+
         QtWidgets.QWidget.__init__(self, parent=parent)
 
         assert (parameter is not None or parent is not None), 'either parent or parameter has to be set'
@@ -54,29 +61,36 @@ class Array_map(QtWidgets.QWidget):
         self.uncertainties = None
         self.autopicks_dict = None
         self.hybrids_dict = None
+        self.picks_dict_ref = None
+        self.autopicks_dict_ref = None
         self.eventLoc = None
         self.parameter = parameter if parameter else parent._inputs
 
         self.picks_rel = {}
+        self.picks_rel_mean_corrected = {}
         self.marked_stations = []
         self.highlighted_stations = []
 
         # call functions to draw everything
+        self.projection = ccrs.PlateCarree()
         self.init_graphics()
+        self.ax = self.canvas.axes
+        self.ax.set_adjustable('datalim')
+
         self.init_stations()
         self.init_crtpyMap()
         self.init_map()
+
         # set original map limits to fall back on when home button is pressed
-        self.org_xlim = self.canvas.axes.get_xlim()
+        self.org_xlim = self.ax.get_xlim()
-        self.org_ylim = self.canvas.axes.get_ylim()
+        self.org_ylim = self.ax.get_ylim()
+
         # initial map without event
-        self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
+        self.ax.set_xlim(self.org_xlim[0], self.org_xlim[1])
-        self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
+        self.ax.set_ylim(self.org_ylim[0], self.org_ylim[1])
 
         self._style = None if not hasattr(parent, '_style') else parent._style
 
-        self.show()
-
     def init_map(self):
         self.init_colormap()
         self.connectSignals()
@@ -89,23 +103,26 @@ class Array_map(QtWidgets.QWidget):
         # initialize figure elements
 
         if self.extern_plot_axes is None:
-            self.canvas = MplCanvas(self)
+            self.canvas = MplCanvas(projection=self.projection)
-            self.plotWidget = FigureCanvas(self.canvas.fig)
         else:
-            self.canvas = MplCanvas(self, extern_axes=self.extern_plot_axes)
+            self.canvas = MplCanvas(extern_axes=self.extern_plot_axes)
-            self.plotWidget = FigureCanvas(self.canvas.fig)
+
+        self.plotWidget = self.canvas
 
         # initialize GUI elements
         self.status_label = QtWidgets.QLabel()
         self.map_reset_button = QtWidgets.QPushButton('Reset Map View')
         self.save_map_button = QtWidgets.QPushButton('Save Map')
+        self.load_reference_picks_button = QtWidgets.QPushButton('Load reference picks.')
         self.go2eq_button = QtWidgets.QPushButton('Go to Event Location')
+        self.subtract_mean_cb = QtWidgets.QCheckBox('Subtract mean')
+        self.subtract_ref_cb = QtWidgets.QCheckBox('Subtract reference onsets')
 
         self.main_box = QtWidgets.QVBoxLayout()
         self.setLayout(self.main_box)
 
         self.top_row = QtWidgets.QHBoxLayout()
-        self.main_box.addLayout(self.top_row, 1)
+        self.main_box.addLayout(self.top_row, 0)
 
         self.comboBox_phase = QtWidgets.QComboBox()
         self.comboBox_phase.insertItem(0, 'P')
@@ -124,8 +141,8 @@ class Array_map(QtWidgets.QWidget):
         self.cmaps_box = QtWidgets.QComboBox()
         self.cmaps_box.setMaxVisibleItems(20)
         [self.cmaps_box.addItem(map_name) for map_name in sorted(plt.colormaps())]
-        # try to set to hsv as default
+        # try to set to viridis as default
-        self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('hsv'))
+        self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('viridis'))
 
         self.top_row.addWidget(QtWidgets.QLabel('Select a phase: '))
         self.top_row.addWidget(self.comboBox_phase)
@@ -138,14 +155,17 @@ class Array_map(QtWidgets.QWidget):
         self.top_row.addWidget(self.auto_refresh_box)
         self.top_row.addWidget(self.refresh_button)
 
-        self.main_box.addWidget(self.plotWidget, 1)
+        self.main_box.addWidget(self.plotWidget, 10)
 
         self.bot_row = QtWidgets.QHBoxLayout()
-        self.main_box.addLayout(self.bot_row, 0.3)
+        self.main_box.addLayout(self.bot_row, 0)
         self.bot_row.addWidget(QtWidgets.QLabel(''), 5)
         self.bot_row.addWidget(self.map_reset_button, 2)
         self.bot_row.addWidget(self.go2eq_button, 2)
         self.bot_row.addWidget(self.save_map_button, 2)
+        self.bot_row.addWidget(self.load_reference_picks_button, 2)
+        self.bot_row.addWidget(self.subtract_mean_cb, 0)
+        self.bot_row.addWidget(self.subtract_ref_cb, 0)
         self.bot_row.addWidget(self.status_label, 5)
 
     def init_colormap(self):
@@ -153,14 +173,12 @@ class Array_map(QtWidgets.QWidget):
         self.init_lat_lon_grid()
 
     def init_crtpyMap(self):
-        self.canvas.axes.cla()
+        self.ax.add_feature(cf.LAND)
-        self.canvas.axes = plt.axes(projection=ccrs.PlateCarree())
+        self.ax.add_feature(cf.OCEAN)
-        self.canvas.axes.add_feature(cf.LAND)
+        self.ax.add_feature(cf.COASTLINE, linewidth=1, edgecolor='gray')
-        self.canvas.axes.add_feature(cf.OCEAN)
+        self.ax.add_feature(cf.BORDERS, alpha=0.7)
-        self.canvas.axes.add_feature(cf.COASTLINE, linewidth=1, edgecolor='gray')
+        self.ax.add_feature(cf.LAKES, alpha=0.7)
-        self.canvas.axes.add_feature(cf.BORDERS, alpha=0.7)
+        self.ax.add_feature(cf.RIVERS, linewidth=1)
-        self.canvas.axes.add_feature(cf.LAKES, alpha=0.7)
-        self.canvas.axes.add_feature(cf.RIVERS, linewidth=1)
 
         # parallels and meridians
         self.add_merid_paral()
@@ -168,12 +186,8 @@ class Array_map(QtWidgets.QWidget):
         self.canvas.fig.tight_layout()
 
     def add_merid_paral(self):
-        self.gridlines = self.canvas.axes.gridlines(draw_labels=False, alpha=0.6, color='gray',
+        self.gridlines = self.ax.gridlines(draw_labels=False, alpha=0.6, color='gray',
-                                                    linewidth=self.linewidth / 2, zorder=7)
+                                           linewidth=self.linewidth / 2, zorder=7, crs=ccrs.PlateCarree())
-        # TODO: current cartopy version does not support label removal. Devs are working on it.
-        #  Should be fixed in coming cartopy versions
-        # self.gridlines.xformatter = LONGITUDE_FORMATTER
-        # self.gridlines.yformatter = LATITUDE_FORMATTER
 
     def remove_merid_paral(self):
         if len(self.gridlines.xline_artists):
@@ -181,24 +195,24 @@ class Array_map(QtWidgets.QWidget):
             self.gridlines.yline_artists[0].remove()
 
     def org_map_view(self):
-        self.canvas.axes.set_xlim(self.org_xlim[0], self.org_xlim[1])
+        self.ax.set_xlim(self.org_xlim[0], self.org_xlim[1])
-        self.canvas.axes.set_ylim(self.org_ylim[0], self.org_ylim[1])
+        self.ax.set_ylim(self.org_ylim[0], self.org_ylim[1])
         # parallels and meridians
-        self.remove_merid_paral()
+        #self.remove_merid_paral()
-        self.add_merid_paral()
+        #self.add_merid_paral()
 
-        self.canvas.axes.figure.canvas.draw_idle()
+        self.canvas.draw_idle()
 
     def go2eq(self):
         if self.eventLoc:
             lats, lons = self.eventLoc
-            self.canvas.axes.set_xlim(lons - 10, lons + 10)
+            self.ax.set_xlim(lons - 10, lons + 10)
-            self.canvas.axes.set_ylim(lats - 5, lats + 5)
+            self.ax.set_ylim(lats - 5, lats + 5)
             # parallels and meridians
-            self.remove_merid_paral()
+            #self.remove_merid_paral()
-            self.add_merid_paral()
+            #self.add_merid_paral()
 
-            self.canvas.axes.figure.canvas.draw_idle()
+            self.canvas.draw_idle()
 
         else:
             self.status_label.setText('No event information available')
@@ -212,6 +226,9 @@ class Array_map(QtWidgets.QWidget):
         self.map_reset_button.clicked.connect(self.org_map_view)
         self.go2eq_button.clicked.connect(self.go2eq)
         self.save_map_button.clicked.connect(self.saveFigure)
+        self.load_reference_picks_button.clicked.connect(self.load_reference_picks)
+        self.subtract_mean_cb.stateChanged.connect(self.toggle_subtract_mean)
+        self.subtract_ref_cb.stateChanged.connect(self.toggle_subtract_ref)
 
         self.plotWidget.mpl_connect('motion_notify_event', self.mouse_moved)
         self.plotWidget.mpl_connect('scroll_event', self.mouse_scroll)
@@ -220,21 +237,32 @@ class Array_map(QtWidgets.QWidget):
 
     # set mouse events -----------------------------------------------------
     def mouse_moved(self, event):
-        if not event.inaxes == self.canvas.axes:
+        if not event.inaxes == self.ax:
             return
+        else:
+            cont, inds = self.sc.contains(event)
         lat = event.ydata
         lon = event.xdata
-        self.status_label.setText('Latitude: {:3.5f}, Longitude: {:3.5f}'.format(lat, lon))
+        text = f'Longitude: {lon:3.3f}, Latitude: {lat:3.3f}'
+
+        if cont:
+            indices = inds['ind']
+            text += ' | Station: ' if len(indices) == 1 else ' | Stations: '
+            text += ' - '.join([self._station_onpick_ids[index] for index in indices[:5]])
+            if len(indices) > 5:
+                text += '...'
+
+        self.status_label.setText(text)
 
     def mouse_scroll(self, event):
-        if not event.inaxes == self.canvas.axes:
+        if not event.inaxes == self.ax:
             return
 
         zoom = {'up': 1. / 2., 'down': 2.}
 
         if event.button in zoom:
-            xlim = self.canvas.axes.get_xlim()
+            xlim = self.ax.get_xlim()
-            ylim = self.canvas.axes.get_ylim()
+            ylim = self.ax.get_ylim()
 
             x, y = event.xdata, event.ydata
 
@@ -246,24 +274,24 @@ class Array_map(QtWidgets.QWidget):
             yb = y - 0.5 * ydiff
             yt = y + 0.5 * ydiff
 
-            self.canvas.axes.set_xlim(xl, xr)
+            self.ax.set_xlim(xl, xr)
-            self.canvas.axes.set_ylim(yb, yt)
+            self.ax.set_ylim(yb, yt)
             # parallels and meridians
-            self.remove_merid_paral()
+            #self.remove_merid_paral()
-            self.add_merid_paral()
+            #self.add_merid_paral()
 
-            self.canvas.axes.figure.canvas.draw_idle()
+            self.ax.figure.canvas.draw_idle()
 
     def mouseLeftPress(self, event):
-        if not event.inaxes == self.canvas.axes:
+        if not event.inaxes == self.ax:
             return
         self.map_x = event.xdata
         self.map_y = event.ydata
-        self.map_xlim = self.canvas.axes.get_xlim()
+        self.map_xlim = self.ax.get_xlim()
-        self.map_ylim = self.canvas.axes.get_ylim()
+        self.map_ylim = self.ax.get_ylim()
 
     def mouseLeftRelease(self, event):
-        if not event.inaxes == self.canvas.axes:
+        if not event.inaxes == self.ax:
             return
         new_x = event.xdata
         new_y = event.ydata
@@ -271,18 +299,21 @@ class Array_map(QtWidgets.QWidget):
         dx = new_x - self.map_x
         dy = new_y - self.map_y
 
-        self.canvas.axes.set_xlim((self.map_xlim[0] - dx, self.map_xlim[1] - dx))
+        self.ax.set_xlim((self.map_xlim[0] - dx, self.map_xlim[1] - dx))
-        self.canvas.axes.set_ylim(self.map_ylim[0] - dy, self.map_ylim[1] - dy)
+        self.ax.set_ylim(self.map_ylim[0] - dy, self.map_ylim[1] - dy)
         # parallels and meridians
-        self.remove_merid_paral()
+        #self.remove_merid_paral()
-        self.add_merid_paral()
+        #self.add_merid_paral()
 
-        self.canvas.axes.figure.canvas.draw_idle()
+        self.ax.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)
@@ -354,14 +385,11 @@ class Array_map(QtWidgets.QWidget):
     def get_max_from_stations(self, key):
         return self._from_dict(max, key)
 
-    def get_min_from_picks(self):
+    def get_selected_pick_type(self):
-        return min(self.picks_rel.values())
+        return self.comboBox_am.currentText().split(' ')[0]
-
-    def get_max_from_picks(self):
-        return max(self.picks_rel.values())
 
     def current_picks_dict(self):
-        picktype = self.comboBox_am.currentText().split(' ')[0]
+        picktype = self.get_selected_pick_type()
         auto_manu = {'auto': self.autopicks_dict,
                      'manual': self.picks_dict,
                      'hybrid': self.hybrids_dict}
@@ -385,7 +413,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
@@ -397,22 +432,37 @@ class Array_map(QtWidgets.QWidget):
                     print('Cannot display pick for station {}. Reason: {}'.format(station_name, e))
             return picks, uncertainties
 
-        def get_picks_rel(picks):
+        def get_picks_rel(picks, func=min):
             picks_rel = {}
             picks_utc = []
             for pick in picks.values():
-                if type(pick) is obspy.core.utcdatetime.UTCDateTime:
+                if type(pick) is UTCDateTime:
-                    picks_utc.append(pick)
+                    picks_utc.append(pick.timestamp)
             if picks_utc:
-                self._earliest_picktime = min(picks_utc)
+                self._reference_picktime = UTCDateTime(func(picks_utc))
                 for st_id, pick in picks.items():
-                    if type(pick) is obspy.core.utcdatetime.UTCDateTime:
+                    if type(pick) is UTCDateTime:
-                        pick -= self._earliest_picktime
+                        pick -= self._reference_picktime
                     picks_rel[st_id] = pick
             return picks_rel
 
+        def get_picks_rel_mean_corr(picks):
+            return get_picks_rel(picks, func=np.nanmean)
+
         self.picks, self.uncertainties = get_picks(self.stations_dict)
         self.picks_rel = get_picks_rel(self.picks)
+        self.picks_rel_mean_corrected = get_picks_rel_mean_corr(self.picks)
+
+    def toggle_subtract_mean(self):
+        if self.subtract_mean_cb.isChecked():
+            cmap = 'seismic'
+        else:
+            cmap = 'viridis'
+        self.cmaps_box.setCurrentIndex(self.cmaps_box.findText(cmap))
+        self._refresh_drawings()
+
+    def toggle_subtract_ref(self):
+        self._refresh_drawings()
 
     def init_lat_lon_dimensions(self):
         # init minimum and maximum lon and lat dimensions
@@ -426,7 +476,10 @@ class Array_map(QtWidgets.QWidget):
         self.longrid, self.latgrid = np.meshgrid(lonaxis, lataxis)
 
     def init_picksgrid(self):
-        picks, uncertainties, lats, lons = self.get_picks_lat_lon()
+        rval = self.get_picks_lat_lon()
+        if not rval:
+            return
+        picks, uncertainties, lats, lons = rval
         try:
             self.picksgrid_active = griddata((lats, lons), picks, (self.latgrid, self.longrid), method='linear')
         except Exception as e:
@@ -443,15 +496,21 @@ class Array_map(QtWidgets.QWidget):
         return stations, latitudes, longitudes
 
     def get_picks_lat_lon(self):
+        picks_rel = self.picks_rel_mean_corrected if self.subtract_mean_cb.isChecked() else self.picks_rel
+        picks_rel = self.picks_rel_ref_corrected if self.subtract_ref_cb.isChecked() else picks_rel
+        if not picks_rel:
+            return
         picks = []
         uncertainties = []
         latitudes = []
         longitudes = []
-        for st_id, pick in self.picks_rel.items():
+        for nwst_id, pick in picks_rel.items():
+            if nwst_id not in self.uncertainties or nwst_id not in self.stations_dict:
+                continue
             picks.append(pick)
-            uncertainties.append(self.uncertainties.get(st_id))
+            uncertainties.append(self.uncertainties.get(nwst_id))
-            latitudes.append(self.stations_dict[st_id]['latitude'])
+            latitudes.append(self.stations_dict[nwst_id]['latitude'])
-            longitudes.append(self.stations_dict[st_id]['longitude'])
+            longitudes.append(self.stations_dict[nwst_id]['longitude'])
         return picks, uncertainties, latitudes, longitudes
 
     # plotting -----------------------------------------------------
@@ -459,22 +518,31 @@ class Array_map(QtWidgets.QWidget):
         stat_dict = self.stations_dict['{}.{}'.format(network, station)]
         lat = stat_dict['latitude']
         lon = stat_dict['longitude']
-        self.highlighted_stations.append(self.canvas.axes.scatter(lon, lat, s=self.pointsize, edgecolors=color,
+        self.highlighted_stations.append(self.ax.scatter(lon, lat, s=self.pointsize, edgecolors=color,
-                                                                  facecolors='none', zorder=12,
+                                                         facecolors='none', zorder=12,
-                                                                  transform=ccrs.PlateCarree(), label='deleted'))
+                                                         transform=ccrs.PlateCarree(), label='deleted'))
 
     def openPickDlg(self, ind):
-        data = self._parent.get_data().getWFData()
+        try:
+            wfdata = self._parent.get_data().getWFData()
+        except AttributeError:
+            QtWidgets.QMessageBox.warning(
+                self, "PyLoT Warning",
+                "No waveform data found. Check if they were already loaded in Waveform plot tab."
+            )
+            return
+        wfdata_comp = self._parent.get_data().getAltWFdata()
         for index in ind:
             network, station = self._station_onpick_ids[index].split('.')[:2]
             pyl_mw = self._parent
             try:
-                data = data.select(station=station)
+                wfdata = wfdata.select(station=station)
-                if not data:
+                wfdata_comp = wfdata_comp.select(station=station)
+                if not wfdata:
                     self._warn('No data for station {}'.format(station))
                     return
                 pickDlg = PickDlg(self._parent, parameter=self.parameter,
-                                  data=data, network=network, station=station,
+                                  data=wfdata.copy(), data_compare=wfdata_comp.copy(), network=network, station=station,
                                   picks=self._parent.get_current_event().getPick(station),
                                   autopicks=self._parent.get_current_event().getAutopick(station),
                                   filteroptions=self._parent.filteroptions, metadata=self.metadata,
@@ -505,12 +573,19 @@ class Array_map(QtWidgets.QWidget):
                 print(message, e)
                 print(traceback.format_exc())
 
-    def draw_contour_filled(self, nlevel=50):
+    def draw_contour_filled(self, nlevel=51):
-        levels = np.linspace(self.get_min_from_picks(), self.get_max_from_picks(), nlevel)
+        if self.subtract_mean_cb.isChecked():
+            abs_max = self.get_residuals_absmax()
+            levels = np.linspace(-abs_max, abs_max, nlevel)
+        else:
+            levels = np.linspace(min(self.picks_rel.values()), max(self.picks_rel.values()), nlevel)
 
-        self.contourf = self.canvas.axes.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels,
+        self.contourf = self.ax.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels,
-                                                  linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(),
+                                         linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(),
-                                                  alpha=0.4, zorder=8, cmap=self.get_colormap())
+                                         alpha=0.4, zorder=8, cmap=self.get_colormap())
+
+    def get_residuals_absmax(self):
+        return np.max(np.absolute(list(self.picks_rel_mean_corrected.values())))
 
     def get_colormap(self):
         return plt.get_cmap(self.cmaps_box.currentText())
@@ -518,21 +593,24 @@ class Array_map(QtWidgets.QWidget):
     def scatter_all_stations(self):
         stations, lats, lons = self.get_st_lat_lon_for_plot()
 
-        self.sc = self.canvas.axes.scatter(lons, lats, s=self.pointsize * 3, facecolor='none', marker='.',
+        self.sc = self.ax.scatter(lons, lats, s=self.pointsize * 3, facecolor='none', marker='.',
-                                           zorder=10, picker=True, edgecolor='0.5', label='Not Picked',
+                                  zorder=10, picker=True, edgecolor='0.5', label='Not Picked',
-                                           transform=ccrs.PlateCarree())
+                                  transform=ccrs.PlateCarree())
 
         self.cid = self.plotWidget.mpl_connect('pick_event', self.onpick)
         self._station_onpick_ids = stations
         if self.eventLoc:
             lats, lons = self.eventLoc
-            self.sc_event = self.canvas.axes.scatter(lons, lats, s=5 * self.pointsize, facecolor='red', zorder=11,
+            self.sc_event = self.ax.scatter(lons, lats, s=5 * self.pointsize, facecolor='red', zorder=11,
-                                                     label='Event (might be outside map region)', marker='*',
+                                            label='Event (might be outside map region)', marker='*',
-                                                     edgecolors='black',
+                                            edgecolors='black',
-                                                     transform=ccrs.PlateCarree())
+                                            transform=ccrs.PlateCarree())
 
     def scatter_picked_stations(self):
-        picks, uncertainties, lats, lons = self.get_picks_lat_lon()
+        rval = self.get_picks_lat_lon()
+        if not rval:
+            return
+        picks, uncertainties, lats, lons = rval
         if len(lons) < 1 and len(lats) < 1:
             return
 
@@ -542,8 +620,13 @@ class Array_map(QtWidgets.QWidget):
                           for uncertainty in uncertainties])
 
         cmap = self.get_colormap()
-        self.sc_picked = self.canvas.axes.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap,
+
-                                                  c=picks, zorder=11, label='Picked', transform=ccrs.PlateCarree())
+        vmin = vmax = None
+        if self.subtract_mean_cb.isChecked():
+            vmin, vmax = -self.get_residuals_absmax(), self.get_residuals_absmax()
+
+        self.sc_picked = self.ax.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap, vmin=vmin, vmax=vmax,
+                                         c=picks, zorder=11, label='Picked', transform=ccrs.PlateCarree())
 
     def annotate_ax(self):
         self.annotations = []
@@ -561,20 +644,20 @@ class Array_map(QtWidgets.QWidget):
             if st in self.marked_stations:
                 color = 'red'
             self.annotations.append(
-                self.canvas.axes.annotate(' %s' % st, xy=(x + 0.003, y + 0.003), fontsize=self.pointsize / 4.,
+                self.ax.annotate(f'{st}', xy=(x + 0.003, y + 0.003), fontsize=self.pointsize / 4.,
-                                          fontweight='semibold', color=color, alpha=0.8,
+                                 fontweight='semibold', color=color, alpha=0.8,
-                                          transform=ccrs.PlateCarree(), zorder=14,
+                                 transform=ccrs.PlateCarree(), zorder=14,
-                                          path_effects=[PathEffects.withStroke(
+                                 path_effects=[PathEffects.withStroke(
-                                              linewidth=self.pointsize / 15., foreground='k')]))
+                                     linewidth=self.pointsize / 15., foreground='k')]))
 
-        self.legend = self.canvas.axes.legend(loc=1, framealpha=1)
+        self.legend = self.ax.legend(loc=1, framealpha=1)
         self.legend.set_zorder(100)
         self.legend.get_frame().set_facecolor((1, 1, 1, 0.95))
 
     def add_cbar(self, label):
-        self.cbax_bg = inset_axes(self.canvas.axes, width="6%", height="75%", loc=5)
+        self.cbax_bg = inset_axes(self.ax, width="6%", height="75%", loc=5)
-        cbax = inset_axes(self.canvas.axes, width='2%', height='70%', loc=5)
+        cbax = inset_axes(self.ax, width='2%', height='70%', loc=5)
-        cbar = self.canvas.axes.figure.colorbar(self.sc_picked, cax=cbax)
+        cbar = self.ax.figure.colorbar(self.sc_picked, cax=cbax)
         cbar.set_label(label)
         cbax.yaxis.tick_left()
         cbax.yaxis.set_label_position('left')
@@ -619,7 +702,9 @@ class Array_map(QtWidgets.QWidget):
         if picks_available:
             self.scatter_picked_stations()
             if hasattr(self, 'sc_picked'):
-                self.cbar = self.add_cbar(label='Time relative to first onset ({}) [s]'.format(self._earliest_picktime))
+                self.cbar = self.add_cbar(
+                    label='Time relative to reference onset ({}) [s]'.format(self._reference_picktime)
+                )
             self.comboBox_phase.setEnabled(True)
         else:
             self.comboBox_phase.setEnabled(False)
@@ -680,6 +765,52 @@ class Array_map(QtWidgets.QWidget):
                 fname += '.png'
             self.canvas.fig.savefig(fname)
 
+    def load_reference_picks(self):
+        fname = self._parent.filename_from_action()
+        data_ref = Data(parent=self._parent, evtdata=str(fname))
+        evt_ref = data_ref.get_evt_data()
+        if not evt_ref:
+            return
+        picks_ref = evt_ref.picks
+        if not picks_ref:
+            return
+        picksdict_ref = picksdict_from_picks(evt_ref, nwst_id=False)
+        self.autopicks_dict = picksdict_ref['manual']
+        self.autopicks_dict_ref = picksdict_ref['auto']
+        return True
+
+    @property
+    def picks_rel_ref_corrected(self):
+        picktype = self.get_selected_pick_type()
+        if picktype in ['auto', 'hybrid']:
+            picks_dict = self.autopicks_dict_ref
+        elif picktype == 'manual':
+            picks_dict = self.autopicks_dict
+        else:
+            return
+        if picks_dict:
+            return self.subtract_picks(picks_dict)
+        else:
+            if self.load_reference_picks():
+                return self.picks_rel_ref_corrected
+
+    def subtract_picks(self, picks_dict):
+        current_picks = self.current_picks_dict()
+        subtracted_picks = {}
+        for station, ps_dict in current_picks.items():
+            for phase, pick in ps_dict.items():
+                pick_ref_ps = picks_dict.get(station)
+                if not pick_ref_ps:
+                    continue
+                pick_ref = pick_ref_ps.get(phase)
+                if not pick_ref:
+                    continue
+                mpp = pick['mpp'] - pick_ref['mpp']
+                nwst_id = f'{pick["network"]}.{station}'
+                subtracted_picks[nwst_id] = mpp
+        return subtracted_picks
+
+
     def _warn(self, message):
         self.qmb = QtWidgets.QMessageBox(QtWidgets.QMessageBox.Icon.Warning, 'Warning', message)
         self.qmb.show()

pylot/core/util/dataprocessing.py

@@ -2,6 +2,7 @@
 # -*- coding: utf-8 -*-
 
 import glob
+import logging
 import os
 import sys
 
@@ -26,6 +27,10 @@ class Metadata(object):
         # saves which metadata files are from obspy dmt
         self.obspy_dmt_invs = []
         if inventory:
+            # make sure that no accidental backslashes mess up the path
+            if isinstance(inventory, str):
+                inventory = inventory.replace('\\', '/')
+                inventory = os.path.abspath(inventory)
             if os.path.isdir(inventory):
                 self.add_inventory(inventory)
             if os.path.isfile(inventory):
@@ -54,6 +59,8 @@ class Metadata(object):
         :type path_to_inventory: str
         :return: None
         """
+        path_to_inventory = path_to_inventory.replace('\\', '/')
+        path_to_inventory = os.path.abspath(path_to_inventory)
         assert (os.path.isdir(path_to_inventory)), '{} is no directory'.format(path_to_inventory)
         if path_to_inventory not in self.inventories:
             self.inventories.append(path_to_inventory)
@@ -189,7 +196,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):
@@ -271,7 +282,7 @@ class Metadata(object):
                     continue
             invtype, robj = self._read_metadata_file(os.path.join(path_to_inventory, fname))
             try:
-                robj.get_coordinates(station_seed_id)
+                # robj.get_coordinates(station_seed_id)  # TODO: Commented out, failed with Parser, is this needed?
                 self.inventory_files[fname] = {'invtype': invtype,
                                                'data': robj}
                 if station_seed_id in self.seed_ids.keys():

pylot/core/util/defaults.py

@@ -26,9 +26,7 @@ elif system_name == "Windows":
 # suffix for phase name if not phase identified by last letter (P, p, etc.)
 ALTSUFFIX = ['diff', 'n', 'g', '1', '2', '3']
 
-FILTERDEFAULTS = readDefaultFilterInformation(os.path.join(os.path.expanduser('~'),
+FILTERDEFAULTS = readDefaultFilterInformation()
-                                                           '.pylot',
-                                                           'pylot.in'))
 
 TIMEERROR_DEFAULTS = os.path.join(os.path.expanduser('~'),
                                   '.pylot',

pylot/core/util/event.py

@@ -22,14 +22,11 @@ class Event(ObsPyEvent):
         :param path: path to event directory
         :type path: str
         """
-        # TODO: remove rootpath and database
         self.pylot_id = path.split('/')[-1]
         # initialize super class
         super(Event, self).__init__(resource_id=ResourceIdentifier('smi:local/' + self.pylot_id))
         self.path = path
-        self.database = path.split('/')[-2]
         self.datapath = os.path.split(path)[0]  # path.split('/')[-3]
-        self.rootpath = '/' + os.path.join(*path.split('/')[:-3])
         self.pylot_autopicks = {}
         self.pylot_picks = {}
         self.notes = ''

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/gui.py

@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 import os
+from functools import lru_cache
 
 try:
     import pyqtgraph as pg
@@ -25,14 +26,14 @@ def pick_linestyle_pg(picktype, key):
     :return: Qt line style parameters
     :rtype:
     """
-    linestyles_manu = {'mpp': (QtCore.Qt.SolidLine, 2.),
+    linestyles_manu = {'mpp': (QtCore.Qt.SolidLine, 2),
-                       'epp': (QtCore.Qt.DashLine, 1.),
+                       'epp': (QtCore.Qt.DashLine, 1),
-                       'lpp': (QtCore.Qt.DashLine, 1.),
+                       'lpp': (QtCore.Qt.DashLine, 1),
-                       'spe': (QtCore.Qt.DashLine, 1.)}
+                       'spe': (QtCore.Qt.DashLine, 1)}
-    linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2.),
+    linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2),
-                       'epp': (QtCore.Qt.DashDotLine, 1.),
+                       'epp': (QtCore.Qt.DashDotLine, 1),
-                       'lpp': (QtCore.Qt.DashDotLine, 1.),
+                       'lpp': (QtCore.Qt.DashDotLine, 1),
-                       'spe': (QtCore.Qt.DashDotLine, 1.)}
+                       'spe': (QtCore.Qt.DashDotLine, 1)}
     linestyles = {'manual': linestyles_manu,
                   'auto': linestyles_auto}
     return linestyles[picktype][key]
@@ -56,7 +57,7 @@ def which(program, parameter):
         nllocpath = ":" + parameter.get('nllocbin')
         os.environ['PATH'] += nllocpath
     except Exception as e:
-        print(e.message)
+        print(e)
 
     def is_exe(fpath):
         return os.path.exists(fpath) and os.access(fpath, os.X_OK)
@@ -80,6 +81,7 @@ def which(program, parameter):
     return None
 
 
+@lru_cache(maxsize=128)
 def make_pen(picktype, phase, key, quality):
     """
     Make  PyQtGraph.QPen

pylot/core/util/thread.py

@@ -22,9 +22,11 @@ class Thread(QThread):
         self.abortButton = abortButton
         self.finished.connect(self.hideProgressbar)
         self.showProgressbar()
+        self.old_stdout = None
 
     def run(self):
         if self.redirect_stdout:
+            self.old_stdout = sys.stdout
             sys.stdout = self
         try:
             if self.arg is not None:
@@ -39,7 +41,8 @@ class Thread(QThread):
             exctype, value = sys.exc_info()[:2]
             self._executedErrorInfo = '{} {} {}'. \
                 format(exctype, value, traceback.format_exc())
-        sys.stdout = sys.__stdout__
+        if self.redirect_stdout:
+            sys.stdout = self.old_stdout
 
     def showProgressbar(self):
         if self.progressText:
@@ -96,10 +99,12 @@ class Worker(QRunnable):
         self.progressText = progressText
         self.pb_widget = pb_widget
         self.redirect_stdout = redirect_stdout
+        self.old_stdout = None
 
     @Slot()
     def run(self):
         if self.redirect_stdout:
+            self.old_stdout = sys.stdout
             sys.stdout = self
 
         try:
@@ -112,7 +117,7 @@ class Worker(QRunnable):
             self.signals.result.emit(result)
         finally:
             self.signals.finished.emit('Done')
-        sys.stdout = sys.__stdout__
+        sys.stdout = self.old_stdout
 
     def write(self, text):
         self.signals.message.emit(text)
@@ -144,16 +149,18 @@ class MultiThread(QThread):
         self.progressText = progressText
         self.pb_widget = pb_widget
         self.redirect_stdout = redirect_stdout
+        self.old_stdout = None
         self.finished.connect(self.hideProgressbar)
         self.showProgressbar()
 
     def run(self):
         if self.redirect_stdout:
+            self.old_stdout = sys.stdout
             sys.stdout = self
         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()
@@ -164,7 +171,7 @@ class MultiThread(QThread):
             exc_type, exc_obj, exc_tb = sys.exc_info()
             fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
             print('Exception: {}, file: {}, line: {}'.format(exc_type, fname, exc_tb.tb_lineno))
-        sys.stdout = sys.__stdout__
+        sys.stdout = self.old_stdout
 
     def showProgressbar(self):
         if self.progressText:

pylot/core/util/utils.py

@@ -1,12 +1,15 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
+import glob
 import hashlib
+import logging
 import os
 import platform
 import re
 import subprocess
 import warnings
+from typing import Literal, Tuple, Type
+from functools import lru_cache
 
 import numpy as np
 from obspy import UTCDateTime, read
@@ -18,6 +21,10 @@ from pylot.core.io.inputs import PylotParameter, FilterOptions
 from pylot.core.util.obspyDMT_interface import check_obspydmt_eventfolder
 from pylot.styles import style_settings
 
+Rgba: Type[tuple] = Tuple[int, int, int, int]
+Mplrgba: Type[tuple] = Tuple[float, float, float, float]
+Mplrgbastr: Type[tuple] = Tuple[str, str, str, str]
+
 
 def _pickle_method(m):
     if m.im_self is None:
@@ -37,15 +44,14 @@ def getAutoFilteroptions(phase, parameter):
     return filteroptions
 
 
-def readDefaultFilterInformation(fname):
+def readDefaultFilterInformation():
     """
     Read default filter information from pylot.in file
-    :param fname: path to pylot.in file
-    :type fname: str
     :return: dictionary containing the defailt filter information
     :rtype: dict
     """
-    pparam = PylotParameter(fname)
+    pparam = PylotParameter()
+    pparam.reset_defaults()
     return readFilterInformation(pparam)
 
 
@@ -82,25 +88,6 @@ def fit_curve(x, y):
     return splev, splrep(x, y)
 
 
-def getindexbounds(f, eta):
-    """
-    Get indices of values closest below and above maximum value in an array
-    :param f: array
-    :type f: `~numpy.ndarray`
-    :param eta: look for value in array that is closes to max_value * eta
-    :type eta: float
-    :return: tuple containing index of max value, index of value closest below max value,
-     index of value closest above max value
-    :rtype: (int, int, int)
-    """
-    mi = f.argmax()  # get indices of max values
-    m = max(f)  # get maximum value
-    b = m * eta  #
-    l = find_nearest(f[:mi], b)  # find closest value below max value
-    u = find_nearest(f[mi:], b) + mi  # find closest value above max value
-    return mi, l, u
-
-
 def gen_Pool(ncores=0):
     """
     Generate mulitprocessing pool object utilizing ncores amount of cores
@@ -120,7 +107,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
 
 
@@ -166,11 +153,11 @@ def clims(lim1, lim2):
     """
     takes two pairs of limits and returns one pair of common limts
     :param lim1: limit 1
-    :type lim1: int
+    :type lim1: List[int]
     :param lim2: limit 2
-    :type lim2: int
+    :type lim2: List[int]
     :return: new upper and lower limit common to both given limits
-    :rtype: [int, int]
+    :rtype: List[int]
 
     >>> clims([0, 4], [1, 3])
     [0, 4]
@@ -302,7 +289,7 @@ def fnConstructor(s):
     if type(s) is str:
         s = s.split(':')[-1]
     else:
-        s = getHash(UTCDateTime())
+        s = get_hash(UTCDateTime())
 
     badchars = re.compile(r'[^A-Za-z0-9_. ]+|^\.|\.$|^ | $|^$')
     badsuffix = re.compile(r'(aux|com[1-9]|con|lpt[1-9]|prn)(\.|$)')
@@ -314,32 +301,75 @@ def fnConstructor(s):
     return fn
 
 
-def get_None(value):
+def get_none(value):
     """
     Convert "None" to None
     :param value:
-    :type value: str, bool
+    :type value: str, NoneType
     :return:
-    :rtype: bool
+    :rtype: type(value) or NoneType
+
+    >>> st = read()
+    >>> print(get_none(st))
+    3 Trace(s) in Stream:
+    BW.RJOB..EHZ | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
+    BW.RJOB..EHN | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
+    BW.RJOB..EHE | 2009-08-24T00:20:03.000000Z - 2009-08-24T00:20:32.990000Z | 100.0 Hz, 3000 samples
+    >>> get_none('Stream')
+    'Stream'
+    >>> get_none(0)
+    0
+    >>> get_none(0.)
+    0.0
+    >>> print(get_none('None'))
+    None
+    >>> print(get_none(None))
+    None
     """
-    if value == 'None':
+    if value is None or (type(value) is str and value == 'None'):
         return None
     else:
         return value
 
 
-def get_Bool(value):
+def get_bool(value):
     """
-    Convert string representations of bools to their true boolean value
+    Convert string representations of bools to their true boolean value. Return value if it cannot be identified as bool.
     :param value:
-    :type value: str, bool
+    :type value: str, bool, int, float
     :return: true boolean value
     :rtype: bool
+
+    >>> get_bool(True)
+    True
+    >>> get_bool(False)
+    False
+    >>> get_bool(0)
+    False
+    >>> get_bool(0.)
+    False
+    >>> get_bool(0.1)
+    True
+    >>> get_bool(2)
+    True
+    >>> get_bool(-1)
+    False
+    >>> get_bool(-0.3)
+    False
+    >>> get_bool(None)
+    None
     """
-    if value in ['True', 'true']:
+    if type(value) is bool:
+        return value
+    elif value in ['True', 'true']:
         return True
     elif value in ['False', 'false']:
         return False
+    elif isinstance(value, float) or isinstance(value, int):
+        if value > 0. or value > 0:
+            return True
+        else:
+            return False
     else:
         return value
 
@@ -353,8 +383,8 @@ def four_digits(year):
     :return: four digit year correspondent
     :rtype: int
 
-    >>> four_digits(20)
+    >>> four_digits(75)
-    1920
+    1975
     >>> four_digits(16)
     2016
     >>> four_digits(00)
@@ -436,36 +466,53 @@ def backtransformFilterString(st):
     return st
 
 
-def getHash(time):
+def get_hash(time):
     """
     takes a time object and returns the corresponding SHA1 hash of the formatted date string
     :param time: time object for which a hash should be calculated
     :type time: `~obspy.core.utcdatetime.UTCDateTime`
     :return: SHA1 hash
     :rtype: str
+
+    >>> time = UTCDateTime(0)
+    >>> get_hash(time)
+    '7627cce3b1b58dd21b005dac008b34d18317dd15'
+    >>> get_hash(0)
+    Traceback (most recent call last):
+    ...
+    AssertionError: 'time' is not an ObsPy UTCDateTime object
     """
+    assert isinstance(time, UTCDateTime), '\'time\' is not an ObsPy UTCDateTime object'
     hg = hashlib.sha1()
-    hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f'))
+    hg.update(time.strftime('%Y-%m-%d %H:%M:%S.%f').encode('utf-8'))
     return hg.hexdigest()
 
 
-def getLogin():
+def get_login():
     """
-    returns the actual user's login ID
+    returns the actual user's name
-    :return: login ID
+    :return: login name
     :rtype: str
     """
     import getpass
     return getpass.getuser()
 
 
-def getOwner(fn):
+def get_owner(fn):
     """
     takes a filename and return the login ID of the actual owner of the file
     :param fn: filename of the file tested
     :type fn: str
     :return: login ID of the file's owner
     :rtype: str
+
+    >>> import tempfile
+    >>> with tempfile.NamedTemporaryFile() as tmpfile:
+    ...     tmpfile.write(b'') and True
+    ...     tmpfile.flush()
+    ...     get_owner(tmpfile.name) == os.path.expanduser('~').split('/')[-1]
+    0
+    True
     """
     system_name = platform.system()
     if system_name in ["Linux", "Darwin"]:
@@ -511,6 +558,11 @@ def is_executable(fn):
     :param fn: path to the file to be tested
     :return: True or False
     :rtype: bool
+
+    >>> is_executable('/bin/ls')
+    True
+    >>> is_executable('/var/log/system.log')
+    False
     """
     return os.path.isfile(fn) and os.access(fn, os.X_OK)
 
@@ -537,24 +589,36 @@ def isSorted(iterable):
     >>> isSorted([2,3,1,4])
     False
     """
-    assert isIterable(iterable), 'object is not iterable; object: {' \
+    assert is_iterable(iterable), "object is not iterable; object: {}".format(iterable)
-                                 '}'.format(iterable)
     if type(iterable) is str:
         iterable = [s for s in iterable]
     return sorted(iterable) == iterable
 
 
-def isIterable(obj):
+def is_iterable(obj):
     """
     takes a python object and returns True is the object is iterable and
     False otherwise
     :param obj: a python object
-    :type obj: object
+    :type obj: obj
     :return: True of False
     :rtype: bool
+
+    >>> is_iterable(1)
+    False
+    >>> is_iterable(True)
+    False
+    >>> is_iterable(0.)
+    False
+    >>> is_iterable((0,1,3,4))
+    True
+    >>> is_iterable([1])
+    True
+    >>> is_iterable('a')
+    True
     """
     try:
-        iterator = iter(obj)
+        iter(obj)
     except TypeError as te:
         return False
     return True
@@ -563,13 +627,19 @@ def isIterable(obj):
 def key_for_set_value(d):
     """
     takes a dictionary and returns the first key for which's value the
-    boolean is True
+    boolean representation is True
     :param d: dictionary containing values
     :type d: dict
     :return: key to the first non-False value found; None if no value's
     boolean equals True
-    :rtype:
+    :rtype: bool or NoneType
+
+    >>> key_for_set_value({'one': 0, 'two': 1})
+    'two'
+    >>> print(key_for_set_value({1: 0, 2: False}))
+    None
     """
+    assert type(d) is dict, "Function only defined for inputs of type 'dict'."
     r = None
     for k, v in d.items():
         if v:
@@ -577,32 +647,53 @@ def key_for_set_value(d):
     return r
 
 
-def prepTimeAxis(stime, trace, verbosity=0):
+def prep_time_axis(offset, trace, verbosity=0):
     """
-    takes a starttime and a trace object and returns a valid time axis for
+    takes an offset and a trace object and returns a valid time axis for
     plotting
-    :param stime: start time of the actual seismogram as UTCDateTime
+    :param offset: offset of the actual seismogram on plotting axis
-    :type stime: `~obspy.core.utcdatetime.UTCDateTime`
+    :type offset: float or int
     :param trace: seismic trace object
     :type trace: `~obspy.core.trace.Trace`
     :param verbosity: if != 0, debug output will be written to console
     :type verbosity: int
     :return: valid numpy array with time stamps for plotting
     :rtype: `~numpy.ndarray`
+
+    >>> tr = read()[0]
+    >>> prep_time_axis(0., tr)
+    array([0.00000000e+00, 1.00033344e-02, 2.00066689e-02, ...,
+           2.99799933e+01, 2.99899967e+01, 3.00000000e+01])
+    >>> prep_time_axis(22.5, tr)
+    array([22.5       , 22.51000333, 22.52000667, ..., 52.47999333,
+           52.48999667, 52.5       ])
+    >>> prep_time_axis(tr.stats.starttime, tr)
+    Traceback (most recent call last):
+    ...
+    AssertionError: 'offset' is not of type 'float' or 'int'; type: <class 'obspy.core.utcdatetime.UTCDateTime'>
+    >>> tr.stats.npts -= 1
+    >>> prep_time_axis(0, tr)
+    array([0.00000000e+00, 1.00033356e-02, 2.00066711e-02, ...,
+           2.99699933e+01, 2.99799967e+01, 2.99900000e+01])
+    >>> tr.stats.npts += 2
+    >>> prep_time_axis(0, tr)
+    array([0.00000000e+00, 1.00033333e-02, 2.00066667e-02, ...,
+           2.99899933e+01, 2.99999967e+01, 3.00100000e+01])
     """
+    assert isinstance(offset, (float, int)), "'offset' is not of type 'float' or 'int'; type: {}".format(type(offset))
     nsamp = trace.stats.npts
     srate = trace.stats.sampling_rate
     tincr = trace.stats.delta
-    etime = stime + nsamp / srate
+    etime = offset + nsamp / srate
-    time_ax = np.linspace(stime, etime, nsamp)
+    time_ax = np.linspace(offset, etime, nsamp)
     if len(time_ax) < nsamp:
         if verbosity:
             print('elongate time axes by one datum')
-        time_ax = np.arange(stime, etime + tincr, tincr)
+        time_ax = np.arange(offset, etime + tincr, tincr)
     elif len(time_ax) > nsamp:
         if verbosity:
             print('shorten time axes by one datum')
-        time_ax = np.arange(stime, etime - tincr, tincr)
+        time_ax = np.arange(offset, etime - tincr, tincr)
     if len(time_ax) != nsamp:
         print('Station {0}, {1} samples of data \n '
               '{2} length of time vector \n'
@@ -618,13 +709,13 @@ def find_horizontals(data):
     :param data: waveform data
     :type data: `obspy.core.stream.Stream`
     :return: components list
-    :rtype: list
+    :rtype: List(str)
 
     ..example::
 
     >>> st = read()
     >>> find_horizontals(st)
-    [u'N', u'E']
+    ['N', 'E']
     """
     rval = []
     for tr in data:
@@ -635,7 +726,7 @@ def find_horizontals(data):
     return rval
 
 
-def pick_color(picktype, phase, quality=0):
+def pick_color(picktype: Literal['manual', 'automatic'], phase: Literal['P', 'S'], quality: int = 0) -> Rgba:
     """
     Create pick color by modifying the base color by the quality.
 
@@ -648,7 +739,7 @@ def pick_color(picktype, phase, quality=0):
     :param quality: quality of pick. Decides the new intensity of the modifier color
     :type quality: int
     :return: tuple containing modified rgba color values
-    :rtype: (int, int, int, int)
+    :rtype: Rgba
     """
     min_quality = 3
     bpc = base_phase_colors(picktype, phase)  # returns dict like {'modifier': 'g', 'rgba': (0, 0, 255, 255)}
@@ -704,17 +795,17 @@ def pick_linestyle_plt(picktype, key):
     return linestyles[picktype][key]
 
 
-def modify_rgba(rgba, modifier, intensity):
+def modify_rgba(rgba: Rgba, modifier: Literal['r', 'g', 'b'], intensity: float) -> Rgba:
     """
     Modify rgba color by adding the given intensity to the modifier color
     :param rgba: tuple containing rgba values
-    :type rgba: (int, int, int, int)
+    :type rgba: Rgba
-    :param modifier: which color should be modified, eg. 'r', 'g', 'b'
+    :param modifier: which color should be modified; options: 'r', 'g', 'b'
-    :type modifier: str
+    :type modifier: Literal['r', 'g', 'b']
     :param intensity: intensity to be added to selected color
     :type intensity: float
     :return: tuple containing rgba values
-    :rtype: (int, int, int, int)
+    :rtype: Rgba
     """
     rgba = list(rgba)
     index = {'r': 0,
@@ -748,18 +839,20 @@ def transform_colors_mpl_str(colors, no_alpha=False):
     Transforms rgba color values to a matplotlib string of color values with a range of [0, 1]
     :param colors: tuple of rgba color values ranging from [0, 255]
     :type colors: (float, float, float, float)
-    :param no_alpha: Wether to return a alpha value in the matplotlib color string
+    :param no_alpha: Whether to return an alpha value in the matplotlib color string
     :type no_alpha: bool
     :return: String containing r, g, b values and alpha value if no_alpha is False (default)
     :rtype: str
+
+    >>> transform_colors_mpl_str((255., 255., 255., 255.), True)
+    '(1.0, 1.0, 1.0)'
+    >>> transform_colors_mpl_str((255., 255., 255., 255.))
+    '(1.0, 1.0, 1.0, 1.0)'
     """
-    colors = list(colors)
-    colors_mpl = tuple([color / 255. for color in colors])
     if no_alpha:
-        colors_mpl = '({}, {}, {})'.format(*colors_mpl)
+        return '({}, {}, {})'.format(*transform_colors_mpl(colors))
     else:
-        colors_mpl = '({}, {}, {}, {})'.format(*colors_mpl)
+        return '({}, {}, {}, {})'.format(*transform_colors_mpl(colors))
-    return colors_mpl
 
 
 def transform_colors_mpl(colors):
@@ -769,27 +862,16 @@ def transform_colors_mpl(colors):
     :type colors: (float, float, float, float)
     :return: tuple of rgba color values ranging from [0, 1]
     :rtype: (float, float, float, float)
+
+    >>> transform_colors_mpl((127.5, 0., 63.75, 255.))
+    (0.5, 0.0, 0.25, 1.0)
+    >>> transform_colors_mpl(())
     """
     colors = list(colors)
     colors_mpl = tuple([color / 255. for color in colors])
     return colors_mpl
 
 
-def remove_underscores(data):
-    """
-    takes a `obspy.core.stream.Stream` object and removes all underscores
-    from station names
-    :param data: stream of seismic data
-    :type data: `~obspy.core.stream.Stream`
-    :return: data stream
-    :rtype: `~obspy.core.stream.Stream`
-    """
-    # for tr in data:
-    #    # remove underscores
-    #    tr.stats.station = tr.stats.station.strip('_')
-    return data
-
-
 def trim_station_components(data, trim_start=True, trim_end=True):
     """
     cut a stream so only the part common to all three traces is kept to avoid dealing with offsets
@@ -818,19 +900,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 +920,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 +936,7 @@ def check4gapsAndMerge(data):
         for merged_station in merged:
             print(merged_station)
 
-    return data
+    return data, gaps
 
 
 def check4doubled(data):
@@ -897,13 +966,53 @@ 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 = event.pylot_id #path.split('/')[-1]  # or event.pylot_id
+    filename = os.path.join(eventpath, 'PyLoT_' + eventname + fext)
+    if os.path.isfile(filename):
+        return filename
+
+
+def get_possible_pylot_eventfile_extensions(event, fext):
+    if hasattr(event, 'path'):
+        eventpath = event.path
+    else:
+        logging.warning('No attribute path found for event.')
+        return []
+    eventname = event.pylot_id
+    filename = os.path.join(eventpath, 'PyLoT_' + eventname + fext)
+    filenames = glob.glob(filename)
+    extensions = [os.path.split(path)[-1].split('PyLoT_' + eventname)[-1] for path in filenames]
+    return extensions
+
+
 def get_stations(data):
     """
-    Get list of all station names in data stream
+    Get list of all station names in data-stream
     :param data: stream containing seismic traces
     :type data: `~obspy.core.stream.Stream`
     :return: list of all station names in data, no duplicates
-    :rtype: list of str
+    :rtype: List(str)
     """
     stations = []
     for tr in data:
@@ -930,66 +1039,88 @@ def check4rotated(data, metadata=None, verbosity=1):
     :rtype: `~obspy.core.stream.Stream`
     """
 
-    def rotate_components(wfstream, metadata=None):
+    def rotation_required(trace_ids):
+        """
+        Derive if any rotation is required from the orientation code of the input.
+
+        :param trace_ids: string identifier of waveform data trace
+        :type trace_ids: List(str)
+        :return: boolean representing if rotation is necessary for any of the traces
+        :rtype: bool
+        """
+        orientations = [trace_id[-1] for trace_id in trace_ids]
+        return any([orientation.isnumeric() for orientation in orientations])
+
+    def rotate_components(wfs_in, metadata=None):
         """
         Rotate components if orientation code is numeric (= non traditional orientation).
 
         Azimut and dip are fetched from metadata. To be rotated, traces of a station have to be cut to the same length.
         Returns unrotated traces of no metadata is provided
-        :param wfstream: stream containing seismic traces of a station
+        :param wfs_in: stream containing seismic traces of a station
-        :type wfstream: `~obspy.core.stream.Stream`
+        :type wfs_in: `~obspy.core.stream.Stream`
         :param metadata: tuple containing metadata type string and metadata parser object
         :type metadata: (str, `~obspy.io.xseed.parser.Parser`)
         :return: stream object with traditionally oriented traces (ZNE)
         :rtype: `~obspy.core.stream.Stream`
         """
 
-        # check if any traces in this station need to be rotated
+        if len(wfs_in) < 3:
-        trace_ids = [trace.id for trace in wfstream]
+            print(f"Stream {wfs_in=}, has not enough components to rotate.")
-        orientations = [trace_id[-1] for trace_id in trace_ids]
+            return wfs_in
-        rotation_required = [orientation.isnumeric() for orientation in orientations]
-        if any(rotation_required):
-            t_start = full_range(wfstream)
-            try:
-                azimuts = []
-                dips = []
-                for tr_id in trace_ids:
-                    azimuts.append(metadata.get_coordinates(tr_id, t_start)['azimuth'])
-                    dips.append(metadata.get_coordinates(tr_id, t_start)['dip'])
-            except (KeyError, TypeError) as e:
-                print('Failed to rotate trace {}, no azimuth or dip available in metadata'.format(tr_id))
-                return wfstream
-            if len(wfstream) < 3:
-                print('Failed to rotate Stream {}, not enough components available.'.format(wfstream))
-                return wfstream
-            # to rotate all traces must have same length, so trim them
-            wfstream = trim_station_components(wfstream, trim_start=True, trim_end=True)
-            try:
-                z, n, e = rotate2zne(wfstream[0], azimuts[0], dips[0],
-                                     wfstream[1], azimuts[1], dips[1],
-                                     wfstream[2], azimuts[2], dips[2])
-                print('check4rotated: rotated trace {} to ZNE'.format(trace_ids))
-                # replace old data with rotated data, change the channel code to ZNE
-                z_index = dips.index(min(
-                    dips))  # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90 being vertical)
-                wfstream[z_index].data = z
-                wfstream[z_index].stats.channel = wfstream[z_index].stats.channel[0:-1] + 'Z'
-                del trace_ids[z_index]
-                for trace_id in trace_ids:
-                    coordinates = metadata.get_coordinates(trace_id, t_start)
-                    dip, az = coordinates['dip'], coordinates['azimuth']
-                    trace = wfstream.select(id=trace_id)[0]
-                    if az > 315 or az <= 45 or az > 135 and az <= 225:
-                        trace.data = n
-                        trace.stats.channel = trace.stats.channel[0:-1] + 'N'
-                    elif az > 45 and az <= 135 or az > 225 and az <= 315:
-                        trace.data = e
-                        trace.stats.channel = trace.stats.channel[0:-1] + 'E'
-            except (ValueError) as e:
-                print(e)
-                return wfstream
 
-        return wfstream
+        # check if any traces in this station need to be rotated
+        trace_ids = [trace.id for trace in wfs_in]
+        if not rotation_required(trace_ids):
+            logging.debug(f"Stream does not need any rotation: Traces are {trace_ids=}")
+            return wfs_in
+
+        # check metadata quality
+        t_start = full_range(wfs_in)[0]
+        try:
+            azimuths = []
+            dips = []
+            for tr_id in trace_ids:
+                azimuths.append(metadata.get_coordinates(tr_id, t_start)['azimuth'])
+                dips.append(metadata.get_coordinates(tr_id, t_start)['dip'])
+        except (KeyError, TypeError) as err:
+            logging.warning(f"Rotating not possible, not all azimuth and dip information "
+                            f"available in metadata. Stream remains unchanged.")
+            logging.debug(f"Rotating not possible, {err=}, {type(err)=}")
+            return wfs_in
+        except Exception as err:
+            print(f"Unexpected {err=}, {type(err)=}")
+            raise
+
+        # to rotate all traces must have same length, so trim them
+        wfs_out = trim_station_components(wfs_in, trim_start=True, trim_end=True)
+        try:
+            z, n, e = rotate2zne(wfs_out[0], azimuths[0], dips[0],
+                                 wfs_out[1], azimuths[1], dips[1],
+                                 wfs_out[2], azimuths[2], dips[2])
+            print('check4rotated: rotated trace {} to ZNE'.format(trace_ids))
+            # replace old data with rotated data, change the channel code to ZNE
+            z_index = dips.index(min(
+                dips))  # get z-trace index, z has minimum dip of -90 (dip is measured from 0 to -90, with -90
+            # being vertical)
+            wfs_out[z_index].data = z
+            wfs_out[z_index].stats.channel = wfs_out[z_index].stats.channel[0:-1] + 'Z'
+            del trace_ids[z_index]
+            for trace_id in trace_ids:
+                coordinates = metadata.get_coordinates(trace_id, t_start)
+                dip, az = coordinates['dip'], coordinates['azimuth']
+                trace = wfs_out.select(id=trace_id)[0]
+                if az > 315 or az <= 45 or 135 < az <= 225:
+                    trace.data = n
+                    trace.stats.channel = trace.stats.channel[0:-1] + 'N'
+                elif 45 < az <= 135 or 225 < az <= 315:
+                    trace.data = e
+                    trace.stats.channel = trace.stats.channel[0:-1] + 'E'
+        except ValueError as err:
+            print(f"{err=} Rotation failed. Stream remains unchanged.")
+            return wfs_in
+
+        return wfs_out
 
     if metadata is None:
         if verbosity:
@@ -1003,38 +1134,6 @@ def check4rotated(data, metadata=None, verbosity=1):
     return data
 
 
-def scaleWFData(data, factor=None, components='all'):
-    """
-    produce scaled waveforms from given waveform data and a scaling factor,
-    waveform may be selected by their components name
-    :param data: waveform data to be scaled
-    :type data: `~obspy.core.stream.Stream` object
-    :param factor: scaling factor
-    :type factor: float
-    :param components: components labels for the traces in data to be scaled by
-     the scaling factor (optional, default: 'all')
-    :type components: tuple
-    :return:  scaled waveform data
-    :rtype: `~obspy.core.stream.Stream` object
-    """
-    if components != 'all':
-        for comp in components:
-            if factor is None:
-                max_val = np.max(np.abs(data.select(component=comp)[0].data))
-                data.select(component=comp)[0].data /= 2 * max_val
-            else:
-                data.select(component=comp)[0].data /= 2 * factor
-    else:
-        for tr in data:
-            if factor is None:
-                max_val = float(np.max(np.abs(tr.data)))
-                tr.data /= 2 * max_val
-            else:
-                tr.data /= 2 * factor
-
-    return data
-
-
 def runProgram(cmd, parameter=None):
     """
     run an external program specified by cmd with parameters input returning the
@@ -1106,6 +1205,7 @@ def identifyPhase(phase):
         return False
 
 
+@lru_cache
 def identifyPhaseID(phase):
     """
     Returns phase id (capital P or S)
@@ -1169,7 +1269,7 @@ def correct_iplot(iplot):
     try:
         iplot = int(iplot)
     except ValueError:
-        if get_Bool(iplot):
+        if get_bool(iplot):
             iplot = 2
         else:
             iplot = 0

pylot/correlation/__init__.py

@@ -0,0 +1,2 @@
+# -*- coding: utf-8 -*-
+#

pylot/correlation/parameters_adriaarray.yaml

@@ -0,0 +1,90 @@
+#############################  correlation parameters  #####################################
+# min_corr_stacking: minimum correlation coefficient for building beam trace
+# min_corr_export: minimum correlation coefficient for pick export
+# min_stack: minimum number of stations for building beam trace
+# t_before: correlation window before pick
+# t_after: correlation window after pick#
+# cc_maxlag: maximum shift for initial correlation
+# cc_maxlag2: maximum shift for second (final) correlation (also for calculating pick uncertainty)
+# initial_pick_outlier_threshold: (hopefully) threshold for excluding large outliers of initial (AIC) picks
+# export_threshold: automatically exclude all onsets which deviate more than this threshold from corrected taup onsets
+# min_picks_export: minimum number of correlated picks for export
+# min_picks_autopylot: minimum number of reference autopicks picks to continue with event
+# check_RMS: do RMS check to search for restitution errors (very experimental)
+# use_taupy_onsets: use taupy onsets as reference picks instead of external picks
+# station_list: use the following stations as reference for stacking
+# use_stacked_trace: use existing stacked trace if found (spare re-computation)
+# data_dir: obspyDMT data subdirectory (e.g. 'raw', 'processed')
+# pickfile_extension: use quakeML files (PyLoT output) with the following extension, e.g. '_autopylot' for pickfiles
+#                     such as 'PyLoT_20170501_141822_autopylot.xml'
+
+logging: info
+pick_phases: ['P', 'S']
+
+# P-phase
+P:
+  min_corr_stacking: 0.8
+  min_corr_export: 0.6
+  min_stack: 20
+  t_before: 30.
+  t_after: 50.
+  cc_maxlag: 50.
+  cc_maxlag2: 5.
+  initial_pick_outlier_threshold: 30.
+  export_threshold: 2.5
+  min_picks_export: 100
+  min_picks_autopylot: 50
+  check_RMS: True
+  use_taupy_onsets: False
+  station_list: ['HU.MORH', 'HU.TIH',  'OX.FUSE', 'OX.BAD']
+  use_stacked_trace: False
+  data_dir: 'processed'
+  pickfile_extension: '_autopylot'
+  dt_stacking: [250, 250]
+
+  # filter for first correlation (rough)
+  filter_options:
+    freqmax: 0.5
+    freqmin: 0.03
+  # filter for second correlation (fine)
+  filter_options_final:
+    freqmax: 0.5
+    freqmin: 0.03
+
+  filter_type: bandpass
+  sampfreq: 20.0
+
+# S-phase
+S:
+  min_corr_stacking: 0.7
+  min_corr_export: 0.6
+  min_stack: 20
+  t_before: 60.
+  t_after: 60.
+  cc_maxlag: 100.
+  cc_maxlag2: 25.
+  initial_pick_outlier_threshold: 30.
+  export_threshold: 5.0
+  min_picks_export: 200
+  min_picks_autopylot: 50
+  check_RMS: True
+  use_taupy_onsets: False
+  station_list: ['HU.MORH','HU.TIH',  'OX.FUSE', 'OX.BAD']
+  use_stacked_trace: False
+  data_dir: 'processed'
+  pickfile_extension: '_autopylot'
+  dt_stacking: [250, 250]
+
+  # filter for first correlation (rough)
+  filter_options:
+    freqmax: 0.1
+    freqmin: 0.01
+
+  # filter for second correlation (fine)
+  filter_options_final:
+    freqmax: 0.2
+    freqmin: 0.01
+
+  filter_type: bandpass
+  sampfreq: 20.0
+

pylot/correlation/submit_pick_corr_correction.sh

@@ -0,0 +1,40 @@
+#!/bin/bash
+
+#ulimit -s 8192
+#ulimit -v $(ulimit -v | awk '{printf("%d",$1*0.95)}')
+#ulimit -v
+
+#655360
+
+source /opt/anaconda3/etc/profile.d/conda.sh
+conda activate pylot_311
+NSLOTS=20
+
+#qsub -l low -cwd -l "os=*stretch" -pe smp 40 submit_pick_corr_correction.sh
+#$ -l low
+#$ -l h_vmem=6G
+#$ -cwd
+#$ -pe smp 20
+#$ -N corr_pick
+
+
+export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot_tools/"
+export PYTHONPATH="$PYTHONPATH:/home/marcel/git/"
+export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot/"
+
+#export MKL_NUM_THREADS=${NSLOTS:=1}
+#export NUMEXPR_NUM_THREADS=${NSLOTS:=1}
+#export OMP_NUM_THREADS=${NSLOTS:=1}
+
+#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 0 -istop 100
+#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 100 -istop 200
+#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M6.0-6.5' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 0 -istop 100
+#python pick_correlation_correction.py '/data/AlpArray_Data/dmt_database_mantle_M5.8-6.0' '/home/marcel/.pylot/pylot_alparray_mantle_corr_stack_0.03-0.5.in' -pd -n ${NSLOTS:=1} -istart 100 -istop 200
+#python pick_correlation_correction.py 'H:\sciebo\dmt_database' 'H:\Sciebo\dmt_database\pylot_alparray_mantle_corr_S_0.01-0.2.in' -pd -n 4 -t
+
+pylot_infile='/home/marcel/.pylot/pylot_alparray_syn_fwi_mk6_it3.in'
+#pylot_infile='/home/marcel/.pylot/pylot_adriaarray_corr_P_and_S.in'
+
+# THIS SCRIPT SHOLD BE CALLED BY "submit_to_grid_engine.py" using the following line:
+python pick_correlation_correction.py $1 $pylot_infile -pd -n ${NSLOTS:=1} -istart $2 --params 'parameters_fwi_mk6_it3.yaml'
+#--event_blacklist eventlist.txt

pylot/correlation/submit_to_grid_engine.py

@@ -0,0 +1,23 @@
+#!/usr/bin/env python
+
+import subprocess
+
+fnames = [
+    ('/data/AlpArray_Data/dmt_database_synth_model_mk6_it3_no_rotation', 0),
+          ]
+
+#fnames = [('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_SKS-phase', 0),
+#          ('/data/AlpArray_Data/dmt_database_mantle_0.01-0.2_S-phase', 0),]
+
+####
+script_location = '/home/marcel/VersionCtrl/git/code_base/correlation_picker/submit_pick_corr_correction.sh'
+####
+
+for fnin, istart in fnames:
+    input_cmds = f'qsub -q low.q@minos15,low.q@minos14,low.q@minos13,low.q@minos12,low.q@minos11 {script_location} {fnin} {istart}'
+
+    print(input_cmds)
+    print(subprocess.check_output(input_cmds.split()))
+
+
+

pylot/correlation/utils.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+
+import os
+import glob
+import json
+
+from obspy import read_events
+
+from pylot.core.util.dataprocessing import Metadata
+from pylot.core.util.obspyDMT_interface import qml_from_obspyDMT
+
+
+def get_event_obspy_dmt(eventdir):
+    event_pkl_file = os.path.join(eventdir, 'info', 'event.pkl')
+    if not os.path.exists(event_pkl_file):
+        raise IOError('Could not find event path for event: {}'.format(eventdir))
+    event = qml_from_obspyDMT(event_pkl_file)
+    return event
+
+
+def get_event_pylot(eventdir, extension=''):
+    event_id = get_event_id(eventdir)
+    filename = os.path.join(eventdir, 'PyLoT_{}{}.xml'.format(event_id, extension))
+    if not os.path.isfile(filename):
+        return
+    cat = read_events(filename)
+    return cat[0]
+
+
+def get_event_id(eventdir):
+    event_id = os.path.split(eventdir)[-1]
+    return event_id
+
+
+def get_picks(eventdir, extension=''):
+    event_id = get_event_id(eventdir)
+    filename = 'PyLoT_{}{}.xml'
+    filename = filename.format(event_id, extension)
+    fpath = os.path.join(eventdir, filename)
+    fpaths = glob.glob(fpath)
+    if len(fpaths) == 1:
+        cat = read_events(fpaths[0])
+        picks = cat[0].picks
+        return picks
+    elif len(fpaths) == 0:
+        print('get_picks: File not found: {}'.format(fpath))
+        return
+    print(f'WARNING: Ambiguous pick file specification. Found the following pick files {fpaths}\nFilemask: {fpath}')
+    return
+
+
+def write_json(object, fname):
+    with open(fname, 'w') as outfile:
+        json.dump(object, outfile, sort_keys=True, indent=4)
+
+
+def get_metadata(eventdir):
+    metadata_path = os.path.join(eventdir, 'resp')
+    metadata = Metadata(inventory=metadata_path, verbosity=0)
+    return metadata

requirements.txt

@@ -1,12 +1,7 @@
-# This file may be used to create an environment using:
+Cartopy==0.23.0
-# $ conda create --name <env> --file <this file>
+joblib==1.4.2
-# platform: win-64
+obspy==1.4.1
-cartopy=0.20.2
+pyaml==24.7.0
-matplotlib-base=3.3.4
+pyqtgraph==0.13.7
-numpy=1.22.3
+PySide2==5.15.8
-obspy=1.3.0
+pytest==8.3.2
-pyqtgraph=0.12.4
-pyside2=5.13.2
-python=3.8.12
-qt=5.12.9
-scipy=1.8.0

tests/test_autopicker/pylot_alparray_mantle_corr_stack_0.03-0.5.in

@@ -0,0 +1,99 @@
+%This is a parameter input file for PyLoT/autoPyLoT.
+%All main and special settings regarding data handling
+%and picking are to be set here!
+%Parameters are optimized for %extent data sets!
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#main settings#
+dmt_database_test                                  #datapath#       %data path
+20171010_063224.a                                  #eventID#        %event ID for single event processing (* for all events found in database)
+                                                   #invdir#         %full path to  inventory or dataless-seed file                                                       
+PILOT                                              #datastructure#  %choose data structure                                                                               
+True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#NLLoc settings#
+None                                               #nllocbin#       %path to NLLoc executable                                                                            
+None                                               #nllocroot#      %root of NLLoc-processing directory                                                                  
+None                                               #phasefile#      %name of autoPyLoT-output phase file for NLLoc                                                       
+None                                               #ctrfile#        %name of autoPyLoT-output control file for NLLoc                                                     
+ttime                                              #ttpatter#       %pattern of NLLoc ttimes from grid                                                                   
+AUTOLOC_nlloc                                      #outpatter#      %pattern of NLLoc-output file                                                                        
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#parameters for seismic moment estimation#
+3530.0                                             #vp#             %average P-wave velocity                                                                             
+2500.0                                             #rho#            %average rock density [kg/m^3]                                                                       
+300.0 0.8                                          #Qp#             %quality factor for P waves (Qp*f^a); list(Qp, a)                                                    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#settings local magnitude#
+1.0 1.0 1.0                                        #WAscaling#      %Scaling relation (log(Ao)+Alog(r)+Br+C) of Wood-Anderson amplitude Ao [nm]                           If zeros are set, original Richter magnitude is calculated!
+1.0 1.0                                            #magscaling#     %Scaling relation for derived local magnitude [a*Ml+b].                            If zeros are set, no scaling of network magnitude is applied!
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#filter settings#
+0.03 0.03                                          #minfreq#        %Lower filter frequency [P, S]                                                                       
+0.5 0.5                                            #maxfreq#        %Upper filter frequency [P, S]                                                                       
+4 4                                                #filter_order#   %filter order [P, S]                                                                                 
+bandpass bandpass                                  #filter_type#    %filter type (bandpass, bandstop, lowpass, highpass) [P, S]                                          
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+#common settings picker#
+global                                             #extent#         %extent of array ("local", "regional" or "global")                                                   
+-100.0                                             #pstart#         %start time [s] for calculating CF for P-picking (if TauPy: seconds relative to estimated onset)     
+50.0                                               #pstop#          %end time [s] for calculating CF for P-picking (if TauPy: seconds relative to estimated onset)       
+-50.0                                              #sstart#         %start time [s] relative to P-onset for calculating CF for S-picking                                 
+50.0                                               #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
+True                                               #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
+ak135                                              #taup_model#     %Define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
+P,Pdiff,S,SKS                                      #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
+0.03 0.5                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
+0.01 0.5                                           #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
+0.03 0.5                                           #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     
+0.01 0.5                                           #bph2#           %lower/upper corner freq. of second band pass filter z-comp. [Hz]                                    
+#special settings for calculating CF#
+%!!Edit the following only if you know what you are doing!!%
+#Z-component#
+HOS                                                #algoP#          %choose algorithm for P-onset determination (HOS, ARZ, or AR3)                                       
+300.0                                              #tlta#           %for HOS-/AR-AIC-picker, length of LTA window [s]                                                    
+4                                                  #hosorder#       %for HOS-picker, order of Higher Order Statistics                                                    
+2                                                  #Parorder#       %for AR-picker, order of AR process of Z-component                                                   
+16.0                                               #tdet1z#         %for AR-picker, length of AR determination window [s] for Z-component, 1st pick                      
+10.0                                               #tpred1z#        %for AR-picker, length of AR prediction window [s] for Z-component, 1st pick                         
+12.0                                               #tdet2z#         %for AR-picker, length of AR determination window [s] for Z-component, 2nd pick                      
+6.0                                                #tpred2z#        %for AR-picker, length of AR prediction window [s] for Z-component, 2nd pick                         
+0.001                                              #addnoise#       %add noise to seismogram for stable AR prediction                                                    
+60.0 5.0 20.0 12.0                                 #tsnrz#          %for HOS/AR, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]     
+50.0                                               #pickwinP#       %for initial AIC pick, length of P-pick window [s]                                                   
+30.0                                               #Precalcwin#     %for HOS/AR, window length [s] for recalculation of CF (relative to 1st pick)                        
+2.0                                                #aictsmooth#     %for HOS/AR, take average of samples for smoothing of AIC-function [s]                               
+2.0                                                #tsmoothP#       %for HOS/AR, take average of samples in this time window for smoothing CF [s]                        
+0.006                                              #ausP#           %for HOS/AR, artificial uplift of samples (aus) of CF (P)                                            
+2.0                                                #nfacP#          %for HOS/AR, noise factor for noise level determination (P)                                          
+#H-components#
+ARH                                                #algoS#          %choose algorithm for S-onset determination (ARH or AR3)                                             
+12.0                                               #tdet1h#         %for HOS/AR, length of AR-determination window [s], H-components, 1st pick                           
+6.0                                                #tpred1h#        %for HOS/AR, length of AR-prediction window [s], H-components, 1st pick                              
+8.0                                                #tdet2h#         %for HOS/AR, length of AR-determinaton window [s], H-components, 2nd pick                            
+4.0                                                #tpred2h#        %for HOS/AR, length of AR-prediction window [s], H-components, 2nd pick                              
+4                                                  #Sarorder#       %for AR-picker, order of AR process of H-components                                                  
+100.0                                              #Srecalcwin#     %for AR-picker, window length [s] for recalculation of CF (2nd pick) (H)                             
+195.0                                              #pickwinS#       %for initial AIC pick, length of S-pick window [s]                                                   
+60.0 10.0 30.0 12.0                                #tsnrh#          %for ARH/AR3, window lengths for SNR-and slope estimation [tnoise, tsafetey, tsignal, tslope] [s]    
+22.0                                               #aictsmoothS#    %for AIC-picker, take average of samples in this time window for smoothing of AIC-function [s]       
+20.0                                               #tsmoothS#       %for AR-picker, take average of samples for smoothing CF [s] (S)                                     
+0.001                                              #ausS#           %for HOS/AR, artificial uplift of samples (aus) of CF (S)                                            
+2.0                                                #nfacS#          %for AR-picker, noise factor for noise level determination (S)                                       
+#first-motion picker#
+1                                                  #minfmweight#    %minimum required P weight for first-motion determination                                            
+3.0                                                #minFMSNR#       %miniumum required SNR for first-motion determination                                                
+10.0                                               #fmpickwin#      %pick window [s] around P onset for calculating zero crossings                                       
+#quality assessment#
+0.1 0.2 0.4 0.8                                    #timeerrorsP#    %discrete time errors [s] corresponding to picking weights [0 1 2 3] for P                           
+4.0 8.0 16.0 32.0                                  #timeerrorsS#    %discrete time errors [s] corresponding to picking weights [0 1 2 3] for S                           
+0.005                                              #minAICPslope#   %below this slope [counts/s] the initial P pick is rejected                                          
+1.1                                                #minAICPSNR#     %below this SNR the initial P pick is rejected                                                       
+0.002                                              #minAICSslope#   %below this slope [counts/s] the initial S pick is rejected                                          
+1.3                                                #minAICSSNR#     %below this SNR the initial S pick is rejected                                                       
+20.0                                               #minsiglength#   %length of signal part for which amplitudes must exceed noiselevel [s]                               
+1.0                                                #noisefactor#    %noiselevel*noisefactor=threshold                                                                    
+10.0                                               #minpercent#     %required percentage of amplitudes exceeding threshold                                               
+0.1                                                #zfac#           %P-amplitude must exceed at least zfac times RMS-S amplitude                                         
+100.0                                              #mdttolerance#   %maximum allowed deviation of P picks from median [s]                                                
+50.0                                               #wdttolerance#   %maximum allowed deviation from Wadati-diagram                                                       
+25.0                                               #jackfactor#     %pick is removed if the variance of the subgroup with the pick removed is larger than the mean variance of all subgroups times safety factor

tests/test_autopicker/test_autopylot.py

@@ -0,0 +1,67 @@
+import os
+import pytest
+
+from obspy import read_events
+
+from autoPyLoT import autoPyLoT
+
+
+class TestAutopickerGlobal():
+    def init(self):
+        self.params_infile = 'pylot_alparray_mantle_corr_stack_0.03-0.5.in'
+        self.test_event_dir = 'dmt_database_test'
+        self.fname_outfile_xml = os.path.join(
+            self.test_event_dir, '20171010_063224.a', 'PyLoT_20171010_063224.a_autopylot.xml'
+        )
+
+        # check if the input files exist
+        if not os.path.isfile(self.params_infile):
+            print(f'Test input file {os.path.abspath(self.params_infile)} not found.')
+            return False
+
+        if not os.path.exists(self.test_event_dir):
+            print(
+                f'Test event directory not found at location "{os.path.abspath(self.test_event_dir)}". '
+                f'Make sure to load it from the website first.'
+            )
+            return False
+
+        return True
+
+    def test_autopicker(self):
+        assert self.init(), 'Initialization failed due to missing input files.'
+        # check for output file in test directory and remove it if necessary
+        if os.path.isfile(self.fname_outfile_xml):
+            os.remove(self.fname_outfile_xml)
+        autoPyLoT(inputfile=self.params_infile, eventid='20171010_063224.a', obspyDMT_wfpath='processed')
+
+        # test for different known output files if they are identical or not
+        compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_autopylot.xml', True)
+        compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_saved_from_GUI.xml', True)
+        compare_pickfiles(self.fname_outfile_xml, 'PyLoT_20171010_063224.a_corrected_taup_times_0.03-0.5_P.xml', False)
+
+
+def compare_pickfiles(pickfile1: str, pickfile2: str, samefile: bool = True) -> None:
+    """
+    Compare the pick times and errors from two pick files.
+
+    Parameters:
+        pickfile1 (str): The path to the first pick file.
+        pickfile2 (str): The path to the second pick file.
+        samefile (bool): A flag indicating whether the two files are expected to be the same. Defaults to True.
+
+    Returns:
+        None
+    """
+    cat1 = read_events(pickfile1)
+    cat2 = read_events(pickfile2)
+    picks1 = sorted(cat1[0].picks, key=lambda pick: str(pick.waveform_id))
+    picks2 = sorted(cat2[0].picks, key=lambda pick: str(pick.waveform_id))
+    pick_times1 = [pick.time for pick in picks1]
+    pick_times2 = [pick.time for pick in picks2]
+    pick_terrs1 = [pick.time_errors for pick in picks1]
+    pick_terrs2 = [pick.time_errors for pick in picks2]
+
+    # check if times and errors are identical or not depending on the samefile flag
+    assert (pick_times1 == pick_times2) is samefile, 'Pick times error'
+    assert (pick_terrs1 == pick_terrs2) is samefile, 'Pick time errors errors'

tests/test_autopickstation/autoPyLoT_global_taupy_false.in

@@ -4,10 +4,8 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-/home/darius                                       #rootpath#       %project path                                                                                        
+/home/darius/alparray/waveforms_used               #datapath#       %data path
-alparray                                           #datapath#       %data path                                                                                           
+e0093.173.16                                       #eventID#        %event ID for single event processing (* for all events found in datapath)
-waveforms_used                                     #database#       %name of data base                                                                                   
-e0093.173.16                                       #eventID#        %event ID for single event processing (* for all events found in database)                           
 /home/darius/alparray/metadata                     #invdir#         %full path to  inventory or dataless-seed file                                                       
 PILOT                                              #datastructure#  %choose data structure                                                                               
 True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
@@ -43,6 +41,7 @@ global                                             #extent#         %extent of a
 875.0                                              #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
 False                                              #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
 IASP91                                             #taup_model#     %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
+P,Pdiff,S,Sdiff                                    #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
 0.01 0.1                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
 0.001 0.5                                          #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
 0.01 0.5                                           #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     

tests/test_autopickstation/autoPyLoT_global_taupy_true.in

@@ -4,10 +4,8 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-/home/darius                                       #rootpath#       %project path                                                                                        
+/home/darius/alparray/waveforms_used               #datapath#       %data path
-alparray                                           #datapath#       %data path                                                                                           
+e0093.173.16                                       #eventID#        %event ID for single event processing (* for all events found in datapath)
-waveforms_used                                     #database#       %name of data base                                                                                   
-e0093.173.16                                       #eventID#        %event ID for single event processing (* for all events found in database)                           
 /home/darius/alparray/metadata                     #invdir#         %full path to  inventory or dataless-seed file                                                       
 PILOT                                              #datastructure#  %choose data structure                                                                               
 True                                               #apverbose#      %choose 'True' or 'False' for terminal output                                                        
@@ -43,6 +41,7 @@ global                                             #extent#         %extent of a
 875.0                                              #sstop#          %end time [s] after P-onset for calculating CF for S-picking                                         
 True                                               #use_taup#       %use estimated traveltimes from TauPy for calculating windows for CF                                 
 IASP91                                             #taup_model#     %define TauPy model for traveltime estimation. Possible values: 1066a, 1066b, ak135, ak135f, herrin, iasp91, jb, prem, pwdk, sp6
+P,Pdiff,S,Sdiff                                    #taup_phases#    %Specify possible phases for TauPy (comma separated). See Obspy TauPy documentation for possible values.
 0.01 0.1                                           #bpz1#           %lower/upper corner freq. of first band pass filter Z-comp. [Hz]                                     
 0.001 0.5                                          #bpz2#           %lower/upper corner freq. of second band pass filter Z-comp. [Hz]                                    
 0.01 0.5                                           #bph1#           %lower/upper corner freq. of first band pass filter H-comp. [Hz]                                     

tests/test_autopickstation/test_autopickstation.py

@@ -1,6 +1,7 @@
 import os
 import sys
 import unittest
+import pytest
 
 import obspy
 from obspy import UTCDateTime
@@ -24,7 +25,7 @@ class HidePrints:
         if self.hide:
             self._original_stdout = sys.stdout
             devnull = open(os.devnull, "w")
-            sys.stdout = devnull
+            #sys.stdout = devnull
 
     def __exit__(self, exc_type, exc_val, exc_tb):
         if self.hide:
@@ -105,7 +106,6 @@ class TestAutopickStation(unittest.TestCase):
         # show complete diff when difference in results dictionaries are found
         self.maxDiff = None
 
-    # @skip("Works")
     def test_autopickstation_taupy_disabled_gra1(self):
         expected = {
             'P': {'picker': 'auto', 'snrdb': 15.405649120980094, 'weight': 0, 'Mo': None, 'marked': [], 'Mw': None,
@@ -121,8 +121,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_disabled,
                                               metadata=(None, None))
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA1', station)
 
     def test_autopickstation_taupy_enabled_gra1(self):
@@ -140,8 +140,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.gra1, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA1', station)
 
     def test_autopickstation_taupy_disabled_gra2(self):
@@ -157,8 +157,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_disabled,
                                               metadata=(None, None))
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA2', station)
 
     def test_autopickstation_taupy_enabled_gra2(self):
@@ -175,8 +175,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.gra2, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA2', station)
 
     def test_autopickstation_taupy_disabled_ech(self):
@@ -190,8 +190,8 @@ class TestAutopickStation(unittest.TestCase):
                           'fm': None, 'spe': None, 'channel': u'LHE'}}
         with HidePrints():
             result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_disabled)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('ECH', station)
 
     def test_autopickstation_taupy_enabled_ech(self):
@@ -208,8 +208,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.ech, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('ECH', station)
 
     def test_autopickstation_taupy_disabled_fiesa(self):
@@ -224,8 +224,8 @@ class TestAutopickStation(unittest.TestCase):
                           'fm': None, 'spe': None, 'channel': u'LHE'}}
         with HidePrints():
             result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_disabled)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('FIESA', station)
 
     def test_autopickstation_taupy_enabled_fiesa(self):
@@ -242,8 +242,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.fiesa, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertDictContainsSubset(expected=expected['P'], actual=result['P'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('FIESA', station)
 
     def test_autopickstation_gra1_z_comp_missing(self):
@@ -271,8 +271,9 @@ class TestAutopickStation(unittest.TestCase):
                   'fm': 'N', 'channel': None}}
         with HidePrints():
             result, station = autopickstation(wfstream=wfstream, pickparam=self.pickparam_taupy_disabled,
-                                              metadata=(None, None))
+                                              metadata=(None, None), iplot=2)
-        self.assertEqual(expected, result)
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA1', station)
 
     def test_autopickstation_a106_taupy_enabled(self):
@@ -290,7 +291,9 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.a106, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertEqual(expected, result)
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
+
 
     def test_autopickstation_station_missing_in_metadata(self):
         """This station is not in the metadata, but Taupy is enabled. Taupy should exit cleanly and modify the starttime
@@ -311,8 +314,37 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=self.a005a, pickparam=self.pickparam_taupy_enabled,
                                               metadata=self.metadata, origin=self.origin)
-        self.assertEqual(expected, result)
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
 
 
+def run_dict_comparison(result, expected):
+    for key, expected_value in expected.items():
+        if isinstance(expected_value, dict):
+            run_dict_comparison(result[key], expected[key])
+        else:
+            res = result[key]
+            if isinstance(res, UTCDateTime) and isinstance(expected_value, UTCDateTime):
+                res = res.timestamp
+                expected_value = expected_value.timestamp
+            assert expected_value == pytest.approx(res), f'{key}: {expected_value} != {res}'
+
+
+def compare_dicts(result, expected, hint=''):
+    try:
+        run_dict_comparison(result, expected)
+    except AssertionError:
+        raise AssertionError(f'{hint}Dictionaries not equal.'
+                             f'\n\n<<Expected>>\n{pretty_print_dict(expected)}'
+                             f'\n\n<<Result>>\n{pretty_print_dict(result)}')
+
+
+def pretty_print_dict(dct):
+    retstr = ''
+    for key, value in sorted(dct.items(), key=lambda x: x[0]):
+        retstr += f"{key} : {value}\n"
+
+    return retstr
+
 if __name__ == '__main__':
     unittest.main()

tests/test_get_qualities_from_xml.py

@@ -0,0 +1,33 @@
+import unittest
+
+from pylot.core.io.phases import getQualitiesfromxml
+
+
+class TestQualityFromXML(unittest.TestCase):
+    def setUp(self):
+        self.path = '.'
+        self.ErrorsP = [0.02, 0.04, 0.08, 0.16]
+        self.ErrorsS = [0.04, 0.08, 0.16, 0.32]
+        self.test0_result = [[0.0136956521739, 0.0126, 0.0101612903226, 0.00734848484849, 0.0135069444444,
+                              0.00649659863946, 0.0129513888889, 0.0122747747748, 0.0119252873563, 0.0103947368421,
+                              0.0092380952381, 0.00916666666667, 0.0104444444444, 0.0125333333333, 0.00904761904762,
+                              0.00885714285714, 0.00911616161616, 0.0164166666667, 0.0128787878788, 0.0122756410256,
+                              0.013653253667966917], [0.0239333333333, 0.0223791578953, 0.0217974304255],
+                             [0.0504861111111, 0.0610833333333], [], [0.171029411765]], [
+                                [0.0195, 0.0203623188406, 0.0212121212121, 0.0345833333333, 0.0196180555556,
+                                 0.0202536231884, 0.0200347222222, 0.0189, 0.0210763888889, 0.018275862069,
+                                 0.0213888888889, 0.0319791666667, 0.0205303030303, 0.0156388888889, 0.0192,
+                                 0.0231349206349, 0.023625, 0.02875, 0.0195512820513, 0.0239393939394, 0.0234166666667,
+                                 0.0174702380952, 0.0204151307995], [0.040314343081226646], [0.148555555556], [], []]
+        self.test1_result = [77.77777777777777, 11.11111111111111, 7.407407407407407, 0, 3.7037037037037037],\
+                            [92.0, 4.0, 4.0, 0, 0]
+
+    def test_result_plotflag0(self):
+        self.assertEqual(getQualitiesfromxml(self.path, self.ErrorsP, self.ErrorsS, 0), self.test0_result)
+
+    def test_result_plotflag1(self):
+        self.assertEqual(getQualitiesfromxml(self.path, self.ErrorsP, self.ErrorsS, 1), self.test1_result)
+
+
+if __name__ == '__main__':
+    unittest.main()

tests/test_pick_correlation_correction/test_pick_correlation_correction.py

@@ -0,0 +1,76 @@
+import pytest
+from obspy import read, Trace, UTCDateTime
+
+from pylot.correlation.pick_correlation_correction import XCorrPickCorrection
+
+
+class TestXCorrPickCorrection():
+    def setup(self):
+        self.make_test_traces()
+        self.make_test_picks()
+        self.t_before = 2.
+        self.t_after = 2.
+        self.cc_maxlag = 0.5
+
+    def make_test_traces(self):
+        # take first trace of test Stream from obspy
+        tr1 = read()[0]
+        # filter trace
+        tr1.filter('bandpass', freqmin=1, freqmax=20)
+        # make a copy and shift the copy by 0.1 s
+        tr2 = tr1.copy()
+        tr2.stats.starttime += 0.1
+
+        self.trace1 = tr1
+        self.trace2 = tr2
+
+    def make_test_picks(self):
+        # create an artificial reference pick on reference trace (trace1) and another one on the 0.1 s shifted trace
+        self.tpick1 = UTCDateTime('2009-08-24T00:20:07.7')
+        # shift the second pick by 0.2 s, the correction should be around 0.1 s now
+        self.tpick2 = self.tpick1 + 0.2
+
+    def test_slice_trace_okay(self):
+
+        self.setup()
+        xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
+                                   t_before=self.t_before, t_after=self.t_after, cc_maxlag=self.cc_maxlag)
+
+        test_trace = self.trace1
+        pick_time = self.tpick2
+
+        sliced_trace = xcpc.slice_trace(test_trace, pick_time)
+        assert ((sliced_trace.stats.starttime == pick_time - self.t_before - self.cc_maxlag / 2)
+                and (sliced_trace.stats.endtime == pick_time + self.t_after + self.cc_maxlag / 2))
+
+    def test_slice_trace_fails(self):
+        self.setup()
+
+        test_trace = self.trace1
+        pick_time = self.tpick1
+
+        with pytest.raises(ValueError):
+            xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
+                                       t_before=self.t_before + 20, t_after=self.t_after, cc_maxlag=self.cc_maxlag)
+            xcpc.slice_trace(test_trace, pick_time)
+
+        with pytest.raises(ValueError):
+            xcpc = XCorrPickCorrection(UTCDateTime(), Trace(), UTCDateTime(), Trace(),
+                                       t_before=self.t_before, t_after=self.t_after + 50, cc_maxlag=self.cc_maxlag)
+            xcpc.slice_trace(test_trace, pick_time)
+
+    def test_cross_correlation(self):
+        self.setup()
+
+        # create XCorrPickCorrection object
+        xcpc = XCorrPickCorrection(self.tpick1, self.trace1, self.tpick2, self.trace2, t_before=self.t_before,
+                                   t_after=self.t_after, cc_maxlag=self.cc_maxlag)
+
+        # execute correlation
+        correction, cc_max, uncert, fwfm = xcpc.cross_correlation(False, '', '')
+
+        # define awaited test result
+        test_result = (-0.09983091718314982, 0.9578431835689154, 0.0015285160561610929, 0.03625786256084631)
+
+        # check results
+        assert pytest.approx(test_result, rel=1e-6) == (correction, cc_max, uncert, fwfm)