feat: scan avec un taux d'échantillonnage différent du template (FDF 1998-2002)

CHANGELOG.md

@@ -4,6 +4,11 @@ Copyright (c) 2021-2026 Claudio Satriano <satriano@ipgp.fr>
 
 ## [unreleased]
 
+- Cross-correlation now handles continuous data whose sampling rate differs
+  from the template (for example the 20 Hz FDF recordings from 1998 to
+  2002): the finer trace is resampled to the common coarser rate instead of
+  aborting the scan, and the high-corner frequency is clamped just below the
+  Nyquist frequency so the bandpass stays valid at low sampling rates.
 - New interactive curses pager for all ``print_`` commands
   (``print_catalog``, ``print_pairs``, ``print_families``).
   Automatically activated when output is a terminal; use ``--no-pager``

requake/waveforms/waveforms.py

@@ -18,7 +18,7 @@
 from obspy.signal.cross_correlation import correlate, xcorr_max
 from obspy.clients.fdsn.header import FDSNException, FDSNNoDataException
 from scipy.stats import median_abs_deviation
-from ..config import config, rq_exit
+from ..config import config
 from ..wfcache import (
     clear_waveform_failure,
     read_cache_meta,
@@ -494,6 +494,22 @@ def process_waveforms(tr_or_st):
         config.args.freq_band[1] if getattr(config.args, 'freq_band', None)
         else config.cc_freq_max
     )
+    # Clamp the high corner just below the Nyquist frequency so the bandpass
+    # stays valid for low sampling-rate data (for example 20 Hz continuous
+    # recordings, where the default 10 Hz corner equals the Nyquist).
+    if isinstance(tr_or_st, Stream):
+        sampling_rate = min(tr.stats.sampling_rate for tr in tr_or_st)
+    else:
+        sampling_rate = tr_or_st.stats.sampling_rate
+    nyquist = 0.5 * sampling_rate
+    if freq_max >= nyquist:
+        clamped_freq_max = 0.99 * nyquist
+        logger.warning(
+            'High corner frequency %g Hz is at or above the Nyquist '
+            'frequency %g Hz; clamping to %g Hz.',
+            freq_max, nyquist, clamped_freq_max
+        )
+        freq_max = clamped_freq_max
     tr_or_st.filter(
         type='bandpass',
         freqmin=freq_min,
@@ -510,28 +526,32 @@ def process_waveforms(tr_or_st):
 
 def cc_waveform_pair(tr1, tr2, mode='events'):
     """Perform cross-correlation."""
-    dt1 = tr1.stats.delta
-    dt2 = tr2.stats.delta
-    if dt1 != dt2:
+    tr1 = process_waveforms(tr1)
+    tr2 = process_waveforms(tr2)
+    # Reconcile different sampling rates (for example a 100 Hz template
+    # against 20 Hz continuous data) by resampling the finer trace down to
+    # the common coarser rate. Both traces are already band-limited by
+    # process_waveforms, so resampling introduces no aliasing.
+    sampling_rate1 = tr1.stats.sampling_rate
+    sampling_rate2 = tr2.stats.sampling_rate
+    if sampling_rate1 != sampling_rate2:
+        target_rate = min(sampling_rate1, sampling_rate2)
         if mode == 'events':
-            evid1 = tr1.stats.evid
-            evid2 = tr2.stats.evid
             logger.warning(
-                f'{evid1} {evid2} - '
-                'The two traces have a different sampling interval. '
-                'Skipping pair.'
+                f'{tr1.stats.evid} {tr2.stats.evid} - '
+                'the two traces have a different sampling interval; '
+                f'resampling to {target_rate:g} Hz.'
             )
-        elif mode == 'scan':
-            logger.error(
-                'The two traces have a different sampling interval.')
-            rq_exit(1)
-    tr1 = process_waveforms(tr1)
-    tr2 = process_waveforms(tr2)
-    shift = int(config.cc_max_shift / dt1)
+        if sampling_rate1 > target_rate:
+            tr1.resample(target_rate)
+        if sampling_rate2 > target_rate:
+            tr2.resample(target_rate)
+    dt = tr1.stats.delta
+    shift = int(config.cc_max_shift / dt)
     cc = correlate(tr1, tr2, shift)
     abs_max = bool(config.cc_allow_negative)
     lag, cc_max = xcorr_max(cc, abs_max)
-    lag_sec = lag * dt1
+    lag_sec = lag * dt
     if mode != 'scan':
         return lag, lag_sec, cc_max
     # compute median absolute deviation for the non-zero portion of cc

tests/unit/test_cc_resampling.py

@@ -0,0 +1,80 @@
+# -*- coding: utf8 -*-
+# SPDX-License-Identifier: GPL-3.0-or-later
+"""
+Unit tests for cross-correlation with mismatched sampling rates.
+
+:copyright:
+    2021-2026 Claudio Satriano <satriano@ipgp.fr>,
+    Marius Yvard
+:license:
+    GNU General Public License v3.0 or later
+    (https://www.gnu.org/licenses/gpl-3.0-standalone.html)
+"""
+
+import unittest
+from argparse import Namespace
+from unittest.mock import patch
+import numpy as np
+from obspy import Trace, UTCDateTime
+from requake.config import config
+from requake.waveforms.waveforms import process_waveforms, cc_waveform_pair
+
+
+class TestCrossCorrelationResampling(unittest.TestCase):
+    """Cross-correlation should handle low and mismatched sampling rates."""
+
+    def setUp(self):
+        """Patch the cross-correlation config parameters."""
+        patcher = patch.multiple(
+            config,
+            create=True,
+            cc_freq_min=2.0,
+            cc_freq_max=10.0,
+            cc_max_shift=1.0,
+            cc_allow_negative=False,
+            args=Namespace(),
+        )
+        patcher.start()
+        self.addCleanup(patcher.stop)
+
+    def _make_trace(self, sampling_rate, duration=6.0):
+        """Build a band-limited synthetic trace at a given sampling rate."""
+        npts = int(duration * sampling_rate)
+        t = np.arange(npts) / sampling_rate
+        data = (
+            np.sin(2 * np.pi * 3.0 * t)
+            + 0.7 * np.sin(2 * np.pi * 5.0 * t)
+            + 0.5 * np.sin(2 * np.pi * 7.0 * t)
+        )
+        tr = Trace(data=data.astype('float64'))
+        tr.stats.sampling_rate = sampling_rate
+        tr.stats.starttime = UTCDateTime(0)
+        return tr
+
+    def test_process_waveforms_clamps_below_nyquist(self):
+        """At 20 Hz the 10 Hz corner is clamped just below the Nyquist."""
+        out = process_waveforms(self._make_trace(20.0))
+        self.assertLess(out.stats.freq_max, 10.0)
+        # A 100 Hz trace keeps the configured 10 Hz corner.
+        out_fast = process_waveforms(self._make_trace(100.0))
+        self.assertEqual(out_fast.stats.freq_max, 10.0)
+
+    def test_cc_mismatched_rates_resamples_instead_of_aborting(self):
+        """A 100 Hz template against 20 Hz data correlates after resampling."""
+        template = self._make_trace(100.0)
+        data = self._make_trace(20.0)
+        result = cc_waveform_pair(data, template, mode='scan')
+        self.assertEqual(len(result), 4)
+        _lag, _lag_sec, cc_max, _cc_mad = result
+        self.assertGreater(cc_max, 0.9)
+
+    def test_cc_equal_rates_still_correlates(self):
+        """Equal sampling rates keep correlating the same signal near 1."""
+        _lag, _lag_sec, cc_max = cc_waveform_pair(
+            self._make_trace(100.0), self._make_trace(100.0), mode='events'
+        )
+        self.assertGreater(cc_max, 0.99)
+
+
+if __name__ == '__main__':
+    unittest.main()