BPMF.spectrum

class BPMF.spectrum.Spectrum(event=None)

Bases: object

Class for handling spectral data and calculations.$a

compute_correction_factor(rho_source, rho_receiver, vp_source, vp_receiver, vs_source, vs_receiver, radiation_S=np.float64(0.6324555320336759), radiation_P=np.float64(0.5163977794943222))

Compute the correction factor and attenuation factor for a seismic event.

Parameters:

• rho_source (float) – Density of the source medium, in kg/m3.

• rho_receiver (float) – Density of the receiver medium, in kg/m3.

• vp_source (float) – P-wave velocity of the source medium, in m/s.

• vp_receiver (float) – P-wave velocity of the receiver medium, in m/s.

• vs_source (float) – S-wave velocity of the source medium, in m/s.

• vs_receiver (float) – S-wave velocity of the receiver medium, in m/s.

• radiation_S (float, optional) – Radiation coefficient for S-wave. Default is sqrt(2/5).

• radiation_P (float, optional) – Radiation coefficient for P-wave. Default is sqrt(4/15).

Returns:

The correction factor and attenuation factor are stored in the object’s attributes geometrical_factor and attenuation_factor.

Return type:

None

Notes

This method requires the object to have an attached BPMF.dataset.Event instance and for the instance to have called the set_source_receiver_dist(network) method.

compute_multi_band_spectrum(traces, phase, buffer_seconds, **kwargs)

Compute spectrum from the maximum amplitude in multiple 1-octave frequency bands.

Parameters:

• traces (list) – List of seismic traces.

• phase (str) – Phase of the seismic event. Should be ‘noise’, ‘p’, or ‘s’.

• buffer_seconds (float) – Buffer duration in seconds to remove from the beginning and end of the trace.

• **kwargs – Additional keyword arguments for filtering.

Returns:

The computed spectrum is stored in the object’s attribute {phase}_spectrum.

Return type:

None

Notes

• The attribute frequency_bands is required for this method.

• The spectrum is computed by finding the maximum amplitude in each 1-octave frequency band.

• The resulting spectrum and frequency values are stored in the spectrum dictionary.

• The relative distance error percentage is calculated and stored for each trace.

• The attribute {phase.lower()}_spectrum is updated with the computed spectrum.

• If the attribute phases exists, the phase is appended to the list. Otherwise, a new list is created.

compute_network_average_spectrum(phase, snr_threshold, average_log=True, min_num_valid_channels_per_freq_bin=0, max_relative_distance_err_pct=25.0, reduce='mean', verbose=0)

Compute the network average spectrum for a given phase.

Parameters:

• phase (str) – Phase of the seismic event. Should be either ‘p’ or ‘s’.

• snr_threshold (float) – Signal-to-noise ratio threshold for valid channels.

• average_log (bool, optional) – Flag indicating whether to average the logarithm of the spectra. Default is True.

• min_num_valid_channels_per_freq_bin (int, optional) – Minimum number of valid channels required per frequency bin. Default is 0.

• max_relative_distance_err_pct (float, optional) – Maximum relative distance error percentage for a valid channel. Default is 25.0.

• verbose (int, optional) – Verbosity level. Set to 0 for no output. Default is 0.

Returns:

The average spectrum and related information are stored in the object’s attributes.

Return type:

None

Notes

This method requires the object to have already computed the spectrum for the specified phase and to have set the target frequencies using the set_target_frequencies method. If correct_propagation is set to True, the method also requires the object to have computed the correction factor using the compute_correction_factor method.

compute_signal_to_noise_ratio(phase)

Compute the signal-to-noise ratio (SNR) for the specified phase.

Parameters:

phase (str) – Phase of the seismic event. Should be ‘p’ or ‘s’.

Returns:

The computed SNR values are stored in the object’s attribute snr_{phase}_spectrum.

Return type:

None

Raises:

AssertionError –

• If phase is not ‘p’ or ‘s’. - If the {phase} spectrum has not been computed. - If the noise spectrum has not been computed.

Notes

• The SNR is calculated as the modulus of the signal spectrum divided by the modulus of the noise spectrum.

• The SNR values and corresponding frequencies are stored in the snr dictionary.

• The attribute snr_{phase}_spectrum is updated with the computed SNR values.

compute_spectrum(traces, phase, taper=None, **taper_kwargs)

Compute spectrum using the Fast Fourier Transform (FFT) on the input traces.

Parameters:

• traces (list) – List of seismic traces.

• phase (str) – Phase of the seismic event. Should be ‘noise’, ‘p’, or ‘s’.

• taper (callable or None, optional) – Tapering function to apply to the traces before computing the spectrum. Default is None.

• **taper_kwargs – Additional keyword arguments for the tapering function.

Returns:

The computed spectrum is stored in the object’s attribute {phase}_spectrum.

Return type:

None

Notes

• The spectrum is computed using the FFT on the input traces.

• The computed spectrum and frequency values are stored in the spectrum dictionary.

• The relative distance error percentage is calculated and stored for each trace.

• The attribute {phase.lower()}_spectrum is updated with the computed spectrum.

• If the attribute phases exists, the phase is appended to the list. Otherwise, a new list is created.

correct_attenuation()

Correct all spectra for attenuation.

correct_geometrical_spreading()

Correct all spectra for geometrical spreading.

differentiate(phase, average=True)

Apply differentiation to the spectrum of the specified phase.

Parameters:

• phase (str) – Phase of the seismic event. Should be ‘p’ or ‘s’.

• average (bool, optional) – Flag indicating whether to differentiate the average spectrum (True) or individual spectra (False). Defaults to True.

Returns:

The spectrum is modified in-place by multiplying it with the corresponding frequency values.

Return type:

None

Raises:

AssertionError –

• If average is True and the average {phase} spectrum has not been computed. - If average is False and the {phase} spectrum has not been computed.

fit_average_spectrum(phase, model='brune', log=True, min_fraction_valid_points_below_fc=0.1, min_fraction_valid_points=0.5, weighted=False, **kwargs)

Fit the average displacement spectrum with a specified model.

Parameters:

• phase (str) – Phase of the seismic event. Should be ‘p’ or ‘s’.

• model (str, optional) – Model to use for fitting the spectrum. Default is ‘brune’.

• log (bool, optional) – Flag indicating whether to fit the logarithm of the spectrum. Default is True.

• min_fraction_valid_points_below_fc (float, optional) – Minimum fraction of valid points required below the corner frequency. Default is 0.10.

• min_fraction_valid_points (float, optional) – Minimum fraction of valid points required overall. Default is 0.50.

• weighted (bool, optional) – Flag indicating whether to apply weighted fitting using sigmoid weights. Default is False.

• **kwargs – Additional keyword arguments to be passed to the curve fitting function.

Returns:

The fitting results are stored as attributes of the object.

Return type:

None

Raises:

AssertionError – If the average {phase} spectrum has not been computed.

Notes

• If the spectrum is below the SNR threshold everywhere, the fitting cannot be performed.

• If there are not enough valid points or valid points below the corner frequency, the fitting is considered unsuccessful.

integrate(phase, average=True)

Integrate the spectrum for a specific phase.

Parameters:

• phase (str) – Phase name. Should be ‘p’ or ‘s’.

• average (bool, optional) – Specifies whether to integrate the average spectrum (if True) or individual spectra (if False). Default is True.

Returns:

The spectrum is integrated in-place.

Return type:

None

Raises:

AssertionError – If the specified phase spectrum has not been computed.

plot_average_spectrum(phase, figname='spectrum', figtitle='', figsize=(10, 10), colors={'noise': 'dimgrey', 'p': 'C3', 's': 'black'}, linestyle={'noise': '--', 'p': '-', 's': '-'}, plot_fit=False, plot_std=False, plot_num_valid_channels=False)

Plot the average spectrum for a given phase.

Parameters:

• phase (str or list of str) – The phase(s) for which to plot the average spectrum.

• figname (str, optional) – The name of the figure. Default is “spectrum”.

• figtitle (str, optional) – The title of the figure. Default is an empty string.

• figsize (tuple, optional) – The size of the figure in inches. Default is (10, 10).

• colors (dict, optional) – A dictionary specifying the colors for different phases. Default is {“noise”: “dimgrey”, “s”: “black”, “p”: “C3”}.

• linestyle (dict, optional) – A dictionary specifying the line styles for different phases. Default is {“noise”: “–”, “s”: “-”, “p”: “-“}.

• plot_fit (bool, optional) – Whether to plot the fitted model spectrum. Default is False.

• plot_std (bool, optional) – Whether to plot the standard deviation range of the average spectrum. Default is False.

• plot_num_valid_channels (bool, optional) – Whether to plot the number of channels above the signal-to-noise ratio (SNR) threshold. Default is False.

Returns:

The generated figure.

Return type:

matplotlib.figure.Figure

plot_spectrum(phase, station=None, component=None, figname='spectrum', figsize=(10, 10), correct_propagation=False, plot_snr=False, colors={'noise': 'dimgrey', 'p': 'C3', 's': 'black'}, linestyle={'noise': '--', 'p': '-', 's': '-'})

Plot the spectrum for the specified phase(s) and trace(s).

Parameters:

• phase (str or list of str) – The phase(s) for which to plot the spectrum.

• station (str or None, optional) – The station code. If None, all stations will be plotted. Default is None.

• component (str or None, optional) – The component code. If None, all components will be plotted. Default is None.

• figname (str, optional) – The name of the figure. Default is “spectrum”.

• figsize (tuple, optional) – The size of the figure in inches. Default is (10, 10).

• correct_propagation (bool, optional) – Whether to correct the spectrum for propagation effects. Default is False.

• plot_snr (bool, optional) – Whether to plot the signal-to-noise ratio (SNR) spectrum if available. Default is False.

• colors (dict, optional) – A dictionary specifying the colors for different phases. Default is {“noise”: “dimgrey”, “s”: “black”, “p”: “C3”}.

• linestyle (dict, optional) – A dictionary specifying the line styles for different phases. Default is {“noise”: “–”, “s”: “-”, “p”: “-“}.

Returns:

The generated figure.

Return type:

matplotlib.figure.Figure

resample(new_frequencies, phase)

Resample the spectrum to new frequencies.

Parameters:

• new_frequencies (array-like) – New frequency values to resample the spectrum to.

• phase (str or list) – Phase(s) of the seismic event to resample. Can be a single phase or a list of phases.

Returns:

The spectrum is resampled and updated in-place.

Return type:

None

set_Q_model(Q, frequencies, Q_phase_prefactor={})

Set the attenuation Q model for P and S phases.

Parameters:

• Q (array-like) – Array of attenuation Q values.

• frequencies (array-like) – Array of corresponding frequencies. These are used to later keep track of which frequencies were used for the Q model.

Returns:

The Q model is stored in the Q0 attribute. The frequencies are stored in the Q0_frequencies attribute. These are later used to computed the Q-model and attenuation factor at arbitrary frequencies.

Return type:

None

set_frequency_bands(frequency_bands)

Set the frequency bands for spectrum analysis.

Parameters:

frequency_bands (dict) – Dictionary specifying the frequency bands. The keys are the names of the bands, and the values are tuples of the form (freqmin, freqmax). freqmin and freqmax represent the minimum and maximum frequencies of each band, respectively.

Returns:

The frequency bands are set and stored as an attribute in the object. The center frequencies of each band are also calculated and stored in the frequencies attribute.

Return type:

None

set_target_frequencies(freq_min, freq_max, num_points)

Set the target frequencies for spectrum analysis.

Parameters:

• freq_min (float) – Minimum frequency.

• freq_max (float) – Maximum frequency.

• num_points (int) – Number of target frequencies to generate.

Returns:

The target frequencies are set and stored as an attribute in the object.

Return type:

None

update_Q_model()

Interpolate the Q-model at the current self.frequencies.

update_attenuation_factor()

Compute attenuation factor at the current self.frequencies.

BPMF.spectrum.approximate_moment_magnitude(spectrum, snr_threshold=10.0, num_averaging_bands=1, low_snr_freq_min_hz=2.0, magnitude_log_moment_scaling=0.6666666666666666, phases=None, snr_based_weights=<function _snr_based_weights>)

Estimate the approximate moment magnitude (Mw*) from displacement spectra.

This function calculates Mw* by identifying valid frequency bands that satisfy a minimum SNR threshold. For high-quality signals, it prioritizes low-frequency bands to reflect the physical seismic moment. For low-quality signals, it applies a weighted average across available frequencies.

The magnitude is calculated using the relation: $$M_w = ext{scaling} imes (log_{10}(M_0) - 9.1)$$

Parameters:

• spectrum (BPMF.spectrum.Spectrum) – The spectrum object containing multi-band displacement and SNR data.

• snr_threshold (float, optional) – The minimum signal-to-noise ratio required to consider a frequency band “valid”. Default is 10.0.

• num_averaging_bands (int, optional) – The number of low-frequency valid bands to average when calculating the plateau. Default is 1.

• low_snr_freq_min_hz (float, optional) – The minimum frequency (in Hz) to consider when the SNR is below the threshold. Default is 2.0.

• magnitude_log_moment_scaling (float, optional) – The scaling factor for the moment-to-magnitude conversion. Default is 2/3 (standard Kanamori scale).

• phases (list of str or None, optional) – Seismic phases to process (e.g., [‘p’, ‘s’]). If None, retrieves phases from the spectrum object. Default is None.

• snr_based_weights (callable, optional) – A function used to calculate station weights based on peak SNR. Default is _snr_based_weights.

Returns:

approx_moment – A dictionary where keys are phase names and values are the calculated approximate moment magnitudes (Mw).

Return type:

dict

Notes

The function performs a station-by-station analysis. If no bands meet the snr_threshold, it falls back to a weighted mean of bands above low_snr_freq_min_hz to provide a best-effort estimate, though this may introduce significant error.

BPMF.spectrum.boatwright(freqs, omega0, fc, log=False)

Boatwright model.

BPMF.spectrum.brune(freqs, omega0, fc, log=False)

Brune model.

BPMF.spectrum.compute_moment_magnitude(event, windows, method='regular', phases=None, freq_min_hz=None, freq_max_hz=None, num_freqs=25, frequency_bands=None, window_buffer_sec=None, snr_threshold=10.0, min_num_valid_channels_per_freq_bin=3, max_relative_distance_err_pct=33.0, medium_properties={'Q_1Hz': None, 'attenuation_n': None, 'rho_receiver_kgm3': None, 'rho_source_kgm3': None, 'vp_receiver_ms': None, 'vp_source_ms': None, 'vs_receiver_ms': None, 'vs_source_ms': None}, approximate_moment_magnitude_args={'low_snr_freq_min_hz': 2.0, 'magnitude_log_moment_scaling': 0.6666666666666666, 'num_averaging_bands': 3}, q_phase_prefactor={'p': 2.25, 's': 1.0}, qc=True, full_output=False, spectral_model='brune', min_fraction_valid_points=0.5, min_fraction_valid_points_below_fc=0.2, num_channel_weighted_fit=True, max_rel_m0_err_pct=33.0, max_rel_fc_err_pct=33.0, stress_drop_mpa_min=0.001, stress_drop_mpa_max=10000.0, plot_above_mw=100.0, plot_above_random=1.0, plot_spectrum=False, figsize=(8, 8))

Apply workflow to estimate moment magnitude and approximate moment magnitude.

Parameters:

• event (BPMF.dataset.Event) – The event for which the moment magnitude must be estimated.

• windows (dict) – Seismic phase windows extracted using, for example, spectrum.extract_windows.

• method (str, optional) –

  Either of ‘regular’ or ‘multiband’. - ‘regular’: Calculates spectra with FFT.

  Must provide freq_min_hz and freq_max_hz.

  • ’multiband’: Calculates spectra following Al-Ismail et al., 2022. Must provide frequency_bands and window_buffer_sec.

• phases (list or None, optional) – List of seismic phase names, including ‘noise’. If None, uses windows to determine the phases. Defaults to None.

• freq_min_hz (float or None, optional) – Minimum frequency, in Hertz, at which spectra are resampled. Used only if method=’regular’. Defaults to None.

• freq_max_hz (float or None, optional) – Maximum frequency, in Hertz, at which spectra are resampled. Used only if method=’regular’. Defaults to None.

• num_freqs (int, optional) – Number of frequency bins used in the log-uniformly resampled spectra. Defaults to 25.

• frequency_bands (dict, optional) – Dictionary with frequency bands used when method=’multiband’. Example: frequency_bands = {‘band1’: [0.4, 0.8], ‘band2’: [0.8, 1.6]}

• window_buffer_sec (float, optional) – The left and right time series buffer, in seconds, to ignore when calculating the spectra with the ‘multiband’ method.

• snr_threshold (float, optional) – Only frequency bins with signal-to-noise ratio (snr) larger than snr_threshold are used in the network average. Defaults to 10.

• min_num_valid_channels_per_freq_bin (int, optional) – If less than min_num_valid_channels_per_freq_bin contributed to the network average at a given frequency, this frequency is masked and does not contribute to evaluating the spectral model. Defaults to 0.20.

• max_relative_distance_err_pct (float, optional) – If the relative distance error, 100 x dr / r, is more than max_relative_distance_err_pct, the station does not contribute to the network average. Defaults to 33%.

• medium_properties (dict, optional) – Medium properties used to relate the displacement spectrum to seismic moment release.

• approximate_moment_magnitude_args (dict, optional) – Arguments passed to BPMF.spectrum.approximate_moment_magnitude.

• qc (bool, optional) – If True, attempts to calculate the moment magnitude. Defaults to True.

• full_output (bool, optional) – If True, returns the propagation-corrected displacement spectra as well as snr spectra. Defaults to True.

• spectral_model (str, optional) – Either of ‘brune’ or ‘boatwright’. Defaults to ‘brune’.

• min_fraction_valid_points_below_fc (float, optional) – Fraction of valid frequency bins of the network-averaged spectrum below the inferred corner frequency required to trust the fit. Defaults to 0.20.

• num_channel_weighted_fit (bool, optional) – If True, fitting the spectral models gives more weight to the frequency bins at which the displacement spectra was observed across more stations.

• max_rel_m0_err_pct (float, optional) – Trust fit only if the relative error on seismic moment is less than max_rel_m0_err_pct. Defaults to 33%.

• max_rel_fc_err_pct (float, optional) – Trust fit only if the relative error on corner frequency is less than max_rel_fc_err_pct. Defaults to 33%.

• stress_drop_mpa_min (float, optional) – Trust fit only if the inferred stress drop is higher than stress_drop_mpa_min. Defaults to 1e-3.

• stress_drop_mpa_max (float, optional) – Trust fit only if the inferred stress drop is lower than stress_drop_mpa_max. Defaults to 1e4.

• plot_above_mw (float, optional) – Plot observed and modeled displacement spectrum if inferred moment magnitude is above plot_above_mw. Defaults to 100.

• plot_above_random (float, optional) – Plot observed and modeled displacement spectrum if random number is above plot_above_random. Defaults to 1.0.

• plot_spectrum (bool, optional) – If False, never plot figures. Defaults to False.

• figsize (tuple, optional) – Size of figure, in inches. Defaults to (8, 8).

Returns:

• spectrum (BPMF.spectrum.Spectrum) – The Spectrum instance built to compute and model the spectra.

• source_parameters (dict) – Dictionary with inferred moment magnitude and approximate moment magnitude.

• prop_corr_disp_spectra (pandas.DataFrame, optional) – Is returned only if full_output=True.

• snr_spectra (pandas.DataFrame, optional) – Is returned only if full_output=True.

• figures (list, optional) – Observed and modeled P- and/or S-wave displacement spectra. Is returned only if plot_spectrum=True and Mw > plot_above_mw or R > plot_above_random.

BPMF.spectrum.extract_windows(event, duration_sec, offset_ot_sec_noise, data_folder, attach_response=True, phase_on_comp_p={'1': 'P', '2': 'P', 'E': 'P', 'N': 'P', 'Z': 'P'}, phase_on_comp_s={'1': 'S', '2': 'S', 'E': 'S', 'N': 'S', 'Z': 'S'}, offset_phase={'P': 0.5, 'S': 0.5}, cleanup_stream=None)

Extract noise, P-wave, and S-wave windows and convert to displacement.

This function performs a three-stage extraction process from raw seismic data: 1. Extracts a noise window relative to the origin time (OT). 2. Extracts P-wave and S-wave windows based on phase arrival picks. 3. Applies detrending, tapering, and instrument response removal to

output displacement (m) seismograms.

Parameters:

• event (BPMF.dataset.Event) – The event object used to read waveforms and access arrival times.

• duration_sec (float) – The length of each extracted window in seconds.

• offset_ot_sec_noise (float) – The time offset from the origin time to start the noise window.

• data_folder (str) – Path to the directory containing the raw waveform files.

• attach_response (bool, optional) – If True, attaches instrument response metadata during reading. Default is True.

• phase_on_comp_p (dict, optional) – Mapping of component labels to the P-phase for windowing.

• phase_on_comp_s (dict, optional) – Mapping of component labels to the S-phase for windowing.

• offset_phase (dict, optional) – Offset in seconds relative to the phase arrival to start the signal window. Default is {“P”: 0.5, “S”: 0.5}.

• cleanup_stream (callable, optional) – A custom function or object (e.g., a filter) applied to the obspy.Stream after each read. Default is None.

Returns:

windows – A dictionary containing the processed obspy.Stream objects: - ‘noise’: The noise window traces. - ‘p’: The P-wave window traces. - ‘s’: The S-wave window traces.

Return type:

dict

Notes

Instrument response removal uses a pre-filter defined by the duration_sec and the Nyquist frequency to avoid numerical instability at very low or high frequencies. All outputs are returned as displacement (DISP).

BPMF.spectrum.fc_circular_crack(Mw, stress_drop_Pa=1000000.0, phase='p', vs_m_per_s=3500.0, vr_vs_ratio=0.9)

Compute the corner frequency assuming a circular crack model (Eshelby).

Parameters:

• Mw (float) – Moment magnitude of the earthquake.

• stress_drop_Pa (float, optional) – Stress drop in Pascals. Default is 1.0e6.

• phase (str, optional) – Seismic phase. Valid values are ‘p’ for P-wave and ‘s’ for S-wave. Default is ‘p’.

• vs_m_per_s (float, optional) – Shear wave velocity in meters per second. Default is 3500.0.

• vr_vs_ratio (float, optional) – Ratio of rupture velocity to shear wave velocity. Default is 0.9.

Returns:

corner_frequency – Corner frequency in Hertz.

Return type:

float

Raises:

AssertionError – If phase is not ‘p’ or ‘s’.

BPMF.spectrum.magnitude_to_moment(Mw)

Convert moment magnitude to seismic moment [N.m].

BPMF.spectrum.moment_to_magnitude(M0)

Convert seismic moment [N.m] to moment magnitude.

BPMF.spectrum.stress_drop_circular_crack(Mw, fc, phase='p', vs_m_per_s=3500.0, vr_vs_ratio=0.9)

Compute the stress drop assuming a circular crack model (Eshelby).

Parameters:

• Mw (float) – Moment magnitude of the earthquake.

• fc (float) – Corner frequency in Hertz.

• phase (str, optional) – Seismic phase. Valid values are ‘p’ for P-wave and ‘s’ for S-wave. Default is ‘p’.

• vs_m_per_s (float, optional) – Shear wave velocity in meters per second. Default is 3500.0.

• vr_vs_ratio (float, optional) – Ratio of rupture velocity to shear wave velocity. Default is 0.9.

Returns:

stress_drop – Stress drop in Pascals.

Return type:

float

Raises:

AssertionError – If phase is not ‘p’ or ‘s’.