Bugfix and improvements in source_residuals.py (#68)

* Removed extrapolation of individual residuals beyond their frequency range, which was leading to incorrect mean residuals at high frequencies.

* Added __main__ block in source_residuals.py.

* Added '--exclude' command-line argument and added 'exclude_subdirs' argument to read_residuals function.

* Added '--yrange' command-line argument and added 'ymin' and 'ymax' arguments to plot_residuals function.

* Save plots in subfolder if runid is specified.

* Added bugfix and new options in source_residuals.py to changelog.

* Added gridlines to plot in plot_residuals function.

* source_residuals: do not filter by runid if not all residuals have a runid attribute

* source_residuals: code linting and updated copyright

* Explain that the runid will be used to create a subdirectory within the output directory

* Default exclude subdirs to None

* source_residuals: accept wildcards in directory name

* Write a line differently for better readability

* Remove two points from the frequency array to avoid border effects

* Added 'runid' argument to plot_residuals function, to be displayed in the title.

* PR reference in CHANGELOG + style improvements

---------

Co-authored-by: Claudio Satriano <satriano@gmail.com>

Author: krisvanneste

CHANGELOG.md

@@ -27,8 +27,12 @@ Copyright (c) 2011-2025 Claudio Satriano <satriano@ipgp.fr>
 
 ### Post-Inversion
 
-- New option for `source_residuals`: `--runid` to select a specific run when
-  multiple runs exist for the same event
+- New options for `source_residuals`:
+  - `--runid` to select a specific run when multiple runs exist for the same
+     event
+  - `--exclude` to exclude 1 or more subfolders when parsing residuals files
+    (see [#68])
+  - `--yrange` to specify a fixed range for the Y axis in the plots (see [#68])
 
 ### Plotting
 
@@ -83,6 +87,9 @@ Copyright (c) 2011-2025 Claudio Satriano <satriano@ipgp.fr>
   `sensitivity` is a numerical value, any file format is accepted
 - Improved estimation of `fc_0` (initial corner frequency) when inverse
   frequency weighting is used (see [#67])
+- Fix in `source_residuals`: removed extrapolation of individual station
+  residuals beyond their valid frequency range, which was leading to incorrect
+  mean residuals at high frequencies (see [#68])
 
 ## v1.8 - 2024-04-07
 
@@ -818,3 +825,4 @@ Initial Python port.
 [#48]: https://github.com/SeismicSource/sourcespec/issues/48
 [#49]: https://github.com/SeismicSource/sourcespec/issues/49
 [#67]: https://github.com/SeismicSource/sourcespec/issues/67
+[#68]: https://github.com/SeismicSource/sourcespec/issues/68

sourcespec/source_residuals.py

@@ -7,13 +7,15 @@
     2013-2014 Claudio Satriano <satriano@ipgp.fr>,
               Agnes Chounet <chounet@ipgp.fr>
 
-    2015-2025 Claudio Satriano <satriano@ipgp.fr>
+    2015-2025 Claudio Satriano <satriano@ipgp.fr>.
+              Kris Vanneste <kris.vanneste@oma.be>
 :license:
     CeCILL Free Software License Agreement v2.1
     (http://www.cecill.info/licences.en.html)
 """
 import sys
 import os
+from glob import glob
 from collections import defaultdict
 from argparse import ArgumentParser
 import numpy as np
@@ -41,26 +43,43 @@ def parse_args():
     parser.add_argument(
         '-r', '--runid', dest='runid', action='store',
         default=None, metavar='RUNID',
-        help='Select a specific run when multiple runs exist for the same '
+        help='select a specific run when multiple runs exist for the same '
              'event. If omitted, residuals will be computed using all '
              'available runs. You can provide a specific RUNID or use '
              '"latest" or "earliest" to select the most recent or the '
-             'earliest run, based on alphabetical or numerical order.'
+             'earliest run, based on alphabetical or numerical order. '
+             'If provided, the runid will be also used to create a '
+             'subdirectory within the output directory, unless the '
+             'provided output directory name (see option "--outdir") is '
+             'already the runid.'
     )
     parser.add_argument(
         '-o', '--outdir', dest='outdir', action='store',
         default='sspec_residuals',
         help='output directory (default="sspec_residuals")'
     )
+    parser.add_argument(
+        '-e', '--exclude', dest='exclude_subdirs', action="append",
+        default=None,
+        help='subfolder to exclude (repeat for multiple subfolders)'
+    )
     parser.add_argument(
         '-p', '--plot', dest='plot', action='store_true',
         default=False, help='save residuals plots to file'
     )
+    parser.add_argument(
+        '-y', '--yrange', dest='yrange', nargs=2, type=float,
+        action='store', default=[-2.5, 2.5],
+        help='min/max range for Y axis (default: -2.5, 2.5)'
+    )
     parser.add_argument(
         'residual_files_dir',
-        help='directory containing source_spec residual files '
-             'in HDF5 format. Residual files can be in subdirectories '
-             '(e.g., a subdirectory for each event).')
+        help='path to a directory containing SourceSpec residual files in '
+             'HDF5 format. Files may be organized in subdirectories '
+             '(e.g., by event). You may use wildcards to match specific '
+             'subdirectories — remember to quote them to avoid shell '
+             'expansion. Example: "sspec_out/*/run1/"'
+    )
     return parser.parse_args()
 
 
@@ -83,8 +102,20 @@ def _filter_by_runid(residual_dict, runid='latest'):
     """
     if runid is None:
         return residual_dict
+    # Residuals generated by older versions of SourceSpec may not have a runid
+    # attribute. In this case, we cannot filter by runid.
+    if not all(
+        hasattr(spec.stats, 'runid')
+        for spec_st in residual_dict.values()
+        for spec in spec_st
+    ):
+        print(
+            'Warning: not all residuals have a "runid" attribute. '
+            'Filtering by runid will not be applied.'
+        )
+        return residual_dict
     if not isinstance(runid, str):
-        raise ValueError("runid must be a string.")
+        raise ValueError('runid must be a string.')
     filt_residual_dict = defaultdict(SpectrumStream)
     if runid not in ['latest', 'earliest']:
         for spec_st in residual_dict.values():
@@ -118,39 +149,81 @@ def _filter_by_runid(residual_dict, runid='latest'):
     return filt_residual_dict
 
 
-def read_residuals(resfiles_dir, runid=None):
+def _read_residuals_from_dir(
+        residual_dict, resfiles_dir, exclude_subdirs=None):
     """
-    Read residuals from HDF5 files in resfiles_dir.
+    Read residuals from HDF5 files from a specified directory and all its
+    subdirectories.
 
     Parameters
     ----------
+    residual_dict : dict
+        Dictionary to store residuals for each station.
     resfiles_dir : str
         Directory containing source_spec residual files in HDF5 format.
         Residual files can be in subdirectories (e.g., a subdirectory for
         each event).
-
-    Returns
-    -------
-    residual_dict : dict
-        Dictionary containing residuals for each station.
+    exclude_subdirs: list of str
+        List of subdirectories (potentially containing residual files)
+        that should not be parsed
+        (default: None, which means all subdirectories are parsed).
     """
+    if exclude_subdirs is None:
+        exclude_subdirs = []
     if not os.path.exists(resfiles_dir):
-        sys.exit(f'Error: directory "{resfiles_dir}" does not exist.')
+        print(f'Skipping directory: {resfiles_dir}: does not exist')
     resfiles = []
     for root, _dirs, files in os.walk(resfiles_dir):
+        _dirs[:] = [d for d in _dirs if d not in exclude_subdirs]
         resfiles.extend(
             os.path.join(root, file)
             for file in files
             if file.endswith('residuals.hdf5')
         )
     if not resfiles:
-        sys.exit(f'No residual file found in directory: {resfiles_dir}')
-    residual_dict = defaultdict(SpectrumStream)
+        print(f'No residual file found in directory: {resfiles_dir}')
     for resfile in sorted(resfiles):
         print(f'Found residual file: {resfile}')
         residual_st = read_spectra(resfile)
         for spec in residual_st:
             residual_dict[spec.id].append(spec)
+
+
+def read_residuals(resfiles_dir, runid=None, exclude_subdirs=None):
+    """
+    Read residuals from HDF5 files from a specified directory and all its
+    subdirectories. The directory can be specified using wildcards.
+
+    Parameters
+    ----------
+    resfiles_dir : str
+        Directory containing source_spec residual files in HDF5 format.
+        Residual files can be in subdirectories (e.g., a subdirectory for
+        each event). The subdirectory names can be specified using wildcards.
+    runid : str
+        Only select residuals with the specified runid. If 'latest' or
+        'earliest', select the most recent or the earliest runid,
+        respectively. If None, select all residuals.
+        (default: None)
+    exclude_subdirs: list of str
+        List of subdirectories (potentially containing residual files)
+        that should not be parsed
+        (default: None, which means all subdirectories are parsed).
+
+    Returns
+    -------
+    residual_dict : dict
+        Dictionary containing residuals for each station.
+    """
+    if exclude_subdirs is None:
+        exclude_subdirs = []
+    residual_dict = defaultdict(SpectrumStream)
+    for _resfiles_dir in glob(resfiles_dir):
+        _read_residuals_from_dir(
+            residual_dict,
+            _resfiles_dir,
+            exclude_subdirs=exclude_subdirs
+        )
     return _filter_by_runid(residual_dict, runid)
 
 
@@ -183,7 +256,9 @@ def compute_mean_residuals(residual_dict, min_spectra=20):
         delta_f_min = min(spec.stats.delta for spec in res)
         freq_min = min(freqs_min)
         freq_max = max(freqs_max)
-        freq_array = np.arange(freq_min, freq_max + delta_f_min, delta_f_min)
+        # stop two steps before the maximum frequency,
+        # to avoid border effects
+        freq_array = np.arange(freq_min, freq_max - delta_f_min, delta_f_min)
 
         spec_mean = None
         for spec in res:
@@ -193,10 +268,12 @@ def compute_mean_residuals(residual_dict, min_spectra=20):
             # spec_slice.data must exist, so we create it as zeros
             spec_interp.data = np.zeros_like(freq_array)
             # interpolate data_mag to the new frequency array
-            f = interp1d(spec.freq, spec.data_mag, fill_value='extrapolate')
+            f = interp1d(spec.freq, spec.data_mag, bounds_error=False)
             spec_interp.data_mag = f(freq_array)
-            # norm is 1 where data_mag is not zero, 0 otherwise
-            norm = (spec_interp.data_mag != 0).astype(float)
+            # norm is 1 where interpolated data_mag is not nan, 0 otherwise
+            norm = (~np.isnan(spec_interp.data_mag)).astype(float)
+            # Replace nan data_mag values with zeros for summation
+            spec_interp.data_mag[np.isnan(spec_interp.data_mag)] = 0
             if spec_mean is None:
                 spec_mean = spec_interp
                 norm_mean = norm
@@ -213,7 +290,8 @@ def compute_mean_residuals(residual_dict, min_spectra=20):
     return residual_mean
 
 
-def plot_residuals(residual_dict, residual_mean, outdir):
+def plot_residuals(residual_dict, residual_mean, outdir,
+                   ymin=None, ymax=None, runid=None):
     """
     Plot residuals.
 
@@ -225,6 +303,15 @@ def plot_residuals(residual_dict, residual_mean, outdir):
         SpectrumStream containing mean residuals for each station.
     outdir : str
         Output directory.
+    ymin : float
+        Lower limit of Y axis
+        (default: None)
+    ymax : float
+        Upper limit of Y axis
+        (default: None)
+    runid : str
+        Run ID, to be shown in title
+        (default: None)
     """
     for spec_mean in residual_mean:
         stat_id = spec_mean.id
@@ -234,9 +321,16 @@ def plot_residuals(residual_dict, residual_mean, outdir):
         for spec in res:
             ax.semilogx(spec.freq, spec.data_mag, 'b-', alpha=0.5)
         ax.semilogx(spec_mean.freq, spec_mean.data_mag, 'r-', linewidth=2)
+        ax.set_ylim([ymin, ymax])
+        ax.grid(
+            True, which='both', linestyle='solid', color='#DDDDDD', zorder=0)
         ax.set_xlabel('frequency (Hz)')
         ax.set_ylabel('residual amplitude (obs - synth) in magnitude units')
-        ax.set_title(f'residuals: {stat_id} – {len(res)} records')
+        if runid:
+            title = f'{stat_id} – runid: {runid} - {len(res)} records'
+        else:
+            title = f'{stat_id} – {len(res)} records'
+        ax.set_title(title)
         fig.savefig(figurefile, bbox_inches='tight')
         plt.close(fig)
         print(f'Residual plot saved to: {figurefile}')
@@ -246,27 +340,47 @@ def main():
     """Main function."""
     args = parse_args()
 
-    residual_dict = read_residuals(args.residual_files_dir, args.runid)
+    runid = args.runid
+    exclude_subdirs = args.exclude_subdirs
+    residual_dict = read_residuals(
+        args.residual_files_dir, runid,
+        exclude_subdirs=exclude_subdirs
+    )
+    if not residual_dict:
+        sys.exit(
+            'No residuals found. Please check the provided input directory.'
+        )
 
     outdir = args.outdir
+    # If runid is not None, create a subdirectory with the runid, unless
+    # the outdir name is already the runid
+    if runid and os.path.split(outdir)[-1] != runid:
+        outdir = os.path.join(outdir, runid)
     if not os.path.exists(outdir):
         os.makedirs(outdir)
 
     min_spectra = int(args.min_spectra)
     residual_mean = compute_mean_residuals(residual_dict, min_spectra)
 
     if args.plot:
-        plot_residuals(residual_dict, residual_mean, outdir)
+        ymin, ymax = args.yrange
+        plot_residuals(
+            residual_dict, residual_mean, outdir,
+            ymin=ymin, ymax=ymax, runid=runid
+        )
 
     # write the mean residuals (the stations corrections)
     res_mean_file = 'residual_mean.hdf5'
     res_mean_file = os.path.join(outdir, res_mean_file)
     if not residual_mean:
-        print(
+        sys.exit(
             'Mean residuals not computed: not enough spectra.\n'
             'Minimum number of spectra ("--min_spectra") currently set to: '
             f'{min_spectra}'
         )
-        sys.exit(0)
     residual_mean.write(res_mean_file, format='HDF5')
     print(f'Mean station residuals saved to: {res_mean_file}')
+
+
+if __name__ == '__main__':
+    main()