[initial] first implementation of residual plotting (WIP)

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

PyLoT.py

@@ -716,14 +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.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.init_spectro_tab()
         self.tabs.setCurrentIndex(0)
 
         self.eventLabel = QLabel()
@@ -1537,8 +1537,8 @@ class MainWindow(QMainWindow):
             return True
         return not bool(os.listdir(wf_path))
 
-    def filename_from_action(self, action):
-        if action.data() is None:
+    def filename_from_action(self, action=None):
+        if not action or action.data() is None:
             filt = "Supported file formats" \
                    " (*.mat *.qml *.xml *.kor *.evt)"
             caption = "Open an event file"
@@ -1976,7 +1976,6 @@ class MainWindow(QMainWindow):
         self.dataPlot.activateObspyDMToptions(self.obspy_dmt)
         if self.obspy_dmt:
             self.prepareObspyDMT_data(eventpath)
-            self.dataPlot.activateCompareOptions(True)
 
     def loadWaveformData(self):
         '''
@@ -2153,10 +2152,11 @@ class MainWindow(QMainWindow):
         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,
-                                                        pen=self.dataPlot.pen_linecolor)
+            self.dataPlot.plotWidget.getPlotItem().plot(np.array(times), np.array(data),
+                                                        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()
@@ -3096,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'))

README.md

@@ -54,33 +54,14 @@ In order to run PyLoT you need to install:
 
 #### Some handwork:
 
-PyLoT needs a properties folder on your system to work. It should be situated in your home directory
-(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
+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
 
@@ -89,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
@@ -97,7 +79,7 @@ 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
 
@@ -110,4 +92,4 @@ Others: A. Bruestle, T. Meier, W. Friederich
 
 [ObsPy]: http://github.com/obspy/obspy/wiki
 
-August 2024
+September 2024

pylot/core/io/phases.py

@@ -218,12 +218,14 @@ def picksdict_from_obs(fn):
     return picks
 
 
-def picksdict_from_picks(evt, parameter=None):
+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 = {
@@ -233,7 +235,10 @@ def picksdict_from_picks(evt, parameter=None):
     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:
@@ -254,7 +259,7 @@ def picksdict_from_picks(evt, parameter=None):
         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 = {}
@@ -301,7 +306,7 @@ def picksdict_from_picks(evt, parameter=None):
         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
 
 

pylot/core/util/array_map.py

@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
-
+import copy
 import traceback
 
 import cartopy.crs as ccrs
@@ -16,6 +16,8 @@ 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
 
@@ -41,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'
@@ -58,6 +61,8 @@ 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
 
@@ -108,8 +113,10 @@ class Array_map(QtWidgets.QWidget):
         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)
@@ -156,7 +163,9 @@ class Array_map(QtWidgets.QWidget):
         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):
@@ -217,7 +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)
@@ -374,8 +385,11 @@ class Array_map(QtWidgets.QWidget):
     def get_max_from_stations(self, key):
         return self._from_dict(max, key)
 
+    def get_selected_pick_type(self):
+        return self.comboBox_am.currentText().split(' ')[0]
+
     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}
@@ -447,6 +461,9 @@ class Array_map(QtWidgets.QWidget):
         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
         self.londim = self.lonmax - self.lonmin
@@ -459,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:
@@ -477,15 +497,20 @@ class Array_map(QtWidgets.QWidget):
 
     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 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))
-            latitudes.append(self.stations_dict[st_id]['latitude'])
-            longitudes.append(self.stations_dict[st_id]['longitude'])
+            uncertainties.append(self.uncertainties.get(nwst_id))
+            latitudes.append(self.stations_dict[nwst_id]['latitude'])
+            longitudes.append(self.stations_dict[nwst_id]['longitude'])
         return picks, uncertainties, latitudes, longitudes
 
     # plotting -----------------------------------------------------
@@ -582,7 +607,10 @@ class Array_map(QtWidgets.QWidget):
                                             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
 
@@ -737,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

@@ -59,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)

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.),
-                       'epp': (QtCore.Qt.DashLine, 1.),
-                       'lpp': (QtCore.Qt.DashLine, 1.),
-                       'spe': (QtCore.Qt.DashLine, 1.)}
-    linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2.),
-                       'epp': (QtCore.Qt.DashDotLine, 1.),
-                       'lpp': (QtCore.Qt.DashDotLine, 1.),
-                       'spe': (QtCore.Qt.DashDotLine, 1.)}
+    linestyles_manu = {'mpp': (QtCore.Qt.SolidLine, 2),
+                       'epp': (QtCore.Qt.DashLine, 1),
+                       'lpp': (QtCore.Qt.DashLine, 1),
+                       'spe': (QtCore.Qt.DashLine, 1)}
+    linestyles_auto = {'mpp': (QtCore.Qt.DotLine, 2),
+                       'epp': (QtCore.Qt.DashDotLine, 1),
+                       'lpp': (QtCore.Qt.DashDotLine, 1),
+                       'spe': (QtCore.Qt.DashDotLine, 1)}
     linestyles = {'manual': linestyles_manu,
                   'auto': linestyles_auto}
     return linestyles[picktype][key]
@@ -80,6 +81,7 @@ def which(program, parameter):
     return None
 
 
+@lru_cache(maxsize=128)
 def make_pen(picktype, phase, key, quality):
     """
     Make  PyQtGraph.QPen

tests/test_autopickstation/test_autopickstation.py

@@ -25,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:
@@ -271,7 +271,7 @@ 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)
         compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
         compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('GRA1', station)