Fix ssp correction

CHANGELOG.md

@@ -24,6 +24,8 @@ Copyright (c) 2011-2025 Claudio Satriano <satriano@ipgp.fr>
 - New config parameter `clipping_min_amplitude_ratio` to set a threshold for
   trace amplitude below which the trace is not checked for clipping
 - Use spectral interpolation to compute and apply station residuals
+- Limit spectrum and residual to common frequency range when applying
+  correction before fitting
 
 ### Post-Inversion
 

sourcespec/ssp_build_spectra.py

@@ -490,7 +490,15 @@ def _build_weight_from_inv_frequency(spec, power=0.25):
 
 def _build_weight_from_ratio(spec, specnoise, smooth_width_decades):
     weight = spec.copy()
-    weight.data /= specnoise.data
+    try:
+        weight.data /= specnoise.data
+    except ValueError as e:
+        logger.error(
+            f'{spec.id}: Error building weight from noise: {str(e)}'
+        )
+        ssp_exit(1)
+    # replace NaN values with a small value
+    weight.data[np.isnan(weight.data)] = 1e-9
     # save signal-to-noise ratio before log10, smoothing, and normalization
     weight.snratio = weight.data.copy()
     # The inversion is done in magnitude units,
@@ -732,7 +740,7 @@ def _build_signal_and_noise_spectral_streams(
     for specnoise in specnoise_st:
         specnoise.data_mag = moment_to_mag(specnoise.data)
     # apply station correction if a residual file is specified in config
-    spec_st = station_correction(spec_st, config)
+    station_correction(spec_st, config)
     return spec_st, specnoise_st
 
 

sourcespec/ssp_correction.py

@@ -21,15 +21,49 @@
 logger = logging.getLogger(__name__.rsplit('.', maxsplit=1)[-1])
 
 
+def _interpolate_spectrum_to_new_freq_range(spec, new_freq):
+    """
+    Interpolate spectrum to a new frequency range.
+
+    Parameters
+    ----------
+    spec : Spectrum
+        Spectrum to be interpolated.
+    new_freq : array_like
+        New frequency range.
+
+    Returns
+    -------
+    Spectrum
+        Interpolated spectrum.
+    """
+    spec_interp = spec.copy()
+    spec_interp.freq = new_freq
+    # spec_interp.data must exist, so we create it with zeros
+    spec_interp.data = np.zeros_like(new_freq)
+    f = interp1d(
+        spec.freq, spec.data_mag, fill_value='extrapolate')
+    spec_interp.data_mag = f(new_freq)
+    return spec_interp
+
+
 def station_correction(spec_st, config):
     """
     Correct spectra using station-average residuals.
 
     Residuals are obtained from a previous run.
+
+    Parameters
+    ----------
+    spec_st : SpectrumStream or list of Spectrum
+        List of spectra to be corrected. Corrected spectra are appended
+        to the list (component code of uncorrected spectra is renamed to 'h').
+    config : Config
+        Configuration object containing the residuals file path.
     """
     res_filepath = config.residuals_filepath
     if res_filepath is None:
-        return spec_st
+        return
     try:
         residual = read_spectra(res_filepath)
     except Exception as msg:
@@ -42,33 +76,42 @@ def station_correction(spec_st, config):
             corr = residual.select(id=spec.id)[0]
         except IndexError:
             continue
-        freq = spec.freq
-        corr_interp = corr.copy()
-        corr_interp.freq = freq
-        # corr_interp.data must exist, so we create it with zeros
-        corr_interp.data = np.zeros_like(freq)
-        # interpolate the correction to the same frequencies as the spectrum
-        f = interp1d(
-            corr.freq, corr.data_mag, fill_value='extrapolate')
-        corr_interp.data_mag = f(freq)
+        # Define common frequency range for the correction and the spectrum
+        freq_min = max(spec.freq.min(), corr.freq.min())
+        freq_max = min(spec.freq.max(), corr.freq.max())
+        # Note that frequency range of corrected spectrum must not change,
+        # otherwise it will be out of sync with the noise spectrum
+        # and with the weight used in the inversion
+        # Instead, we use NaN values outside the common frequency range
+        # Interpolate residual to frequency range of spectrum,
+        # and set it to NaN outside the common frequency range
+        corr_interp = _interpolate_spectrum_to_new_freq_range(corr, spec.freq)
+        corr_interp.data_mag[corr_interp.freq < freq_min] = np.nan
+        corr_interp.data_mag[corr_interp.freq > freq_max] = np.nan
+        # Copy spectrum before correction
         spec_corr = spec.copy()
-        # uncorrected spectrum will have component name 'h'
+        # Uncorrected spectrum will have component name 'h',
+        # while corrected spectrum will have component name 'H'
         spec.stats.channel = f'{spec.stats.channel[:-1]}h'
+        # Apply correction
         try:
             spec_corr.data_mag -= corr_interp.data_mag
         except ValueError as msg:
             logger.error(f'Cannot correct spectrum {spec.id}: {msg}')
             continue
-        # interpolate the corrected data_mag to logspaced frequencies
-        f = interp1d(freq, spec_corr.data_mag, fill_value='extrapolate')
+        # Interpolate the corrected data_mag to logspaced frequencies
+        # We don't allow extrapolation, to make sure logspaced spectrum
+        # is also NaN outside the common frequency range
+        nan_idxs = np.isnan(spec_corr.data_mag)
+        f = interp1d(spec_corr.freq[~nan_idxs], spec_corr.data_mag[~nan_idxs],
+                     bounds_error=False)
         spec_corr.data_mag_logspaced = f(spec_corr.freq_logspaced)
-        # convert mag to moment
+        # Convert mag to moment
         spec_corr.data = mag_to_moment(spec_corr.data_mag)
         spec_corr.data_logspaced = mag_to_moment(spec_corr.data_mag_logspaced)
         spec_st.append(spec_corr)
-        fmin = freq.min()
-        fmax = freq.max()
+        fmin = spec_corr.freq[~nan_idxs].min()
+        fmax = spec_corr.freq[~nan_idxs].max()
         logger.info(
             f'{spec_corr.id}: corrected, frequency range is: '
             f'{fmin:.2f} {fmax:.2f} Hz')
-    return spec_st

sourcespec/ssp_inversion.py

@@ -255,6 +255,15 @@ def _spec_inversion(config, spec, spec_weight):
     # Initial values need to be printed here because Bounds can modify them
     logger.info(f'{statId}: initial values: {initial_values}')
     logger.info(f'{statId}: bounds: {bounds}')
+    # Remove NaN values from log-spaced spectrum
+    # (in case residual correction has been applied)
+    isnan = np.isnan(spec.data_mag_logspaced)
+    if np.sum(isnan) > 0:
+        spec.freq_logspaced = spec.freq_logspaced[~isnan]
+        spec.data_logspaced = spec.data_logspaced[~isnan]
+        spec.data_mag_logspaced = spec.data_mag_logspaced[~isnan]
+        weight = weight[~isnan]
+        yerr = yerr[~isnan]
     try:
         params_opt, params_err, misfit = _curve_fit(
             config, spec, weight, yerr, initial_values, bounds)

sourcespec/ssp_plot_stacked_spectra.py

@@ -264,8 +264,8 @@ def plot_stacked_spectra(config, spec_st, weight_st, sspec_output):
         # store min/max values for axes limits
         fmins.append(freqs.min())
         fmaxs.append(freqs.max())
-        specmins.append(spec.data.min())
-        specmaxs.append(spec.data.max())
+        specmins.append(np.nanmin(spec.data))
+        specmaxs.append(np.nanmax(spec.data))
     fmin = min(fmins)
     fmax = max(fmaxs)
     ax.set_xlim(fmin, fmax)

sourcespec/ssp_radiated_energy.py

@@ -45,8 +45,10 @@ def _spectral_integral(spec, t_star, fmin=None, fmax=None):
         data[freq < fmin] = 0.
     if fmax is not None:
         data[freq > fmax] = 0.
+    # Compute the integral and return it. Make sure to ignore
+    # NaN values in the data array.
     # Returned value has units of (m^2)^2 * Hz = m^4/s
-    return np.sum(data**2) * deltaf
+    return np.sum(data[~np.isnan(data)]**2) * deltaf
 
 
 def _radiated_energy_coefficient(rho, vel, free_surf_ampl, rp, average_rp):