Merge branch 'refs/heads/develop' into correlation_picker

# Conflicts:
#	pylot/core/pick/charfuns.py
#	tests/test_autopicker/pylot_alparray_mantle_corr_stack_0.03-0.5.in
#	tests/test_autopicker/test_autopylot.py

PyLoT.py

@@ -1031,11 +1031,14 @@ class MainWindow(QMainWindow):
                 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 this data?',)
+                              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:
@@ -1045,8 +1048,6 @@ class MainWindow(QMainWindow):
             else:
                 return
 
-        if not event:
-            event = self.get_current_event()
         data = Data(self, event)
         try:
             data_new = Data(self, evtdata=str(fname))

autoPyLoT.py

@@ -184,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']:
@@ -234,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:
+                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':
-                data.setWFData(glob.glob(os.path.join(datapath, event_datapath, '*')))
+            if fnames in [None, 'None']:
+                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

pylot/core/pick/charfuns.py

@@ -231,7 +231,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)

pylot/core/util/array_map.py

@@ -14,6 +14,8 @@ import obspy
 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.util.utils import identifyPhaseID
 from scipy.interpolate import griddata
 
@@ -60,6 +62,7 @@ class Array_map(QtWidgets.QWidget):
         self.parameter = parameter if parameter else parent._inputs
 
         self.picks_rel = {}
+        self.picks_rel_mean_corrected = {}
         self.marked_stations = []
         self.highlighted_stations = []
 
@@ -106,6 +109,7 @@ class Array_map(QtWidgets.QWidget):
         self.map_reset_button = QtWidgets.QPushButton('Reset Map View')
         self.save_map_button = QtWidgets.QPushButton('Save Map')
         self.go2eq_button = QtWidgets.QPushButton('Go to Event Location')
+        self.subtract_mean_cb = QtWidgets.QCheckBox('Subtract mean')
 
         self.main_box = QtWidgets.QVBoxLayout()
         self.setLayout(self.main_box)
@@ -152,6 +156,7 @@ 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.subtract_mean_cb, 0)
         self.bot_row.addWidget(self.status_label, 5)
 
     def init_colormap(self):
@@ -212,6 +217,7 @@ 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.subtract_mean_cb.stateChanged.connect(self.toggle_subtract_mean)
 
         self.plotWidget.mpl_connect('motion_notify_event', self.mouse_moved)
         self.plotWidget.mpl_connect('scroll_event', self.mouse_scroll)
@@ -368,12 +374,6 @@ class Array_map(QtWidgets.QWidget):
     def get_max_from_stations(self, key):
         return self._from_dict(max, key)
 
-    def get_min_from_picks(self):
-        return min(self.picks_rel.values())
-
-    def get_max_from_picks(self):
-        return max(self.picks_rel.values())
-
     def current_picks_dict(self):
         picktype = self.comboBox_am.currentText().split(' ')[0]
         auto_manu = {'auto': self.autopicks_dict,
@@ -418,17 +418,17 @@ 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:
-                    picks_utc.append(pick)
+                if type(pick) is UTCDateTime:
+                    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:
-                        pick -= self._earliest_picktime
+                    if type(pick) is UTCDateTime:
+                        pick -= self._reference_picktime
                     picks_rel[st_id] = pick
             return picks_rel
 
@@ -437,6 +437,15 @@ class Array_map(QtWidgets.QWidget):
 
         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 init_lat_lon_dimensions(self):
         # init minimum and maximum lon and lat dimensions
@@ -467,11 +476,12 @@ 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 = []
         uncertainties = []
         latitudes = []
         longitudes = []
-        for st_id, pick in self.picks_rel.items():
+        for st_id, pick in picks_rel.items():
             picks.append(pick)
             uncertainties.append(self.uncertainties.get(st_id))
             latitudes.append(self.stations_dict[st_id]['latitude'])
@@ -538,13 +548,20 @@ class Array_map(QtWidgets.QWidget):
                 print(message, e)
                 print(traceback.format_exc())
 
-    def draw_contour_filled(self, nlevel=50):
-        levels = np.linspace(self.get_min_from_picks(), self.get_max_from_picks(), nlevel)
+    def draw_contour_filled(self, nlevel=51):
+        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.ax.contourf(self.longrid, self.latgrid, self.picksgrid_active, levels,
                                          linewidths=self.linewidth * 5, transform=ccrs.PlateCarree(),
                                          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())
 
@@ -575,7 +592,12 @@ class Array_map(QtWidgets.QWidget):
                           for uncertainty in uncertainties])
 
         cmap = self.get_colormap()
-        self.sc_picked = self.ax.scatter(lons, lats, s=sizes, edgecolors='white', cmap=cmap,
+
+        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):
@@ -652,7 +674,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)

pylot/core/util/dataprocessing.py

@@ -27,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):

tests/test_autopicker/pylot_alparray_mantle_corr_stack_0.03-0.5.in

@@ -4,9 +4,7 @@
 %Parameters are optimized for %extent data sets!
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 #main settings#
-                                                   #rootpath#       %project path                                                                                        
-                                                   #datapath#       %data path                                                                                           
-                                                   #database#       %name of data base                                                                                   
+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                                                                               

tests/test_autopicker/test_autopylot.py

@@ -1,6 +1,8 @@
 import os
 import pytest
 
+from obspy import read_events
+
 from autoPyLoT import autoPyLoT
 
 
@@ -8,7 +10,11 @@ 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
@@ -24,4 +30,38 @@ class TestAutopickerGlobal():
 
     def test_autopicker(self):
         assert self.init(), 'Initialization failed due to missing input files.'
-        #autoPyLoT(inputfile=self.params_infile, eventid='20171010_063224.a')
+        # 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'