[update] increased robustness of correlation picker. If autoPyLoT fails on the stacked trace it tries to pick other stacked traces. Ignoring failed autoPyLoT picks can manually be set if they are not important (e.g. when only pick differences are important)

[minor] Commented out search for network name in metadata filename considered being unsafe

.gitignore

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

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()
@@ -1011,7 +1011,7 @@ class MainWindow(QMainWindow):
         for event in events:
             for filename in filenames:
                 if os.path.isfile(filename) and event.pylot_id in filename:
-                    self.load_data(filename, draw=False, event=event, ask_user=True, merge_strategy=sld.merge_strategy)
+                    self.load_data(filename, draw=False, event=event, ask_user=False, merge_strategy=sld.merge_strategy)
                     refresh = True
         if not refresh:
             return
@@ -1020,8 +1020,8 @@ class MainWindow(QMainWindow):
         self.fill_eventbox()
         self.setDirty(True)
 
-    def load_data(self, fname=None, loc=False, draw=True, event=None, ask_user=False, merge_strategy='Overwrite'):
-        if not ask_user:
+    def load_data(self, fname=None, loc=False, draw=True, event=None, ask_user=True, merge_strategy='Overwrite',):
+        if ask_user:
             if not self.okToContinue():
                 return
         if fname is None:
@@ -1034,7 +1034,7 @@ class MainWindow(QMainWindow):
         if not event:
             event = self.get_current_event()
 
-        if event.picks:
+        if event.picks and ask_user:
             qmb = QMessageBox(self, icon=QMessageBox.Question,
                               text='Do you want to overwrite the data?',)
             overwrite_button = qmb.addButton('Overwrite', QMessageBox.YesRole)
@@ -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'))
@@ -3590,7 +3590,7 @@ class MainWindow(QMainWindow):
     def calc_magnitude(self):
         self.init_metadata()
         if not self.metadata:
-            return None
+            return []
 
         wf_copy = self.get_data().getWFData().copy()
 
@@ -3599,6 +3599,10 @@ class MainWindow(QMainWindow):
         for station in np.unique(list(self.getPicks('manual').keys()) + list(self.getPicks('auto').keys())):
             wf_select += wf_copy.select(station=station)
 
+        if not wf_select:
+            logging.warning('Empty Stream in calc_magnitude. Return.')
+            return []
+
         corr_wf = restitute_data(wf_select, self.metadata)
         # calculate moment magnitude
         moment_mag = MomentMagnitude(corr_wf, self.get_data().get_evt_data(), self.inputs.get('vp'),

README.md

@@ -54,33 +54,17 @@ 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.
+
+## Example Dataset
+An example dataset with waveform data, metadata and automatic picks in the obspy-dmt dataset format for testing the teleseismic picking can be found at https://zenodo.org/doi/10.5281/zenodo.13759803
 
 ## Release notes
 
@@ -89,6 +73,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,17 +82,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): M. Rische, S. Wehling-Benatelli, L. Kueperkoch, M. Bischoff (PILOT)
+Developer(s): M. Paffrath, S. Wehling-Benatelli, L. Kueperkoch, D. Arnold, K. Cökerim, K. Olbert, M. Bischoff, C. Wollin, M. Rische, S. Zimmermann
 
-Developer(s): S. Wehling-Benatelli, M. Paffrath, L. Kueperkoch, K. Olbert, M. Bischoff, C. Wollin, M. Rische, D. Arnold, K. Cökerim, S. Zimmermann
+Original author(s): M. Rische, S. Wehling-Benatelli, L. Kueperkoch, M. Bischoff (PILOT)
 
 Others: A. Bruestle, T. Meier, W. Friederich
 
 
 [ObsPy]: http://github.com/obspy/obspy/wiki
 
-August 2024
+March 2025

docs/correlation.md

@@ -0,0 +1,77 @@
+# Pick-Correlation Correction
+
+## Introduction
+Currently, the pick-correlation correction algorithm is not accessible from they PyLoT GUI. The main file *pick_correlation_correction.py* is located in the directory *pylot\correlation*.
+The program only works for an obspy dmt database structure.
+
+The basic workflow of the algorithm is shown in the following diagram. The first step **(1)** is the normal (automatic) picking procedure in PyLoT. Everything from step **(2)** to **(5)** is part of the correlation correction algorithm.
+
+*Note: The first step is not required in case theoretical onsets are used instead of external picks when the parameter use_taupy_onsets is set to True. However, an existing event quakeML (.xml) file generated by PyLoT might be required for each event in case not external picks are used.*
+
+![images/workflow_stacking.png](images/workflow_stacking.png)
+
+A detailed description of the algorithm can be found in the corresponding publication:
+
+*Paffrath, M., Friederich, W., and the AlpArray and AlpArray-SWATH D Working Groups: Teleseismic P waves at the AlpArray seismic network: wave fronts, absolute travel times and travel-time residuals, Solid Earth, 12, 1635–1660, https://doi.org/10.5194/se-12-1635-2021, 2021.*
+
+## How to use
+To use the program you have to call the main program providing two mandatory arguments: a path to the obspy dmt database folder *dmt_database_path* and the path to the PyLoT infile *pylot.in* for picking of the beam trace:
+
+```python pick_correlation_correction.py dmt_database_path pylot.in```
+
+By default, the parameter file *parameters.yaml* is used. You can use the command line option *--params* to specify a different parameter file and other optional arguments such as *-pd* for plotting detailed information or *-n 4* to use 4 cores for parallel processing:
+
+```python pick_correlation_correction.py dmt_database_path pylot.in --params parameters_adriaarray.yaml -pd -n 4```
+
+## Cross-Correlation Parameters
+
+The program uses the parameters in the file *parameters.yaml* by default. You can use the command line option *--params* to specify a different parameter file. An example of the parameter file is provided in the *correlation\parameters.yaml* file.
+
+In the top level of the parameter file the logging level *logging* can be set, as well as a list of pick phases *pick_phases* (e.g. ['P', 'S']).
+
+For each pick phase the different parameters can be set in the first sub-level of the parameter file, e.g.:
+
+```yaml
+logging: info
+pick_phases: ['P', 'S']
+
+P:
+ min_corr_stacking: 0.8
+ min_corr_export: 0.6
+ [...]
+
+S:
+ min_corr_stacking: 0.7
+ [...]
+```
+
+The following parameters are available:
+
+
+| Parameter Name                 | Description                                                                                        | Parameter Type |
+|--------------------------------|----------------------------------------------------------------------------------------------------|----------------|
+| min_corr_stacking              | Minimum correlation coefficient for building beam trace                                            | float          |
+| min_corr_export                | Minimum correlation coefficient for pick export                                                    | float          |
+| min_stack                      | Minimum number of stations for building beam trace                                                 | int            |
+| t_before                       | Correlation window before reference pick                                                           | float          |
+| t_after                        | Correlation window after reference pick                                                            | float          |
+| cc_maxlag                      | Maximum shift for initial correlation                                                              | float          |
+| cc_maxlag2                     | Maximum shift for second (final) correlation (also for calculating pick uncertainty)               | float          |
+| initial_pick_outlier_threshold | Threshold for excluding large outliers of initial (AIC) picks                                      | float          |
+| export_threshold               | Automatically exclude all onsets which deviate more than this threshold from corrected taup onsets | float          |
+| min_picks_export               | Minimum number of correlated picks for export                                                      | int            |
+| min_picks_autopylot            | Minimum number of reference auto picks to continue with event                                      | int            |
+| check_RMS                      | Do RMS check to search for restitution errors (very experimental)                                  | bool           |
+| use_taupy_onsets               | Use taupy onsets as reference picks instead of external picks                                      | bool           |
+| station_list                   | Use the following stations as reference for stacking                                               | list[str]      |
+| use_stacked_trace              | Use existing stacked trace if found (spare re-computation)                                         | bool           |
+| data_dir                       | obspyDMT data subdirectory (e.g. 'raw', 'processed')                                               | str            |
+| pickfile_extension             | Use quakeML files (PyLoT output) with the following extension                                      | str            |
+| dt_stacking                    | Time difference for stacking window (in seconds)                                                   | list[float]    |
+| filter_options                 | Filter for first correlation (rough)                                                               | dict           |
+| filter_options_final           | Filter for second correlation (fine)                                                               | dict           |
+| filter_type                    | Filter type (e.g. bandpass)                                                                        | str            |
+| sampfreq                       | Sampling frequency (in Hz)                                                                         | float          |
+
+## Example Dataset
+An example dataset with waveform data, metadata and automatic picks in the obspy-dmt dataset format for testing can be found at https://zenodo.org/doi/10.5281/zenodo.13759803

docs/tuning.md

@@ -7,121 +7,122 @@ A description of the parameters used for determining automatic picks.
 Parameters applied to the traces before picking algorithm starts.
 
 | Name          | Description                                                                                                                                                                                                             |
-|---------------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| *P Start*, *P
-Stop* | Define time interval relative to trace start time for CF calculation on vertical trace. Value is relative to theoretical onset time if 'Use TauPy' option is enabled in main settings of 'Tune Autopicker' dialogue.    |
-| *S Start*, *S
-Stop* | Define time interval relative to trace start time for CF calculation on horizontal traces. Value is relative to theoretical onset time if 'Use TauPy' option is enabled in main settings of 'Tune Autopicker' dialogue. |
-| *Bandpass
-Z1*       | Filter settings for Butterworth bandpass applied to vertical trace for calculation of initial P pick.                                                                                                                   |
-| *Bandpass
-Z2*       | Filter settings for Butterworth bandpass applied to vertical trace for calculation of precise P pick.                                                                                                                   |
-| *Bandpass
-H1*       | Filter settings for Butterworth bandpass applied to horizontal traces for calculation of initial S pick.                                                                                                                |
-| *Bandpass
-H2*       | Filter settings for Butterworth bandpass applied to horizontal traces for calculation of precise S pick.                                                                                                                |
+|---------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
+| *P Start*, *P |                                                                                                                                                                                                                         |
+| Stop*         | Define time interval relative to trace start time for CF calculation on vertical trace. Value is relative to theoretical onset time if 'Use TauPy' option is enabled in main settings of 'Tune Autopicker' dialogue.    |
+| *S Start*, *S |                                                                                                                                                                                                                         |
+| Stop*         | Define time interval relative to trace start time for CF calculation on horizontal traces. Value is relative to theoretical onset time if 'Use TauPy' option is enabled in main settings of 'Tune Autopicker' dialogue. |
+| *Bandpass     |                                                                                                                                                                                                                         |
+| Z1*           | Filter settings for Butterworth bandpass applied to vertical trace for calculation of initial P pick.                                                                                                                   |
+| *Bandpass     |                                                                                                                                                                                                                         |
+| Z2*           | Filter settings for Butterworth bandpass applied to vertical trace for calculation of precise P pick.                                                                                                                   |
+| *Bandpass     |                                                                                                                                                                                                                         |
+| H1*           | Filter settings for Butterworth bandpass applied to horizontal traces for calculation of initial S pick.                                                                                                                |
+| *Bandpass     |                                                                                                                                                                                                                         |
+| H2*           | Filter settings for Butterworth bandpass applied to horizontal traces for calculation of precise S pick.                                                                                                                |
 
 ## Inital P pick
 
 Parameters used for determination of initial P pick.
 
 | Name                                                                                                              | Description                                                                                                            |
-|--------------|------------------------------------------------------------------------------------------------------------------------------|
-| *
-tLTA*       | Size of gliding LTA window in seconds used for calculation of HOS-CF.                                                        |
-| *pickwin
-P*  | Size of time window in seconds in which the minimum of the AIC-CF in front of the maximum of the HOS-CF is determined.       |
-| *
-AICtsmooth* | Average of samples in this time window will be used for smoothing of the AIC-CF.                                             |
-| *
-checkwinP*  | Time in front of the global maximum of the HOS-CF in which to search for a second local extrema.                             |
-| *minfactorP* | Used with *
-checkwinP*. If a second local maximum is found, it has to be at least as big as the first maximum * *minfactorP*. |
-| *
-tsignal*    | Time window in seconds after the initial P pick used for determining signal amplitude.                                       |
-| *
-tnoise*     | Time window in seconds in front of initial P pick used for determining noise amplitude.                                      |
-| *tsafetey*   | Time in seconds between *tsignal* and *
-tnoise*.                                                                              |
-| *
-tslope*     | Time window in seconds after initial P pick in which the slope of the onset is calculated.                                   |
+|-------------------------------------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------|
+| *                                                                                                                 |                                                                                                                        |
+| tLTA*                                                                                                             | Size of gliding LTA window in seconds used for calculation of HOS-CF.                                                  |
+| *pickwin                                                                                                          |                                                                                                                        |
+| P*                                                                                                                | Size of time window in seconds in which the minimum of the AIC-CF in front of the maximum of the HOS-CF is determined. |
+| *                                                                                                                 |                                                                                                                        |
+| AICtsmooth*                                                                                                       | Average of samples in this time window will be used for smoothing of the AIC-CF.                                       |
+| *                                                                                                                 |                                                                                                                        |
+| checkwinP*                                                                                                        | Time in front of the global maximum of the HOS-CF in which to search for a second local extrema.                       |
+| *minfactorP*                                                                                                      | Used with *                                                                                                            |
+| checkwinP*. If a second local maximum is found, it has to be at least as big as the first maximum * *minfactorP*. |                                                                                                                        |
+| *                                                                                                                 |                                                                                                                        |
+| tsignal*                                                                                                          | Time window in seconds after the initial P pick used for determining signal amplitude.                                 |
+| *                                                                                                                 |                                                                                                                        |
+| tnoise*                                                                                                           | Time window in seconds in front of initial P pick used for determining noise amplitude.                                |
+| *tsafetey*                                                                                                        | Time in seconds between *tsignal* and *                                                                                |
+| tnoise*.                                                                                                          |                                                                                                                        |
+| *                                                                                                                 |                                                                                                                        |
+| tslope*                                                                                                           | Time window in seconds after initial P pick in which the slope of the onset is calculated.                             |
 
 ## Inital S pick
 
 Parameters used for determination of initial S pick
 
 | Name                                                                                                              | Description                                                                                         |
-|---------------|------------------------------------------------------------------------------------------------------------------------------|
-| *
-tdet1h*      | Length of time window in seconds in which AR params of the waveform are determined.                                          |
-| *
-tpred1h*     | Length of time window in seconds in which the waveform is predicted using the AR model.                                      |
-| *
-AICtsmoothS* | Average of samples in this time window is used for smoothing the AIC-CF.                                                     |
-| *
-pickwinS*    | Time window in which the minimum in the AIC-CF in front of the maximum in the ARH-CF is determined.                          |
-| *
-checkwinS*   | Time in front of the global maximum of the ARH-CF in which to search for a second local extrema.                             |
-| *minfactorP*  | Used with *
-checkwinS*. If a second local maximum is found, it has to be at least as big as the first maximum * *minfactorS*. |
-| *
-tsignal*     | Time window in seconds after the initial P pick used for determining signal amplitude.                                       |
-| *
-tnoise*      | Time window in seconds in front of initial P pick used for determining noise amplitude.                                      |
-| *tsafetey*    | Time in seconds between *tsignal* and *
-tnoise*.                                                                              |
-| *
-tslope*      | Time window in seconds after initial P pick in which the slope of the onset is calculated.                                   |
+|-------------------------------------------------------------------------------------------------------------------|-----------------------------------------------------------------------------------------------------|
+| *                                                                                                                 |                                                                                                     |
+| tdet1h*                                                                                                           | Length of time window in seconds in which AR params of the waveform are determined.                 |
+| *                                                                                                                 |                                                                                                     |
+| tpred1h*                                                                                                          | Length of time window in seconds in which the waveform is predicted using the AR model.             |
+| *                                                                                                                 |                                                                                                     |
+| AICtsmoothS*                                                                                                      | Average of samples in this time window is used for smoothing the AIC-CF.                            |
+| *                                                                                                                 |                                                                                                     |
+| pickwinS*                                                                                                         | Time window in which the minimum in the AIC-CF in front of the maximum in the ARH-CF is determined. |
+| *                                                                                                                 |                                                                                                     |
+| checkwinS*                                                                                                        | Time in front of the global maximum of the ARH-CF in which to search for a second local extrema.    |
+| *minfactorP*                                                                                                      | Used with *                                                                                         |
+| checkwinS*. If a second local maximum is found, it has to be at least as big as the first maximum * *minfactorS*. |                                                                                                     |
+| *                                                                                                                 |                                                                                                     |
+| tsignal*                                                                                                          | Time window in seconds after the initial P pick used for determining signal amplitude.              |
+| *                                                                                                                 |                                                                                                     |
+| tnoise*                                                                                                           | Time window in seconds in front of initial P pick used for determining noise amplitude.             |
+| *tsafetey*                                                                                                        | Time in seconds between *tsignal* and *                                                             |
+| tnoise*.                                                                                                          |                                                                                                     |
+| *                                                                                                                 |                                                                                                     |
+| tslope*                                                                                                           | Time window in seconds after initial P pick in which the slope of the onset is calculated.          |
 
 ## Precise P pick
 
 Parameters used for determination of precise P pick.
 
 | Name                                                                                | Description                                                                                                                                                                                                      |
-|--------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| *Precalcwin* | Time window in seconds for recalculation of the HOS-CF. The new CF will be two times the size of *
-Precalcwin*, since it will be calculated from the initial pick to +/- *Precalcwin*.                                   |
-| *
-tsmoothP*   | Average of samples in this time window will be used for smoothing the second HOS-CF.                                                                                                                                    |
-| *
-ausP*       | Controls artificial uplift of samples during precise picking. A common local minimum of the smoothed and unsmoothed HOS-CF is found when the previous sample is larger or equal to the current sample times (1+*
-ausP*). |
+|-------------------------------------------------------------------------------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
+| *Precalcwin*                                                                        | Time window in seconds for recalculation of the HOS-CF. The new CF will be two times the size of *                                                                                                               |
+| Precalcwin*, since it will be calculated from the initial pick to +/- *Precalcwin*. |                                                                                                                                                                                                                  |
+| *                                                                                   |                                                                                                                                                                                                                  |
+| tsmoothP*                                                                           | Average of samples in this time window will be used for smoothing the second HOS-CF.                                                                                                                             |
+| *                                                                                   |                                                                                                                                                                                                                  |
+| ausP*                                                                               | Controls artificial uplift of samples during precise picking. A common local minimum of the smoothed and unsmoothed HOS-CF is found when the previous sample is larger or equal to the current sample times (1+* |
+| ausP*).                                                                             |                                                                                                                                                                                                                  |
 
 ## Precise S pick
 
 Parameters used for determination of precise S pick.
 
 | Name         | Description                                                                                                                                                                                                      |
-|--------------|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| *
-tdet2h*     | Time window for determination of AR coefficients.                                                                                                                                                                       |
-| *
-tpred2h*    | Time window in which the waveform is predicted using the determined AR parameters.                                                                                                                                      |
-| *Srecalcwin* | Time window for recalculation of ARH-CF. New CF will be calculated from initial pick +/- *
-Srecalcwin*.                                                                                                                  |
-| *
-tsmoothS*   | Average of samples in this time window will be used for smoothing the second ARH-CF.                                                                                                                                    |
-| *
-ausS*       | Controls artificial uplift of samples during precise picking. A common local minimum of the smoothed and unsmoothed ARH-CF is found when the previous sample is larger or equal to the current sample times (1+*
-ausS*). |
-| *
-pickwinS*   | Time window around initial pick in which to look for a precise pick.                                                                                                                                                    |
+|--------------|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
+| *            |                                                                                                                                                                                                                  |
+| tdet2h*      | Time window for determination of AR coefficients.                                                                                                                                                                |
+| *            |                                                                                                                                                                                                                  |
+| tpred2h*     | Time window in which the waveform is predicted using the determined AR parameters.                                                                                                                               |
+| *Srecalcwin* | Time window for recalculation of ARH-CF. New CF will be calculated from initial pick +/- *                                                                                                                       |
+| Srecalcwin*. |                                                                                                                                                                                                                  |
+| *            |                                                                                                                                                                                                                  |
+| tsmoothS*    | Average of samples in this time window will be used for smoothing the second ARH-CF.                                                                                                                             |
+| *            |                                                                                                                                                                                                                  |
+| ausS*        | Controls artificial uplift of samples during precise picking. A common local minimum of the smoothed and unsmoothed ARH-CF is found when the previous sample is larger or equal to the current sample times (1+* |
+| ausS*).      |                                                                                                                                                                                                                  |
+| *            |                                                                                                                                                                                                                  |
+| pickwinS*    | Time window around initial pick in which to look for a precise pick.                                                                                                                                             |
 
 ## Pick quality control
 
 Parameters used for checking quality and integrity of automatic picks.
 
 | Name                                       | Description                                                           |
-|--------------------------------------------|--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| *
-minAICPslope*                             | Initial P picks with a slope lower than this value will be discared.                                                                                                                                                                                                           |
-| *
-minAICPSNR*                               | Initial P picks with a SNR below this value will be discarded.                                                                                                                                                                                                                 |
-| *
-minAICSslope*                             | Initial S picks with a slope lower than this value will be discarded.                                                                                                                                                                                                          |
-| *
-minAICSSNR*                               | Initial S picks with a SNR below this value will be discarded.                                                                                                                                                                                                                 |
-| *minsiglength*, *noisefacor*. *minpercent* | Parameters for checking signal length. In the time window of size *
+|--------------------------------------------|-----------------------------------------------------------------------|
+| *                                          |                                                                       |
+| minAICPslope*                              | Initial P picks with a slope lower than this value will be discared.  |
+| *                                          |                                                                       |
+| minAICPSNR*                                | Initial P picks with a SNR below this value will be discarded.        |
+| *                                          |                                                                       |
+| minAICSslope*                              | Initial S picks with a slope lower than this value will be discarded. |
+| *                                          |                                                                       |
+| minAICSSNR*                                | Initial S picks with a SNR below this value will be discarded.        |
+| *minsiglength*, *noisefacor*. *minpercent* | Parameters for checking signal length. In the time window of size *   |
+
 minsiglength* after the initial P pick *
 minpercent* of samples have to be larger than the RMS value.                                                                                                        |
 | *
@@ -140,11 +141,11 @@ wdttolerance*                             | Maximum allowed deviation of S onset
 Parameters for discrete quality classes.
 
 | Name                                                                                                               | Description                                                                                                                                 |
-|------------------|-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
-| *
-timeerrorsP*    | Width of the time windows in seconds between earliest and latest possible pick which represent the quality classes 0, 1, 2, 3 for P onsets.                                                                           |
-| *
-timeerrorsS*    | Width of the time windows in seconds between earliest and latest possible pick which represent the quality classes 0, 1, 2, 3 for S onsets.                                                                           |
-| *nfacP*, *nfacS* | For determination of latest possible onset time. The time when the signal reaches an amplitude of *
-nfac* * mean value of the RMS amplitude in the time window *tnoise* corresponds to the latest possible onset time. |
+|--------------------------------------------------------------------------------------------------------------------|---------------------------------------------------------------------------------------------------------------------------------------------|
+| *                                                                                                                  |                                                                                                                                             |
+| timeerrorsP*                                                                                                       | Width of the time windows in seconds between earliest and latest possible pick which represent the quality classes 0, 1, 2, 3 for P onsets. |
+| *                                                                                                                  |                                                                                                                                             |
+| timeerrorsS*                                                                                                       | Width of the time windows in seconds between earliest and latest possible pick which represent the quality classes 0, 1, 2, 3 for S onsets. |
+| *nfacP*, *nfacS*                                                                                                   | For determination of latest possible onset time. The time when the signal reaches an amplitude of *                                         |
+| nfac* * mean value of the RMS amplitude in the time window *tnoise* corresponds to the latest possible onset time. |                                                                                                                                             |
 

inputs/pylot_global.in

@@ -41,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

@@ -41,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

@@ -41,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/core/io/data.py

@@ -36,8 +36,17 @@ class Data(object):
     loaded event. Container object holding, e.g. phase arrivals, etc.
     """
 
-    def __init__(self, parent=None, evtdata=None):
+    def __init__(self, parent=None, evtdata=None, picking_parameter=None):
         self._parent = parent
+
+        if not picking_parameter:
+            if hasattr(parent, '_inputs'):
+                picking_parameter = parent._inputs
+            else:
+                logging.warning('No picking parameters found! Using default input parameters!!!')
+                picking_parameter = PylotParameter()
+        self.picking_parameter = picking_parameter
+
         if self.getParent():
             self.comp = parent.getComponent()
         else:
@@ -403,23 +412,19 @@ class Data(object):
                                      not implemented: {1}'''.format(evtformat, e))
             if fnext == '.cnv':
                 try:
-                    velest.export(picks_copy, fnout + fnext, eventinfo=self.get_evt_data())
+                    velest.export(picks_copy, fnout + fnext, self.picking_parameter, eventinfo=self.get_evt_data())
                 except KeyError as e:
                     raise KeyError('''{0} export format
                                      not implemented: {1}'''.format(evtformat, e))
             if fnext == '_focmec.in':
                 try:
-                    parameter = PylotParameter()
-                    logging.warning('Using default input parameter')
-                    focmec.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
+                    focmec.export(picks_copy, fnout + fnext, self.picking_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:
-                    parameter = PylotParameter()
-                    logging.warning('Using default input parameter')
-                    hypodd.export(picks_copy, fnout + fnext, parameter, eventinfo=self.get_evt_data())
+                    hypodd.export(picks_copy, fnout + fnext, self.picking_parameter, eventinfo=self.get_evt_data())
                 except KeyError as e:
                     raise KeyError('''{0} export format
                                      not implemented: {1}'''.format(evtformat, e))

pylot/core/io/inputs.py

@@ -53,10 +53,16 @@ class PylotParameter(object):
         self.__parameter = {}
         self._verbosity = verbosity
         self._parFileCont = {}
+
         # io from parsed arguments alternatively
         for key, val in kwargs.items():
             self._parFileCont[key] = val
         self.from_file()
+
+        # if no filename or kwargs given, use default values
+        if not fnin and not kwargs:
+            self.reset_defaults()
+
         if fnout:
             self.export2File(fnout)
 

pylot/core/io/phases.py

@@ -278,7 +278,6 @@ def picksdict_from_picks(evt, parameter=None):
         weight = phase.get('weight')
         if not weight:
             if not parameter:
-                logging.warning('Using ')
                 logging.warning('Using default input parameter')
                 parameter = PylotParameter()
             pick.phase_hint = identifyPhase(pick.phase_hint)
@@ -513,7 +512,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         # 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('datapath'), parameter.get('eventID')))
-        arrivals = chooseArrivals(arrivals)  # MP MP what is chooseArrivals? It is not defined anywhere
+        arrivals = chooseArrivals(arrivals)
         for key in arrivals:
             # P onsets
             if 'P' in arrivals[key]:
@@ -666,7 +665,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
         fid = open("%s" % filename, 'w')
         # write header
         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
+        arrivals = chooseArrivals(arrivals)
         for key in arrivals:
             # P onsets
             if 'P' in  arrivals[key] and arrivals[key]['P']['mpp'] is not None:
@@ -757,11 +756,11 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
             cns, eventsource['longitude'], cew, eventsource['depth'], eventinfo.magnitudes[0]['mag'], ifx))
         n = 0
         # check whether arrivals are dictionaries (autoPyLoT) or pick object (PyLoT)
-        if isinstance(arrivals, dict) == False:
+        if isinstance(arrivals, dict) is False:
             # convert pick object (PyLoT) into dictionary
             evt = ope.Event(resource_id=eventinfo['resource_id'])
             evt.picks = arrivals
-            arrivals = picksdict_from_picks(evt)
+            arrivals = picksdict_from_picks(evt, parameter=parameter)
         # check for automatic and manual picks
         # prefer manual picks
         usedarrivals = chooseArrivals(arrivals)
@@ -822,7 +821,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
             # convert pick object (PyLoT) into dictionary
             evt = ope.Event(resource_id=eventinfo['resource_id'])
             evt.picks = arrivals
-            arrivals = picksdict_from_picks(evt)
+            arrivals = picksdict_from_picks(evt, parameter=parameter)
         # check for automatic and manual picks
         # prefer manual picks
         usedarrivals = chooseArrivals(arrivals)
@@ -873,7 +872,7 @@ def writephases(arrivals, fformat, filename, parameter=None, eventinfo=None):
             # convert pick object (PyLoT) into dictionary
             evt = ope.Event(resource_id=eventinfo['resource_id'])
             evt.picks = arrivals
-            arrivals = picksdict_from_picks(evt)
+            arrivals = picksdict_from_picks(evt, parameter=parameter)
         # check for automatic and manual picks
         # prefer manual picks
         usedarrivals = chooseArrivals(arrivals)

pylot/core/pick/autopick.py

@@ -262,6 +262,10 @@ class AutopickStation(object):
         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()
@@ -443,15 +447,15 @@ 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
-            estFirstS = 0
+            arrival_time_p = 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))
-            return estFirstP, estFirstS
+                  ' origin time using TauPy'.format(arrival_time_p, arrival_time_s))
+            return arrival_time_p, arrival_time_s
 
         def exit_taupy():
             """If taupy failed to calculate theoretical starttimes, picking continues.
@@ -477,10 +481,13 @@ 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["pstop"] += (self.origin[0].time + estFirstP) - self.ztrace.stats.starttime
+        self.pickparams["pstart"] += self.estFirstP
+        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
@@ -491,9 +498,10 @@ class AutopickStation(object):
             # 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.origin[0].time + estFirstS) - trace_s_start
-            self.pickparams["sstop"] += (self.origin[0].time + estFirstS) - trace_s_start
+            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
@@ -609,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')
@@ -713,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])
@@ -835,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(z_copy, cuttimes)
+        aiccf = AICcf(cf_stream, 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
@@ -855,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
@@ -894,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()
@@ -1146,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

@@ -60,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([])
 
@@ -212,17 +212,15 @@ class CharacteristicFunction(object):
             data = self.orig_data.copy()
             return data
 
-    def calcCF(self, data=None):
-        self.cf = data
+    def calcCF(self):
+        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:
         """
@@ -260,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), or order 2, i.e. STA/LTA.
-        :param data: data, time series (whether seismogram or CF)
-        :type data: tuple
         :return: HOS cf
         :rtype:
         """
@@ -281,47 +277,28 @@ class HOScf(CharacteristicFunction):
         elif self.getOrder() == 4:  # this is kurtosis
             y = np.power(xnp, 4)
             y1 = np.power(xnp, 2)
-        elif self.getOrder() == 2:  # this is variance, used for STA/LTA processing
-            y = np.power(xnp, 2)
-            y1 = np.power(xnp, 2)
 
         # Initialisation
         # t2: long term moving window
         ilta = int(round(self.getTime2() / self.getIncrement()))
-        ista = int(round((self.getTime2() / 10) / self.getIncrement()))  # TODO: still hard coded!!
         lta = y[0]
         lta1 = y1[0]
-        sta = y[0]
         # moving windows
         LTA = np.zeros(len(xnp))
-        STA = np.zeros(len(xnp))
         for j in range(0, len(xnp)):
             if j < 4:
                 LTA[j] = 0
-                STA[j] = 0
-            elif j <= ista and self.getOrder() == 2:
-                lta = (y[j] + lta * (j - 1)) / j
-                if self.getOrder() == 2:
-                    sta = (y[j] + sta * (j - 1)) / j
-            # elif j < 4:
             elif j <= ilta:
                 lta = (y[j] + lta * (j - 1)) / j
                 lta1 = (y1[j] + lta1 * (j - 1)) / j
-                if self.getOrder() == 2:
-                    sta = (y[j] - y[j - ista]) / ista + sta
             else:
                 lta = (y[j] - y[j - ilta]) / ilta + lta
                 lta1 = (y1[j] - y1[j - ilta]) / ilta + lta1
-                if self.getOrder() == 2:
-                    sta = (y[j] - y[j - ista]) / ista + sta
             # define LTA
             if self.getOrder() == 3:
                 LTA[j] = lta / np.power(lta1, 1.5)
             elif self.getOrder() == 4:
                 LTA[j] = lta / np.power(lta1, 2)
-            else:
-                LTA[j] = lta
-                STA[j] = sta
 
         # remove NaN's with first not-NaN-value,
         # so autopicker doesnt pick discontinuity at start of the trace
@@ -330,10 +307,7 @@ class HOScf(CharacteristicFunction):
             first = ind[0]
             LTA[:first] = LTA[first]
 
-        if self.getOrder() > 2:
         self.cf = LTA
-        else:  # order 2 means STA/LTA!
-            self.cf = STA / LTA
         self.xcf = x
 
 
@@ -343,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:
         """
@@ -479,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:
         """
@@ -635,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,7 +177,7 @@ 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
@@ -316,16 +320,7 @@ class AICPicker(AutoPicker):
                             plt.close(fig)
                     return
                 iislope = islope[0][0:imax + 1]
-            # MP MP change slope calculation
-            # 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]
+            dataslope = self.Data[0].data[iislope]
             # calculate slope as polynomal fit of order 1
             xslope = np.arange(0, len(dataslope), 1)
             try:
@@ -336,7 +331,7 @@ 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 Exception as e:
                 print("AICPicker: Problems with data fitting! {}".format(e))
 
@@ -356,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)
@@ -376,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'):
@@ -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)])
-        x1 = X[0][0:ilen]
-        x2 = X[1][0:ilen]
-        x3 = X[2][0:ilen]
+        earliest_starttime = min(tr.stats.starttime for tr in X)
+        cuttimes = common_range(X)
+        X = X.slice(cuttimes[0], cuttimes[1])
+        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)
 
@@ -874,6 +882,10 @@ def checksignallength(X, pick, minsiglength, pickparams, iplot=0, fig=None, line
         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()

pylot/core/util/array_map.py

@@ -134,8 +134,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 viridis as default
-        self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('viridis'))
+        # try to set to plasma as default
+        self.cmaps_box.setCurrentIndex(self.cmaps_box.findText('plasma'))
 
         self.top_row.addWidget(QtWidgets.QLabel('Select a phase: '))
         self.top_row.addWidget(self.comboBox_phase)

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)
@@ -218,9 +220,9 @@ class Metadata(object):
                     network_name = network.code
                     if not station_name in self.stations_dict.keys():
                         st_id = '{}.{}'.format(network_name, station_name)
-                        self.stations_dict[st_id] = {'latitude': station[0].latitude,
-                                                     'longitude': station[0].longitude,
-                                                     'elevation': station[0].elevation}
+                        self.stations_dict[st_id] = {'latitude': station.latitude,
+                                                     'longitude': station.longitude,
+                                                     'elevation': station.elevation}
 
         read_stat = {'xml': stat_info_from_inventory,
                      'dless': stat_info_from_parser}
@@ -266,9 +268,6 @@ class Metadata(object):
         if not fnames:
             # search for station name in filename
             fnames = glob.glob(os.path.join(path_to_inventory, '*' + station + '*'))
-        if not fnames:
-            # search for network name in filename
-            fnames = glob.glob(os.path.join(path_to_inventory, '*' + network + '*'))
         if not fnames:
             if self.verbosity:
                 print('Could not find filenames matching station name, network name or seed id')
@@ -280,7 +279,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)  # TODO: Commented out, failed with Parser, is this needed?
+                robj.get_coordinates(station_seed_id)
                 self.inventory_files[fname] = {'invtype': invtype,
                                                'data': robj}
                 if station_seed_id in self.seed_ids.keys():
@@ -288,6 +287,7 @@ class Metadata(object):
                 self.seed_ids[station_seed_id] = fname
                 return True
             except Exception as e:
+                logging.warning(e)
                 continue
         print('Could not find metadata for station_seed_id {} in path {}'.format(station_seed_id, path_to_inventory))
 
@@ -652,6 +652,8 @@ def restitute_data(data, metadata, unit='VEL', force=False, ncores=0):
     """
 
     # data = remove_underscores(data)
+    if not data:
+        return
 
     # loop over traces
     input_tuples = []
@@ -659,6 +661,11 @@ def restitute_data(data, metadata, unit='VEL', force=False, ncores=0):
         input_tuples.append((tr, metadata, unit, force))
         data.remove(tr)
 
+    if ncores == 0:
+        result = []
+        for input_tuple in input_tuples:
+            result.append(restitute_trace(input_tuple))
+    else:
         pool = gen_Pool(ncores)
         result = pool.imap_unordered(restitute_trace, input_tuples)
         pool.close()

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

pylot/core/util/utils.py

@@ -51,7 +51,6 @@ def readDefaultFilterInformation():
     :rtype: dict
     """
     pparam = PylotParameter()
-    pparam.reset_defaults()
     return readFilterInformation(pparam)
 
 
@@ -1076,7 +1075,7 @@ def check4rotated(data, metadata=None, verbosity=1):
             return wfs_in
 
         # check metadata quality
-        t_start = full_range(wfs_in)
+        t_start = full_range(wfs_in)[0]
         try:
             azimuths = []
             dips = []

pylot/correlation/parameters_adriaarray.yaml

@@ -9,7 +9,7 @@
 # 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
+# min_picks_autopylot: minimum number of reference auto 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
@@ -17,6 +17,11 @@
 # 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'
+# dt_stacking: time shift for stacking (e.g. [0, 250] for 0 and 250 seconds shift)
+# filter_options: filter for first correlation (rough)
+# filter_options_final: filter for second correlation (fine)
+# filter_type: e.g. 'bandpass'
+# sampfreq: sampling frequency of the data
 
 logging: info
 pick_phases: ['P', 'S']
@@ -54,6 +59,9 @@ P:
   filter_type: bandpass
   sampfreq: 20.0
 
+  # ignore if autopylot fails to pick master-trace (not recommended if absolute onset times matter)
+  ignore_autopylot_fail_on_master: True
+
 # S-phase
 S:
   min_corr_stacking: 0.7
@@ -88,3 +96,6 @@ S:
   filter_type: bandpass
   sampfreq: 20.0
 
+  # ignore if autopylot fails to pick master-trace (not recommended if absolute onset times matter)
+  ignore_autopylot_fail_on_master: True
+

pylot/correlation/pick_correlation_correction.py

@@ -436,6 +436,7 @@ def correlation_main(database_path_dmt: str, pylot_infile_path: str, params: dic
 
     # iterate over all events in "database_path_dmt"
     for eventindex, eventdir in enumerate(eventdirs):
+        eventindex += 1
         if not istart <= eventindex < istop:
             continue
 
@@ -443,12 +444,16 @@ def correlation_main(database_path_dmt: str, pylot_infile_path: str, params: dic
             continue
 
         logging.info('\n' + 100 * '#')
-        logging.info('Working on event {} ({}/{})'.format(eventdir, eventindex + 1, len(eventdirs)))
+        logging.info('Working on event {} ({}/{})'.format(eventdir, eventindex, len(eventdirs)))
         if event_blacklist and get_event_id(eventdir) in event_blacklist:
             logging.info('Event on blacklist. Continue')
 
+        try:
             correlate_event(eventdir, pylot_parameter, params=params, channel_config=channel_config,
                             update=update)
+        except Exception as e:
+            logging.error(f'Could not correlate event {eventindex}: {e}')
+            continue
 
     logging.info('Finished script after {} at {}'.format(datetime.now() - tstart, datetime.now()))
 
@@ -522,7 +527,7 @@ def cut_stream_to_same_length(wfdata: Stream) -> None:
         st.trim(tstart, tend, pad=True, fill_value=0.)
         # check for stream length
         if len(st) < 3:
-            logging.info('Not enough traces in stream, remove it from dataset:', str(st))
+            logging.info(f'Not enough traces in stream, remove it from dataset: {st}')
             remove(wfdata, st)
             continue
         if st[0].stats.starttime != st[1].stats.starttime or st[1].stats.starttime != st[2].stats.starttime:
@@ -635,10 +640,19 @@ def correlate_event(eventdir: str, pylot_parameter: PylotParameter, params: dict
             if len(picks) < params[phase_type]['min_picks_autopylot']:
                 logging.info('Not enough automatic picks for correlation. Continue!')
                 continue
-            # calculate corrected taupy picks and remove strong outliers
-            taupypicks_corr_initial, median_diff = get_corrected_taupy_picks(picks, taupypicks_orig)
-            picks = remove_outliers(picks, taupypicks_corr_initial,
-                                    params[phase_type]['initial_pick_outlier_threshold'])
+            ## calculate corrected taupy picks and remove strong outliers - part commented out for now. Maybe make this
+            ## another option but I think it is not worth it
+            # taupypicks_corr_initial, median_diff = get_corrected_taupy_picks(picks, taupypicks_orig)
+            # picks = remove_outliers(picks, taupypicks_corr_initial,
+            #                        params[phase_type]['initial_pick_outlier_threshold'])
+
+            # use mean/median corrected picks from taupy as the new reference picks. This largely increases the yield
+            # of final picks since every station has a decent reference onset now for a correlation
+            taupypicks, time_shift = get_corrected_taupy_picks(picks, taupypicks_orig, all_available=True)
+
+            picks = copy.deepcopy(taupypicks)
+            for pick in picks:
+                pick.method_id.id = 'auto'
 
         if phase_type == 'S':
             # check whether rotation to ZNE is possible (to get rid of horizontal channel 1,2,3)
@@ -671,19 +685,56 @@ def correlate_event(eventdir: str, pylot_parameter: PylotParameter, params: dict
         if not params[phase_type]['use_stacked_trace']:
             # first stack mastertrace by using stations with high correlation on that station in station list with the
             # highest correlation coefficient
-            stack_result = stack_mastertrace(wfdata_lowf, wfdata_highf, wfdata, picks, params=params[phase_type],
-                                             channels=channels_list, method=method, fig_dir=fig_dir)
-        else:
-            stack_result = load_stacked_trace(eventdir, params[phase_type]['min_corr_stacking'])
+            logging.info('Searching for master trace. ' + 20 * '-')
+            best_result = None
+            for stack_result in stack_mastertrace(wfdata_lowf, wfdata_highf, wfdata, picks, params=params[phase_type],
+                                                  channels=channels_list, method=method, fig_dir=fig_dir):
 
                 if not stack_result:
-            logging.info('No stack result. Continue.')
                     continue
 
+                # save best result in case it is still needed
+                if not best_result:
+                    best_result = stack_result
+
+                # extract stack result
+                correlations_dict, nwst_id, trace_master, nstack = stack_result
+
+                # now pick stacked trace with PyLoT for a more precise pick (use raw trace, gets filtered by autoPyLoT)
+                pick_stacked = repick_master_trace(wfdata_lowf, trace_master, pylot_parameter, event, event_id, metadata,
+                                                   phase_type, correlation_out_dir)
+
+                if not pick_stacked or not pick_stacked.time_errors.uncertainty:
+                    logging.info(f'Invalid autoPyLoT pick on master trace. Try next one.')
+                    continue
+                else:
+                    break
+
+            else:
+                logging.info('Did not find autoPyLoT pick for any stacked trace.')
+
+                if params[phase_type]['ignore_autopylot_fail_on_master'] is True:
+                    # Fallback in case autopylot failed
+                    logging.warning(f'Could not pick stacked trace. Using reference pick instead. If you need '
+                                    f'absolute onsets you need to account for a time-shift.')
+
+                    correlations_dict, nwst_id, trace_master, nstack = best_result
+                    pick_stacked = get_pick4station(picks, network_code=trace_master.stats.network,
+                                                    station_code=trace_master.stats.station, method='auto')
+                else:
+                    pick_stacked = None
+
+                if not pick_stacked:
+                    logging.info('Could not find reference pick also. Continue with next phase.')
+                    continue
+
+        else:
+            raise NotImplementedError('Loading stacked trace currently not implemented')
+            # stack_result = load_stacked_trace(eventdir, params[phase_type]['min_corr_stacking'])
+
+
         #############################################
         # NOW DO THE FINAL CORRELATION
-        # extract stack result
-        correlations_dict, nwst_id, trace_master, nstack = stack_result
 
         if params[phase_type]['plot']:
             # plot correlations of traces used to generate stacked trace
@@ -700,11 +751,6 @@ def correlate_event(eventdir: str, pylot_parameter: PylotParameter, params: dict
         # write unfiltered trace
         trace_master.write(os.path.join(correlation_out_dir, '{}_stacked.mseed'.format(trace_master.id)))
 
-        # now pick stacked trace with PyLoT for a more precise pick (use raw trace, gets filtered by autoPyLoT)
-        pick_stacked = repick_master_trace(wfdata_lowf, trace_master, pylot_parameter, event, event_id, metadata,
-                                           phase_type, correlation_out_dir)
-        if not pick_stacked:
-            continue
 
         # correlate stations with repicked and stacked master trace
         fig_dir_traces = make_figure_dirs(fig_dir, trace_master.id)
@@ -826,7 +872,7 @@ def get_picks_mean(picks: list) -> UTCDateTime:
 
 
 def get_corrected_taupy_picks(picks: list, taupypicks: list, all_available: bool = False) -> tuple:
-    """ get mean/median from picks taupy picks, correct latter for the difference """
+    """ get mean/median from picks relative to taupy picks, correct latter for the difference """
 
     def nwst_id_from_wfid(wfid):
         return '{}.{}'.format(wfid.network_code if wfid.network_code else '',
@@ -1165,26 +1211,26 @@ def stack_mastertrace(wfdata_lowf: Stream, wfdata_highf: Stream, wfdata_raw: Str
     A master trace will be created by stacking well correlating traces onto this station.
     """
 
-    def get_best_station4stack(sta_result):
-        """ return station with maximum mean_ccc"""
+    def get_best_stations4stack(sta_result, n_max=4):
+        """ return stations sorted after maximum mean_ccc"""
         ccc_means = {nwst_id: value['mean_ccc'] for nwst_id, value in sta_result.items() if
                      not np.isnan(value['mean_ccc'])}
         if len(ccc_means) < 1:
             logging.warning('No valid station found for stacking! Return.')
             return
-        best_station_id = max(ccc_means, key=ccc_means.get)
-        logging.info(
-            'Found highest mean correlation for station {} ({})'.format(best_station_id, max(ccc_means.values())))
-        return best_station_id
+        best_station_ids = sorted(ccc_means, key=ccc_means.get, reverse=True)[:n_max]
+        logging.info(f'Found mean correlations for stations: {ccc_means}')
+        return best_station_ids, ccc_means
 
     station_results = iterate_correlation(wfdata_lowf, wfdata_highf, channels, picks, method, params, fig_dir=fig_dir)
-    nwst_id_master = get_best_station4stack(station_results)
+    nwst_ids_master, ccc_means = get_best_stations4stack(station_results)
 
     # in case no stream with a valid pick is found
-    if not nwst_id_master:
-        logging.info('No mastertrace found! Will skip this event.')
+    if not nwst_ids_master:
+        logging.info('No mastertrace found! Continue with next.')
         return None
 
+    for nwst_id_master in nwst_ids_master:
         trace_master = station_results[nwst_id_master]['trace']
         stations4stack = station_results[nwst_id_master]['stations4stack']
         correlations_dict = station_results[nwst_id_master]['correlations_dict']
@@ -1196,7 +1242,7 @@ def stack_mastertrace(wfdata_lowf: Stream, wfdata_highf: Stream, wfdata_raw: Str
                                               do_rms_check=params['check_RMS'], plot=params['plot'], fig_dir=fig_dir,
                                               dt_pre=dt_pre, dt_post=dt_post)
 
-    return correlations_dict, nwst_id_master, trace_master, nstack
+        yield correlations_dict, nwst_id_master, trace_master, nstack
 
 
 def iterate_correlation(wfdata_lowf: Stream, wfdata_highf: Stream, channels: list, picks: list, method: str,
@@ -1984,4 +2030,5 @@ if __name__ == "__main__":
     # MAIN +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
     correlation_main(ARGS.dmt_path, ARGS.pylot_infile, params=CORR_PARAMS, istart=int(ARGS.istart),
                      istop=int(ARGS.istop),
-                     channel_config=CHANNELS, update=ARGS.update, event_blacklist=ARGS.blacklist)
+                     channel_config=CHANNELS, update=ARGS.update, event_blacklist=ARGS.blacklist,
+                     select_events=None) # can add a list of event ids if only specific events shall be picked

pylot/correlation/submit_pick_corr_correction.sh

@@ -32,9 +32,10 @@ export PYTHONPATH="$PYTHONPATH:/home/marcel/git/pylot/"
 #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'
+#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'
+# use -pd for detailed plots in eventdir/correlation_XX_XX/figures
+python pick_correlation_correction.py $1 $pylot_infile -n ${NSLOTS:=1} -istart $2 --params 'parameters_adriaarray.yaml' # -pd
 #--event_blacklist eventlist.txt

pylot/correlation/submit_to_grid_engine.py

@@ -3,14 +3,22 @@
 import subprocess
 
 fnames = [
-    ('/data/AlpArray_Data/dmt_database_synth_model_mk6_it3_no_rotation', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.0-5.4', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.4-5.7', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M5.7-6.0', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M6.0-6.3', 0),
+    ('/data/AdriaArray_Data/dmt_database_mantle_M6.3-10.0', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.0-5.4', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.4-5.7', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M5.7-6.0', 0),
+    # ('/data/AdriaArray_Data/dmt_database_ISC_mantle_M6.0-10.0', 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'
+script_location = '/home/marcel/VersionCtrl/git/pylot/pylot/correlation/submit_pick_corr_correction.sh'
 ####
 
 for fnin, istart in fnames:
@@ -18,6 +26,3 @@ for fnin, istart in fnames:
 
     print(input_cmds)
     print(subprocess.check_output(input_cmds.split()))
-
-
-

tests/test_autopickstation/autoPyLoT_global_taupy_false.in

@@ -41,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

@@ -41,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
@@ -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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        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_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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        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_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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        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'])
-        self.assertDictContainsSubset(expected=expected['S'], actual=result['S'])
+        compare_dicts(expected=expected['P'], result=result['P'], hint='P-')
+        compare_dicts(expected=expected['S'], result=result['S'], hint='S-')
         self.assertEqual('FIESA', station)
 
     def test_autopickstation_gra1_z_comp_missing(self):
@@ -272,7 +272,8 @@ class TestAutopickStation(unittest.TestCase):
         with HidePrints():
             result, station = autopickstation(wfstream=wfstream, pickparam=self.pickparam_taupy_disabled,
                                               metadata=(None, None))
-        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()