[initial] first implementation of residual plotting (WIP)

PyLoT.py

@@ -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"

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

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)