API reference#

Grouped by what a run does in order: get the data, cut it, transform it, fit it, write it out.

If you are writing a package on top of SpecMod, start with specmod.api instead. It is a small, frozen subset of what follows, and the only part that carries a compatibility promise — the rest of this page documents internals that move between 0.x releases. See CONTRIBUTING.md for the exact scope of that promise.

The stable surface#

Stable public surface for downstream packages.

Anything not exported here is internal and may change without notice. Anything exported here follows the deprecation policy in CONTRIBUTING.md: one minor version of DeprecationWarning before a removal or a signature change, even while SpecMod is 0.x.

The point of the module is containment. SpecMod’s internals are still being refactored and its own documentation warns of breaking changes at every 0.x release; downstream packages import this and nothing else, so an internal rename costs a line here instead of a release there.

Five properties hold for everything below, and they are what make the surface usable from a service that owns its own IO and has to be able to replay a job:

Path-free. Every function takes in-memory data — arrays, or ObsPy objects. None of them opens a file. Convenience wrappers that take paths live elsewhere in the package.
Deterministic. The same inputs and the same explicit arguments produce the same outputs. Nothing here reads the working directory or the environment, and nothing draws random numbers. See the caveat on fit_spectrum().
Non-mutating. Inputs are left as they were found; results are new objects.
Quiet. Nothing prints. Diagnostics go through logging and warnings.
Typed errors. Failures are SpecModError subclasses — see specmod.exceptions for which of the three, and why the distinction is the useful part.

Examples

>>> import numpy as np
>>> from specmod import api
>>> rng = np.random.default_rng(0)
>>> signal = api.estimate_spectrum(rng.normal(size=2048), 0.01,
...                                estimator="multitaper")
>>> noise = api.estimate_spectrum(rng.normal(size=1024), 0.01,
...                               estimator="multitaper")
>>> pair = api.compare_spectra(signal, noise)
>>> pair.snr.shape == pair.binned_signal.freq.shape
True

class specmod.api.SpectrumFit(params, stderr, covariance, names, chisqr, redchi, n_points, success)[source]#

Bases: object

The result of fitting a source model to one spectrum.

Frozen, and holding plain numbers rather than the fitter’s own objects, so it can be serialised and compared without depending on lmfit’s API.

Parameters:

params (Mapping[str, float])
stderr (Mapping[str, float])
covariance (ndarray[tuple[Any, ...], dtype[float64]] | None)
names (tuple[str, ...])
chisqr (float)
redchi (float)
n_points (int)
success (bool)

params#

Fitted values, keyed by name — llpsp (the long-period spectral level Ω₀, as its base-10 logarithm), fc, ts (t*).

Type:: collections.abc.Mapping[str, float]

stderr#

One standard error per parameter, where the fitter could estimate one. Empty under some minimisers — see the note on fit_spectrum(). Absent means not measured, and is left absent rather than filled with a zero that would read as “certain”.

Type:: collections.abc.Mapping[str, float]

covariance#

The covariance matrix, with names giving its row and column order, or None when the minimiser produced none. The fc-t* correlation lives here, and reporting either parameter without it overstates both.

Type:: numpy.ndarray[tuple[Any, …], numpy.dtype[numpy.float64]] | None

chisqr, redchi: Misfit, and misfit per degree of freedom.

n_points#

How many spectral samples the fit actually used.

Type:: int

success#

Whether the minimiser reported convergence.

Type:: bool

correlation(a, b)[source]#

Correlation between two fitted parameters, or None.

None when there is no covariance matrix, or when either parameter has no variance to correlate — not zero, which would read as “independent” rather than “not measured”.

Parameters:

a (str)
b (str)

Return type:

float | None

specmod.api.available_estimators()[source]#

The estimators that can actually run in this environment, sorted.

SpecMod installs without its optional backends, so the registry is not the same question as what will work. Ask this before offering a choice to a user, rather than discovering the answer as a failed job.

Returns:: Names accepted by estimator= on estimate_spectrum().
Return type:: tuple of str

Examples

>>> "fft" in available_estimators()
True

specmod.api.compare_spectra(signal, noise, **settings)[source]#

Judge a signal spectrum against its noise window.

Returns the pair, including the per-bin signal-to-noise curve rather than only the band derived from it: pair.snr is an array aligned with pair.binned_signal.freq, and pair.band is one summary of it. A consumer that needs a different threshold, or that admits data bin by bin rather than over a contiguous interval, needs the curve — and a curve cannot be recovered from a stored interval.

Parameters:

signal (Spectrum) – Spectra from estimate_spectrum(). Neither is modified.
noise (Spectrum) – Spectra from estimate_spectrum(). Neither is modified.
**settings (Any) – threshold, f_min, f_max, n_bins, noise_model, bandwidth and the rest of specmod.core.collection.SpectrumPair.compare(). All have explicit defaults; none is read from configuration.

Returns:

With binned_signal, binned_noise, snr, band and resolution_floor.

Return type:

SpectrumPair

Raises:

InvalidInputError – The two spectra do not describe the same record geometry — a frequency axis above its own Nyquist, most often from pairing windows that came from different sampling rates.

specmod.api.config_to_toml(config, *, header=None)[source]#

Serialise a configuration to TOML, as specmod config freeze does.

Returns the text rather than writing it, so the caller decides where it goes — which for anything but a local filesystem is the only workable arrangement.

Parameters:

config (Config)
header (str | None)

Return type:

str

specmod.api.estimate_spectrum(data, dt, *, estimator, motion=Motion.VELOCITY, meta=None, **options)[source]#

Estimate the amplitude spectrum of one in-memory record.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str]) – The record, as a 1-D array of samples. Not modified.
dt (float) – Sample interval in seconds.
estimator (str) – Which backend, from available_estimators(). Required rather than defaulted: the configured default is a property of a study, and a service that resolves it silently cannot replay a job it did not record.
motion (Motion | str) – The ground-motion domain the record is in. Carried on the result, and what makes converting between domains a typed operation later.
meta (Mapping[str, Any] | None) – Extra metadata to attach to the spectrum. Copied, not held.
**options (Any) – Passed to the estimator’s constructor — n_tapers, time_bandwidth and so on. Backend-specific.

Returns:

Frequency axis, amplitude, and the metadata needed to interpret both.

Return type:

Spectrum

Raises:

InvalidInputError – The record is empty, not 1-D, contains non-finite values, or the estimator name is not known.
MissingBackendError – The estimator needs an optional extra that is not installed.

specmod.api.fit_spectrum(pair, *, id='', model=None, guess=None, fit_bins=False, method=None, weight_method=None, **fit_options)[source]#

Fit a source model to one spectrum, with its uncertainty.

This is the per-spectrum fit. SpecMod does not do a joint per-event inversion: specmod.fitting.FitSpectra loops over stations and fits each independently, sharing no parameters between them. A joint solver belongs to whoever needs one, on top of this.

Parameters:

pair (SpectrumPair) – From compare_spectra(). Its selected band is what gets fitted. Not modified.
id (str) – Label carried into the result’s metadata.
model (Any) – A model object, or None for the configured default.
guess (Mapping[str, float] | None) – Starting values for the fitted parameters. None derives them from the spectrum with specmod.fitting.initial_guess(), which is what specmod.fitting.FitSpectra does. Do not skip it: without a starting corner frequency the minimiser walks fc to zero and the model evaluates to NaN, so an unguessed fit does not merely fit worse, it raises.
fit_bins (bool) – Fit the log-binned spectrum rather than the full-resolution one.
method (str | None) – Minimiser name, passed to lmfit. None takes [fitting] method from configuration, which is what specmod.fitting.FitSpectra does — so a single-spectrum fit here matches the same station’s fit in an event run. Naming it explicitly is what makes the call reproducible somewhere else, and the default matters: on the 28 PNR windows lmfit’s own default returns a negative corner frequency on one station where the configured powell does not.
weight_method (str | None) – "log" weights residuals by 1/f; "none" does not. None takes [fitting] weight_method from configuration.
**fit_options (Any) – Anything else lmfit’s fit accepts.

Returns:

Point estimates and their errors and covariance. Frozen.

Return type:

SpectrumFit

Raises:

InvalidInputError – The pair has no usable band, or the spectrum is missing an attribute the model needs.

Notes

Uncertainty depends on the minimiser, and the configured default does not provide it. Only the least-squares family produces a covariance matrix. Measured on one synthetic station, all four agree on the corner frequency and only two report an error for it:

`method`	`fc`	`fc` error	`fc`-`t*` corr.
powell	7.925	–	–
nelder	7.925	–	–
leastsq	7.925	0.129	0.837

[fitting] method ships as powell, so a default fit returns point estimates with an empty stderr and no covariance. Pass method="leastsq" when the uncertainty is the point. That correlation is not incidental: 0.84 between fc and t* is why neither should be quoted alone.

One determinism caveat, and it is the only one on this surface. The initial guess and the default minimiser are read from configuration by internals, through specmod.config.load_config(), which resolves against the current working directory and the environment. Two runs in the same process with the same working directory agree exactly; two runs in different directories may not, if a specmod.toml differs between them.

Pass model and the minimiser options explicitly to close that gap, and record config_hash() alongside any result you intend to replay.

The exception hierarchy specmod.api raises.

Three kinds, because a caller does three different things with them:

InvalidInputError — the caller’s data or arguments are wrong. A program shows a form error; a person fixes the input.
MissingBackendError — the code is fine and the environment is not: an optional extra is not installed. Fixed by installing something, and avoidable up front with specmod.api.available_estimators().
InternalError — an invariant inside SpecMod is broken. Nothing the caller can do; it is a bug report.

Each also inherits the builtin exception the corresponding internal code raises today, so except ValueError keeps working and internals can migrate one at a time without a flag day.

Internals still raise the builtins. specmod.api translates at its own boundary, so the guarantee is specific: functions reached through `specmod.api` raise this hierarchy. Reaching around it gets the builtins.

exception specmod.exceptions.InternalError[source]#

Bases: SpecModError, RuntimeError

An invariant inside SpecMod does not hold.

Not caused by the caller and not fixable by them. If one of these reaches you, it is a bug in SpecMod.

exception specmod.exceptions.InvalidInputError[source]#

Bases: SpecModError, ValueError

The data or arguments given to SpecMod are not usable.

A record containing NaN, a frequency axis that does not belong to its record, an unknown estimator name, a band with no samples in it.

exception specmod.exceptions.MissingBackendError[source]#

Bases: SpecModError, ImportError

An optional backend is not installed.

Raised at call time rather than import time, so a default install stays importable. specmod.api.available_estimators() answers the same question without provoking the error.

exception specmod.exceptions.SpecModError[source]#

Bases: Exception

Base class for every error SpecMod raises deliberately.

The names excluded above are re-exports, documented at the path they are defined — Spectrum and SpectrumPair under Spectra, Config and load_config under Configuration, make_window and window_correction under Transforms. Documenting them twice gives every cross-reference to them two targets and makes all of them ambiguous, which is the same trap package-level automodule set earlier on this page. specmod.api.__all__ is the authoritative list, and tests/test_api_surface.py asserts it.

Packages are documented at the path you import from — specmod.picks.PickSet, not specmod.picks.base.PickSet. Documenting both the package and its submodules gave every re-exported name two targets and made every cross-reference to it ambiguous.

Getting data#

Where an event’s data lives on disk, and the events shipped with the repo.

The layout follows ObsPy’s mass_downloader, which writes waveforms/ and stations/ beneath a per-event directory. This adds the three things a spectral workflow needs alongside them:

tutorial/data/events/<origin>/
    event.xml                  # QuakeML: origin, magnitudes, uncertainties
    waveforms/                 # one miniSEED file per channel
    stations/inventory.xml     # StationXML for those channels
    picks/*.xml                # QuakeML picks (*.picks: Snuffler markers)
    spectra/*.h5               # computed spectra
    spectra/flatfiles/*.csv    # and their tabular export

EventDirectory resolves those paths and Event carries the hypocentre needed to set source-station geometry. Tests, tools/ and the documentation notebooks all read the layout from here.

Datasets that ship with the repository have their own loader — load_pnr_2019() — and need no network. Published ones are fetched by load(), cached by pooch and pinned by hash. specmod.acquire produces both.

specmod.datasets.EVENTS = PosixPath('tutorial/data/events')#: Event directories, relative to the repository root.

specmod.datasets.PNR_2019 = Event(origin='2019-08-26T07:30:47.000000Z', latitude=53.785021, longitude=-2.97078, depth_km=2.04, catalogue_magnitude=2.9, catalogue_magnitude_type='Mw')#

Preston New Road, 26 August 2019 — the induced event the tutorial and both golden references are built around, and the largest of the PNR-2 sequence. The origin time doubles as the directory name.

Mw 2.9 is from the PNR-2 catalogue published with Cuadrilla’s hydraulic-fracture monitoring (NGDC, 709cbc2f-af5c-4d09-a4ea-6deb5aa8c5d8), which gives surface_ML, surface_Mw and corrected_Mw all as 2.9.

The hypocentre is the catalogue’s, converted from its British National Grid easting/northing (336135.0, 432515.0; EPSG:27700) to WGS84. The catalogue gives depth as an elevation of -2040 m.

specmod.datasets.REGISTRY: dict[str, DatasetSpec] = {}#

they ship with the package and need no download.

Type:: Published datasets by name. Local datasets are not listed

class specmod.datasets.Dataset(event, paths, manifest=None)[source]#

Bases: object

One event’s data, wherever it came from.

The readers are methods rather than eager attributes because a dataset is often opened for its metadata alone, and reading a stream costs real time.

Parameters:

event (Event)
paths (EventDirectory)
manifest (dict[str, Any] | None)

manifest: dict[str, Any] | None = None#: The acquisition manifest, where the dataset was produced by specmod.acquire. None for data that ships with the package.

stream(pattern='*')[source]#

Read the waveforms. Raw counts — the response is not removed.

Parameters:: pattern (str)
Return type:: Any

inventory()[source]#

Read the station metadata, including responses.

Return type:: Any

catalog()[source]#

Read the event’s QuakeML.

Everything the catalogue said, rather than the four numbers Event keeps: origin uncertainties, every magnitude rather than the preferred one, agency and evaluation status.

Return type:: Any

class specmod.datasets.DatasetSpec(name, url, sha256, event, member='')[source]#

Bases: object

A published dataset: where to get it, and what it should hash to.

The hash is what makes a regression test mean anything. A config records intent and makes a dataset regenerable, but FDSN is not content-addressed, so re-running the config is not guaranteed to return the same bytes — see §5.2.2 of docs/REFACTOR_PLAN.md.

Versioning is by name. magna_2020_v1 and magna_2020_v2 are separate entries, so a result pinned to v1 keeps fetching v1 after v2 exists.

Parameters:

name (str)
url (str)
sha256 (str)
event (Event)
member (str)

sha256: str#: sha256:... of the archive, as pooch expects it.

member: str = ''#: Path within the unpacked archive holding the event directory.

class specmod.datasets.Event(origin, latitude, longitude, depth_km, catalogue_magnitude=None, catalogue_magnitude_type=None)[source]#

Bases: object

An earthquake: where its data sits, and the hypocentre it happened at.

origin, latitude, longitude and depth_km are the four values specmod.preprocess.set_stream_distance() takes.

Parameters:

origin (str)
latitude (float)
longitude (float)
depth_km (float)
catalogue_magnitude (float | None)
catalogue_magnitude_type (str | None)

catalogue_magnitude: float | None = None#: The published magnitude and the scale it is on, e.g. 2.9 and "Mw". Kept as a pair: ML and Mw diverge below about magnitude 3, so a bare number cannot be compared against a computed one.

directory(project_root)[source]#

Locate this event beneath project_root, the repository root.

Parameters:: project_root (Path | str)
Return type:: EventDirectory

class specmod.datasets.EventDirectory(root)[source]#

Bases: object

The paths beneath one event directory.

Reads nothing on construction, so it is safe to build at import time.

Parameters:: root (Path)

property inventory: Path#: The StationXML covering the channels in waveforms.

waveform_glob(pattern='*')[source]#

A glob over waveforms, as a string for obspy.read.

Parameters:: pattern (str)
Return type:: str

property quakeml: Path#

origin, magnitudes and their uncertainties.

Type:: The event’s QuakeML

picks_file()[source]#

The pick file for this event, QuakeML for preference.

QuakeML is the standard and carries what a marker file cannot — polarity, uncertainty, evaluation status, the full SEED id. Snuffler markers are still read where that is all there is.

Raises FileNotFoundError when neither is present.

Return type:: Path

is_present()[source]#

Whether the waveforms and station metadata are both present.

Does not check spectra/, which is generated rather than shipped.

Return type:: bool

specmod.datasets.data_dir()[source]#

Where downloaded datasets are cached.

SPECMOD_DATA_DIR overrides the platform cache directory, which matters on a cluster where $HOME is small or not writable from a compute node.

Return type:: Path

specmod.datasets.load(name, *, downloader=None)[source]#

Fetch a published dataset by name, from the cache after the first call.

Downloads are hash-checked by pooch: a corrupted or substituted archive fails here rather than quietly becoming a new expected answer.

downloader is passed through to pooch.retrieve(). It exists so the caching, hash check and unpacking can be exercised without a network — pooch has no file:// support — and so an operator behind an authenticating proxy can supply their own.

Parameters:

name (str)
downloader (Any)

Return type:

Dataset

specmod.datasets.load_pnr_2019()[source]#

The Preston New Road event committed to this repository.

Needs no download and no network: the waveforms and inventory are in the checkout. Raises when they are not, rather than reaching for a URL, since there is no published artefact for this one.

Return type:: Dataset

Fetch an event from an FDSN data centre into the layout tests and users read.

The request is declared in TOML and the response is written as an specmod.datasets.EventDirectory, beside a manifest recording what was asked for and what came back:

from specmod.acquire import fetch
fetch("datasets/pnr_2019.toml", out="build/pnr_2019")

or specmod fetch datasets/pnr_2019.toml -o build/pnr_2019.

Waveforms are stored raw. Counts and the response, never a deconvolved trace: baking remove_response into the artefact takes it out of test coverage and freezes one ObsPy version’s behaviour into the fixture.

A config makes the request reproducible, not the response. FDSN is not content-addressed — responses are corrected retroactively, archives are backfilled, catalogue solutions revised. That is what verify() and the manifest are for, and why published artefacts are pinned by hash rather than re-fetched. See §5.2.2 of docs/REFACTOR_PLAN.md.

Every network call goes through the client argument, which defaults to an ObsPy FDSN client and is injected in tests. Nothing here calls the network on import.

class specmod.acquire.AcquisitionConfig(name, data_centre='IRIS', event=<factory>, stations=<factory>, window=<factory>, source_toml='')[source]#

Bases: object

A complete, declarative description of one fetch.

Parameters:

name (str)
data_centre (str)
event (EventSpec)
stations (StationSpec)
window (WindowSpec)
source_toml (str)

data_centre: str = 'IRIS'#: FDSN data centre, by short name ("IRIS") or base URL. Recorded in the manifest because different centres serve different holdings for the same event.

source_toml: str = ''#: The TOML this was parsed from, kept verbatim for the manifest.

class specmod.acquire.EventSpec(eventid=None, catalogue=None, origin=None, latitude=None, longitude=None, depth_km=None, catalogue_magnitude=None, catalogue_magnitude_type=None)[source]#

Bases: object

Which earthquake, and where its parameters come from.

eventid resolves the hypocentre from the data centre’s catalogue, which is preferable to retyping it: a retyped origin is a second source of truth that can disagree with the catalogue silently. The explicit fields are for events the catalogue does not carry — induced sequences monitored privately, most often — and one of the two must be given.

Parameters:

eventid (str | None)
catalogue (str | None)
origin (str | None)
latitude (float | None)
longitude (float | None)
depth_km (float | None)
catalogue_magnitude (float | None)
catalogue_magnitude_type (str | None)

catalogue: str | None = None#: FDSN service to resolve eventid against, when it is not the one serving the waveforms. Event ids are issued per catalogue — a USGS ComCat id means nothing to IRIS — so the two are genuinely separable and the config has to be able to say so.

resolved()[source]#

The Event these fields describe.

Return type:: Event

class specmod.acquire.StationSpec(network='*', station='*', location='*', channel='*', max_radius_km=None, min_radius_km=None)[source]#

Bases: object

Which channels to ask for.

The patterns are FDSN wildcards, so the config alone does not say what you got — which is why the manifest records the channel list after expansion.

Parameters:

network (str)
station (str)
location (str)
channel (str)
max_radius_km (float | None)
min_radius_km (float | None)

max_radius_km: float | None = None#: Kilometres from the epicentre. None means no limit.

class specmod.acquire.WindowSpec(before_origin_s=10.0, after_origin_s=120.0)[source]#

Bases: object

How much record to take, relative to the origin time.

Parameters:

before_origin_s (float)
after_origin_s (float)

specmod.acquire.fetch(config, out, *, client=None, event_client=None)[source]#

Fetch one event and write it as an EventDirectory.

Returns the manifest, which is also written to manifest.json beside the data. client accepts anything with the ObsPy FDSN client’s get_events, get_stations and get_waveforms methods; tests pass a fake so that no test touches the network.

Parameters:

config (str | Path | AcquisitionConfig)
out (str | Path)
client (Any)
event_client (Any)

Return type:

dict[str, Any]

specmod.acquire.read_config(path)[source]#

Parse an acquisition config, keeping the text for the manifest.

Parameters:: path (str | Path)
Return type:: AcquisitionConfig

specmod.acquire.verify(out)[source]#

Re-hash what is on disk and report anything that no longer matches.

Integrity only: it says whether the files changed since they were written, not whether the data centre has revised its holdings. That needs a re-fetch and a diff against the manifest, which is the fuller --verify §5.2.2 describes and which needs the network.

Parameters:: out (str | Path)
Return type:: list[str]

Picks#

Phase arrivals, and matching them to the sensors a stream carries.

read() returns one PickSet per event, choosing a reader from PICK_READERS. select_event() narrows a multi-event source to one. resolve() matches those picks against a set of SensorID, returning a Resolution.

Three conditions have no single right answer and are therefore explicit: a source holding several events, a pick whose identity fits several sensors, and several picks for one sensor and phase. Each raises by default or takes a named policy. Design notes are in §4.9 of docs/REFACTOR_PLAN.md.

specmod.picks.BUILTIN_READERS: dict[str, PickReader] = {'obspy_events': ObsPyEventsReader(), 'snuffler': SnufflerReader()}#: Readers that ship with the package. Cannot be replaced by a plugin.

specmod.picks.ENTRY_POINT_GROUP = 'specmod.pick_readers'#: Entry point group a third party registers a reader under. A format that is an event file should register with ObsPy’s obspy.plugin.event instead, where it serves every ObsPy-based tool and reaches specmod through ObsPyEventsReader with no registration here at all.

specmod.picks.PICK_READERS: dict[str, PickReader] = {'obspy_events': ObsPyEventsReader(), 'snuffler': SnufflerReader()}#: Registered readers, by the name read(format=...) refers to them by. Plugins are added on first use; see load_plugins().

class specmod.picks.CSVPickReader(columns, reader_name='csv', delimiter=',', file_suffixes=('.csv',), skip_unknown_phases=True)[source]#

Bases: DelimitedPickReader

Comma-separated, with the quoting rules of csv.

Parameters:

columns (Mapping[str, str])
reader_name (str)
delimiter (str | None)
file_suffixes (tuple[str, ...])
skip_unknown_phases (bool)

class specmod.picks.DelimitedPickReader(columns, reader_name='delimited', delimiter=',', file_suffixes=('.csv', '.tsv', '.txt', '.dat'), skip_unknown_phases=True)[source]#

Bases: object

A delimited table of picks, one row per arrival.

columns maps this reader’s field names to the column headings in the file: {"station": "sta", "phase": "phase_type", "time": "arrival_time"}. station, phase and time are required; FIELDS lists what else may be mapped.

delimiter is a single character, or None to split on runs of whitespace — the latter takes a separate parse path, since quoting has no meaning in a whitespace-aligned table.

reader_name is what specmod.picks.read() refers to it by and must be unique once registered. file_suffixes is a hint for error messages and nothing else; detection is by header, so a reader claims a file only when every mapped column is present. Two readers configured for different schemas therefore do not collide, and one whose columns are a subset of another’s collides visibly — the ambiguity specmod.picks.detect_reader() reports.

Parameters:

columns (Mapping[str, str])
reader_name (str)
delimiter (str | None)
file_suffixes (tuple[str, ...])
skip_unknown_phases (bool)

skip_unknown_phases: bool = True#: Rows whose phase does not fold to P or S are skipped rather than raising.

class specmod.picks.ObsPyEventsReader[source]#

Bases: object

Anything obspy.read_events() parses.

class specmod.picks.Pick(sensor, phase, time, raw_phase=None, uncertainty=None, polarity=None, weight=None, automatic=None, reviewed=None, channel=None, author=None)[source]#

Bases: object

One phase arrival, with whatever provenance its format carried.

Only sensor, phase and time are always present. The rest is absent wherever the source format has no field for it, and is what a DuplicatePolicy uses to choose between competing picks.

Parameters:

sensor (SensorID)
phase (str)
time (UTCDateTime)
raw_phase (str | None)
uncertainty (float | None)
polarity (str | None)
weight (float | None)
automatic (bool | None)
reviewed (bool | None)
channel (str | None)
author (str | None)

class specmod.picks.PickReader(*args, **kwargs)[source]#

Bases: Protocol

One pick format.

suffixes is a hint for error messages and nothing else: a suffix does not identify a format, so specmod.picks.detect_reader() selects on can_read() alone.

name and suffixes are read-only properties so that a frozen dataclass satisfies the protocol; mutable attributes would not.

can_read(source)[source]#

Whether this reader recognises source.

Must not raise — on a missing, empty, truncated or binary file, and on every other format’s files, it returns False. Must be cheap: read a header, not the whole file.

Parameters:: source (str | PathLike[str])
Return type:: bool

read(source)[source]#

Every event in source, in file order.

Parameters:: source (str | PathLike[str])
Return type:: list[PickSet]

class specmod.picks.PickSet(picks=(), event_id=None, origin=None)[source]#

Bases: object

The picks of a single event.

Parameters:

picks (tuple[Pick, ...])
event_id (str | None)
origin (UTCDateTime | None)

sensors()[source]#

The distinct sensor identities picked, in sorted order.

Return type:: tuple[SensorID, …]

mapping(*, duplicates='prefer_reviewed')[source]#

The picks as {"NET.STA.LOC": {"P": UTCDateTime}}.

Keyed on each pick’s own identity, so an unstated field appears as * and matches no trace. Use resolve() to match against a stream.

Parameters:: duplicates (DuplicatePolicy)
Return type:: dict[str, dict[str, UTCDateTime]]

class specmod.picks.Resolution(attached=<factory>, unused=(), ambiguous=(), duplicated=())[source]#

Bases: object

The outcome of matching a PickSet against a set of sensors.

Parameters:

attached (dict[str, dict[str, Pick]])
unused (tuple[Pick, ...])
ambiguous (tuple[Pick, ...])
duplicated (tuple[tuple[str, str], ...])

attached: dict[str, dict[str, Pick]]#: Attached picks, keyed by the sensor’s complete id, then by phase.

unused: tuple[Pick, ...] = ()#: Picks that matched no sensor.

ambiguous: tuple[Pick, ...] = ()#: Picks skipped because they matched several sensors.

duplicated: tuple[tuple[str, str], ...] = ()#: (sensor, phase) where a policy chose between competing picks.

summary()[source]#

One line, suitable for printing or asserting on.

Return type:: str

class specmod.picks.SensorID(network, station, location)[source]#

Bases: object

A sensor identity, possibly partial.

None means the source did not state the field, and matches anything. "" means stated and empty, and matches only an empty field.

Readers apply that distinction asymmetrically: an empty network code is not a valid SEED network and becomes None, while an empty location code is the ordinary single-sensor case and is kept as "".

Parameters:

network (str | None)
station (str)
location (str | None)

classmethod parse(text)[source]#

Build from NET.STA.LOC, reading -- as an empty location.

Parameters:: text (str)
Return type:: SensorID

property is_complete: bool#: Whether every field is stated, so this names exactly one sensor.

matches(other)[source]#

Whether this identity is consistent with other.

Every field this one specifies must agree. Fields left None do not constrain the match, so a partial identity can match several sensors — which is a condition resolve() reports rather than resolves.

Parameters:: other (SensorID)
Return type:: bool

class specmod.picks.SnufflerReader[source]#

Bases: object

Marker files, as saved from Snuffler.

Not an ObsPy format: markers are an editor’s working notes rather than an event, and carry no origin.

class specmod.picks.TSVPickReader(columns, reader_name='tsv', delimiter='\t', file_suffixes=('.tsv', '.tab'), skip_unknown_phases=True)[source]#

Bases: DelimitedPickReader

Tab-separated.

Parameters:

columns (Mapping[str, str])
reader_name (str)
delimiter (str | None)
file_suffixes (tuple[str, ...])
skip_unknown_phases (bool)

class specmod.picks.WhitespacePickReader(columns, reader_name='whitespace', delimiter=None, file_suffixes=('.txt', '.dat', '.lst'), skip_unknown_phases=True)[source]#

Bases: DelimitedPickReader

Columns separated by runs of spaces or tabs.

The shape most hand-written and Fortran-era arrival tables come in. Cells cannot contain spaces, and quoting is not honoured.

This subsumes TSVPickReader — splitting on whitespace splits on tabs — so registering both against the same column names makes every tab-separated file ambiguous. Register one, or pass format=.

Parameters:

columns (Mapping[str, str])
reader_name (str)
delimiter (str | None)
file_suffixes (tuple[str, ...])
skip_unknown_phases (bool)

delimiter: str | None = None#

specmod.picks.detect_reader(source)[source]#

The one registered reader that recognises source.

Every reader is offered the file and exactly one must claim it. Both no claim and several are errors: a tie means two readers sniff too loosely, which is a bug in them rather than something to settle by priority.

Parameters:: source (str | PathLike[str])
Return type:: PickReader

specmod.picks.from_catalog(catalog)[source]#

Every event in an ObsPy catalogue, as one PickSet each.

Phase hints are folded to P or S on their first letter, with the original kept as raw_phase. Picks with no P/S hint, with an evaluation_status of rejected, or with no station code are dropped.

Parameters:: catalog (Catalog | ObsPyEvent)
Return type:: list[PickSet]

specmod.picks.get_reader(name)[source]#

Look a reader up by name.

Parameters:: name (str)
Return type:: PickReader

specmod.picks.load_plugins()[source]#

Discover readers advertised under ENTRY_POINT_GROUP.

Called on first use of the registry, so a plugin costs nothing to a caller who never reads a pick. A plugin that fails to import, or that claims a built-in name, warns naming its distribution and is skipped: a broken third-party reader must not make the built-in formats unreadable.

Return type:: None

specmod.picks.read(source, *, format=None)[source]#

Every event in source, in file order.

format names a reader in PICK_READERS and skips detection. Without it, detect_reader() selects one. An ObsPy Catalog is converted directly.

Parameters:

source (str | PathLike[str] | Catalog)
format (str | None)

Return type:

list[PickSet]

specmod.picks.register_reader(reader, *, replace=False)[source]#

Add a reader to PICK_READERS.

For a reader defined in a notebook or a script, where there is no installed distribution to hang an entry point on. A name already registered raises unless replace; a name in BUILTIN_READERS raises regardless.

Parameters:

reader (PickReader)
replace (bool)

Return type:

None

specmod.picks.resolve(picks, sensors, *, on_ambiguous='error', duplicates='prefer_reviewed')[source]#

Match a set of picks against the sensors actually present.

sensors are complete identities, as built from a stream. A pick reaches exactly one of them, none — in which case it is unused — or several, which is on_ambiguous:

error: Raise, naming the candidates and the fields the pick left unstated. The default.
skip: Leave the pick unattached, counted in Resolution.ambiguous.
broadcast: Attach to every match.

duplicates chooses between several picks for one sensor and phase: prefer_reviewed (then earliest), earliest, highest_weight, or error.

Parameters:

picks (PickSet)
sensors (Iterable[SensorID])
on_ambiguous (AmbiguousPolicy)
duplicates (DuplicatePolicy)

Return type:

Resolution

specmod.picks.select_event(sets, *, event_id=None, near=None, tolerance_s=60.0)[source]#

Choose one PickSet from a source that may hold several.

Selects by event_id, or by origin time within tolerance_s of near. With neither, the source must hold exactly one event. Raises unless exactly one event is selected, naming those available.

Parameters:

sets (Sequence[PickSet])
event_id (str | None)
near (UTCDateTime | None)
tolerance_s (float)

Return type:

PickSet

Preparing waveforms#

Geometry, picks and window cutting, on ObsPy streams.

This module and specmod.pipeline are the only two that know what a Trace is. ObsPy ships no type information and no stub package is published, so stubs/obspy in this repository declares the surface used here; see stubs/README.md. Stats is deliberately open, so tr.stats.delta is checked and tr.stats["p_time"] — one of the fields this module sets — is not.

specmod.preprocess.STREAM_DISTANCE_METHODS = ['mseed', 'sac', 'list', 'none']#: How station coordinates are supplied to set_stream_distance(). "none" is a deprecated alias for "list", kept because it is the spelling the function used to test for internally — see the note there.

specmod.preprocess.set_stream_distance(st, olat, olon, odep, ot, stlats=None, stlons=None, stelvs=None, inventory=None, dtype='sac')[source]#

Set the origin and source-receiver geometry on every trace in a stream.

dtype selects where the station coordinates come from:

"sac": the SAC header already on the trace.
"mseed": an ObsPy inventory, which is then required.
"list": the stlats, stlons and stelvs sequences, indexed positionally against the stream. "none" is a deprecated alias.

Anything else raises. It used to print invalid method choice and carry on, leaving traces with an origin but no distance and deferring the failure to whatever first asked for repi.

Parameters:

st (Stream)
olat (float)
olon (float)
odep (float)
ot (UTCDateTime)
stlats (Sequence[float] | None)
stlons (Sequence[float] | None)
stelvs (Sequence[float] | None)
inventory (Inventory | None)
dtype (str)

Return type:

None

specmod.preprocess.sensor_id(tr)[source]#

NET.STA.LOC — the sensor a pick belongs to.

Not the channel: an arrival is one sensor’s observation and is shared by its components. Not the bare station either: a borehole and a surface instrument differ only by location code, and they do not see the same arrival.

An empty location code is written --, matching how Snuffler marker files spell it.

Parameters:: tr (Trace)
Return type:: str

specmod.preprocess.read_picks(source, *, format=None)[source]#

Read picks as {"NET.STA.LOC": {"P": UTCDateTime, ...}}.

Accepts a Snuffler marker file or anything obspy.read_events() parses, so a caller holding a path from picks_file() does not have to know which format it got.

Keyed on each pick’s own identity, so a format supplying no network code keys on *.STA.* and matches no trace. set_picks() resolves against the sensors present instead.

Parameters:

source (str | PathLike[str])
format (str | None)

Return type:

dict[str, dict[str, Any]]

specmod.preprocess.set_picks(st, source, emergency_ratio=1.7, *, format=None, event_id=None, on_ambiguous='error', duplicates='prefer_reviewed', report=None)[source]#

Attach p_time and s_time to every trace with a pick.

source is a Snuffler marker file, a registered plugin’s format, or anything obspy.read_events() parses — QuakeML, SEISAN Nordic, HypoDD, NonLinLoc, a bulletin. See docs/pick-formats.md.

Picks are matched per sensor rather than per channel, so a pick made on one component reaches that sensor’s others, and a pick stating only part of an identity matches on the fields it does state. on_ambiguous and duplicates are passed to specmod.picks.resolve(); event_id to specmod.picks.select_event(), and is required for a multi-event source.

A trace whose P pick has no matching S gets one extrapolated at p + (p - otime) * emergency_ratio, unless that would place S before P, in which case s_time is left unset and a warning is issued.

Pass report=[] to receive the Resolution.

Parameters:

st (Stream)
source (str | PathLike[str])
emergency_ratio (float)
format (str | None)
event_id (str | None)
on_ambiguous (pk.AmbiguousPolicy)
duplicates (pk.DuplicatePolicy)
report (list[pk.Resolution] | None)

Return type:

None

specmod.preprocess.set_picks_from_pyrocko(st, pyrock_file, emergency_ratio=1.7)[source]#

Deprecated alias for set_picks().

Renamed because it no longer reads only Pyrocko: QuakeML is now the preferred format and the old name says the opposite of what the function does.

Parameters:

st (Stream)
pyrock_file (str | PathLike[str])
emergency_ratio (float)

Return type:

None

specmod.preprocess.basic_set_theoreticals(st, otime, p=5.9, s=2.9, dmetric='repi')[source]#

basic_set_theoreticals uses average propagation velocities [km/s] to set the arrival times for P and S waves. This assumes epicentral and/or and hypocentral distances have already been calculated and are set in the trace stats dictionary as tr.stats[‘repi’] or tr.stats[‘rhyp’] in units of kilometres.

Parameters:

st (Stream)
otime (UTCDateTime)
p (float)
s (float)
dmetric (str)

Return type:

None

specmod.preprocess.rstfl(fnames, wild='*', ext='sac')[source]#

rstfl reads create an obspy stream by reading each trace from an arbitrary list of paths.

Parameters:

fnames (Iterable[str])
wild (str)
ext (str)

Return type:

Stream

specmod.preprocess.link_window_to_trace(tr, start, end)[source]#

Record a window on a trace, as asked for and as delivered.

Two pairs, because they are not the same thing. trim gives back whatever the record actually holds, so a window that runs off either end comes back short — which is the normal case for noise windows, not the exception. Recording only the request is how a truncated noise trace ends up claiming a duration of data it does not have.

wstart/wend are what the trace holds. wstart_requested/ wend_requested are what was asked for. Callers that want a window length should use the former; callers reporting on the cut want the latter.

Must be called after the trim, or the two pairs are the same.

Parameters:

tr (Trace)
start (UTCDateTime)
end (UTCDateTime)

Return type:

None

specmod.preprocess.get_sta_shift(sta, sta_shift)[source]#

The per-station timing correction for sta, or zero.

sta_shift maps station name to a shift in seconds, e.g. {"STA": 0.5}.

Parameters:

sta (str)
sta_shift (Mapping[str, float] | None)

Return type:

float

specmod.preprocess.cut_p(st, bf=0, tafp=0.8, time_after='relative_time', sta_shift=None, refine_window=False)[source]#

Function to cut a p wave window from an Obspy trace obeject

bf (int/float) time shift in seconds before the P-wave arrival time

raf (int/float) ratio of p-s time to fix the end of the P-window

sta_shift (dict) dictionary of station names and station specific time shifts in seconds

refine_window (bool) True if you want to use squared intergral percentiles to refine the signal window.

Parameters:

st (Stream)
bf (float)
tafp (float)
time_after (str)
sta_shift (Mapping[str, float] | None)
refine_window (bool)

Return type:

None

specmod.preprocess.cut_s(st, rafp=0.8, tafs=20, time_after='absolute_time', sta_shift=None, refine_window=True)[source]#

Function to cut a s wave window from an Obspy trace obeject.

bf (int/float) time shift in seconds before the P-wave arrival time

rafp (int/float) ratio of p-s time to fix the start the of S-window

tafs (int/float) window length in seconds or scaling factor of relative p-s time

time_after (str) can be set to ‘absolute_time’ or ‘relative_ps’

if time_after == ‘absolute_time’ the window length is given as a value in seconds

if time_after == ‘relative_ps’ the value should be some number that scales with the p-s differential time

sta_shift (dict) dictionary of station names and station specific time shifts in seconds

Modified by Pungky Suroyo.

Parameters:

st (Stream)
rafp (float)
tafs (float)
time_after (str)
sta_shift (Mapping[str, float] | None)
refine_window (bool)

Return type:

None

specmod.preprocess.pad_traces(st, pad_len=1, pad_val=0)[source]#

Util to pad waveforms with zeros before and after the start and endtime of trace.

Parameters:

st (Stream)
pad_len (float)
pad_val (float)

Return type:

None

specmod.preprocess.cut_c(st, bf=2, raf=0.8, tafp=1.4, sta_shift=None)[source]#

Function to cut a coda wave window from an Obspy trace object

Written by Pungky Suroyo.

Parameters:

st (Stream)
bf (float)
raf (float)
tafp (float)
sta_shift (Mapping[str, float] | None)

Return type:

None

Spectra#

Waveforms in, SpectrumSet out.

The replacement for spectral.Spectra.from_streams, and the seam that spectral becomes a shell over. Given a signal stream and the noise cut to match it, this produces the immutable containers directly — trace to Spectrum to SpectrumPair — with none of spectral’s mutable classes in between.

Why it lives here rather than in ``core``. core and transforms know nothing about ObsPy: they operate on arrays, a duration and a sampling rate. That is worth keeping, because it is what makes them testable without constructing a Stream and usable on data that never came from one. This module is the single place that knows what a Trace is, so the dependency stops at one file rather than spreading through the containers.

The route is shorter than the legacy one, and identical in value. spectral.Spectrum converts the estimator’s output to a PSD, copies the arrays out, converts back to magnitude, then SNP wraps them in a core.Spectrum again — a round trip through two mutable objects that exists only because the legacy call sequence had it. Here the estimator’s spectrum is converted once, to the unfolded magnitude convention the pipeline reads Omega in, and handed straight to SpectrumPair.compare().

That the two agree is not asserted, it is measured: tests/test_pipeline.py runs both over the same 28 windows and all five estimators and requires the numbers to match to 1 part in 1e12.

specmod.pipeline.estimate_spectrum(data, delta, *, motion=Motion.VELOCITY, **kwargs)[source]#

Transform a record with whichever estimator the configuration names.

The bridge to specmod.transforms. Every estimator — FFT, Welch, multitaper, Prieto, quadratic, CWT — becomes available through specmod.config.TransformConfig, where before there was a hardcoded mtspec(data, delta, 3).

Returns a specmod.core.Spectrum, which carries its own units, so the caller no longer has to track how many times the record has been integrated or what amplitude convention is in force. The convention is the folded one, FAS; spectrum_from_trace() is what converts to the unfolded magnitude the rest of the pipeline reads Omega in.

Keyword arguments override the configured estimator’s parameters, which is how the **kwargs passthrough from the stream entry points works.

Lived in spectral until the direct pipeline existed, which put the one typed bridge to transforms inside the untyped legacy shell and made every caller of it untyped too.

Parameters:

data (NDArray[np.float64])
delta (float)
motion (Motion | str)
kwargs (Any)

Return type:

Spectrum

specmod.pipeline.pair_from_traces(signal_trace, noise_trace, *, compare=None, **kwargs)[source]#

Transform a signal and its noise, and pair them.

kwargs go to the estimator; compare overrides individual arguments to SpectrumPair.compare(), which otherwise come from configuration.

Parameters:

signal_trace (Any)
noise_trace (Any)
compare (Mapping[str, Any] | None)
kwargs (Any)

Return type:

SpectrumPair

specmod.pipeline.spectrum_from_trace(trace, **kwargs)[source]#

Transform one trace into a spectrum on the pipeline’s convention.

The convention is MAGNITUDE, the unfolded transform magnitude |X(f)|, not the folded 2|X| the estimators return. This is not a preference: Omega is defined on the unfolded spectrum — the long-period displacement plateau is |X(f -> 0)| = |int u dt| and M0 is proportional to it — so reading it off a folded spectrum puts every moment out by two, which is 0.2 magnitude units. See specmod.core.Spectrum.to_kind() for the conversion.

kwargs override the configured estimator’s parameters, including estimator itself.

Parameters:

trace (Any)
kwargs (Any)

Return type:

Spectrum

specmod.pipeline.spectrum_set_from_streams(signal, noise, *, event='', compare=None, **kwargs)[source]#

Pair two streams trace by trace.

The two are matched by position, as the legacy path matched them, since that is how preprocess.get_noise_p builds the noise: one trace per signal trace, in order. Mismatched ids are an error rather than a warning — a signal compared against another station’s noise is not a degraded measurement, it is a meaningless one.

event defaults to the origin time carried on the first trace, which is what the legacy container used to name the group.

Parameters:

signal (Iterable[Any])
noise (Iterable[Any])
event (str)
compare (Mapping[str, Any] | None)
kwargs (Any)

Return type:

SpectrumSet

The Spectrum container.

Immutable, self-describing, and normalisation-aware. Every operation returns a new object rather than mutating in place, so a spectrum cannot be silently integrated twice — the pre-refactor Spectrum.integrate() mutated, and its only inverse was differentiate(), which is neither exact nor recorded.

class specmod.core.spectrum.Spectrum(freq, amp, motion, kind, duration, sampling_rate, meta=<factory>)[source]#

Bases: object

A one-sided spectrum that knows its own units.

Parameters:

freq (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Frequency axis in Hz, strictly increasing, excluding DC by default.
amp (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Amplitude in whatever kind declares.
motion (specmod.core.units.Motion) – Ground-motion domain.
kind (specmod.core.units.AmplitudeKind) – What amp represents.
duration (float) – Physical record duration in seconds, n_samples * dt. Carried explicitly because every conversion between kinds needs it and it cannot be recovered from len(freq) once padding is involved.
sampling_rate (float) – Samples per second of the source record, in Hz.
meta (collections.abc.Mapping[str, Any]) – Arbitrary trace metadata. Stored read-only so a shared mapping cannot be mutated through one spectrum and observed through another.

property unit: str#: Unit string, e.g. m/s*s for a velocity FAS.

property n_samples: int#

Samples in the source record, duration * sampling_rate.

The third of the triple, derived rather than stored so it cannot disagree with the other two. Validated on construction — see _validate_record_geometry() for why that matters.

Note this is the record length, not len(freq). Zero-padding changes the second and not the first, and confusing them is the §2.2 bug.

property frequency_resolution: float#

1/T — the narrowest frequency difference the record can resolve.

This is the low-frequency floor the SNR bandwidth search must respect; nothing enforced it before, so a short window could report usable bandwidth below what it could physically resolve.

to_kind(kind)[source]#

Convert between FAS, MAGNITUDE, PSD and ASD.

Conversions go via FAS rather than being enumerated pairwise, so the factors of 2T and of the fold each live in exactly one place.

MAGNITUDE is the conversion to reach for when reading a long-period level: Omega is defined on |X|, not on the folded FAS, and the two differ by two. Asking for it by name is the point — the factor is easy to apply by hand and easy to apply twice, or not at all.

Parameters:: kind (AmplitudeKind | str)
Return type:: Spectrum

to_motion(motion)[source]#

Integrate or differentiate to another ground-motion domain.

Multiplies by (2*pi*f) per order of differentiation. Only valid on an amplitude-like kind, so a PSD is converted to FAS, transformed, and converted back — squaring the frequency factor would otherwise be silently wrong.

Parameters:: motion (Motion | str)
Return type:: Spectrum

band(fmin=None, fmax=None)[source]#

Restrict to a frequency band, inclusive of both bounds.

Parameters:

fmin (float | None)
fmax (float | None)

Return type:

Spectrum

energy()[source]#

Total signal energy, sum(x^2) * dt, recovered from the spectrum.

This is the quantity Parseval’s theorem ties to the time domain, and it is what the cross-estimator normalisation test asserts. For a one-sided FAS the two-sided integral folds to integral of A^2 / 2 df.

Return type:: float

Pairing a signal against its noise, and the band that survives the comparison.

This is the typed replacement for spectral.SNP and spectral.Spectra. The numerics are identical — tests/test_golden_reference.py holds both paths to the same 140 window-estimator results — but three structural properties change, and they are the reason the rewrite is worth doing.

Configuration is an argument, not an import-time global. spectral.py binds every setting at module import (BW_METHOD, ROT_METHOD and eight more). That is why a Brune and a Boatwright model cannot be fitted in one session, why tests cannot vary configuration without reimporting, and why they cannot run in parallel. Everything here takes its settings as parameters.

Nothing mutates. The legacy classes rescale, rotate, interpolate and integrate in place, which is what made core.Spectrum’s read-only arrays break the pipeline when the estimators were rewired: the containers were mutating arrays they did not own. Each step here returns a new object, so a spectrum cannot change under a reference someone else is holding.

The pieces are separable. The binning, the Parseval rescale, the interpolation and the band search are module-level functions over arrays. They were private methods reachable only by constructing a full pair from two obspy traces, so the only way to test the band search was to run the whole pipeline.

class specmod.core.collection.BinnedSpectrum(freq, amp)[source]#

Bases: object

A spectrum averaged into log-spaced bins.

Separate from Spectrum because it is not one: the bin centres are geometric midpoints of the edges rather than Fourier frequencies, so record geometry (duration, sampling_rate) no longer determines the axis and the Parseval contract does not hold on it. Conflating the two is how a binned spectrum ends up being handed to something that assumes an FFT grid.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
amp (ndarray[tuple[Any, ...], dtype[float64]])

class specmod.core.collection.FittableView(pair, id='')[source]#

Bases: object

A pair presented as the flat thing a fitter reads.

SpectrumPair keeps the unbinned spectrum and its binned form as separate objects, which is right for the comparison — they are different kinds of thing, and conflating them is how a binned axis ends up somewhere that assumes an FFT grid. A fitter wants them side by side, so this is the view that puts them there.

A view rather than a conversion: it holds the pair and reads through, so there is one copy of the arrays and no question of which is authoritative.

Parameters:

pair (SpectrumPair)
id (str)

property motion: Motion#

The ground-motion domain the arrays are in.

Read through like the rest, because a fitter has to know it: the initial guess for fc is the frequency of the spectral peak, and that is the corner only in velocity.

class specmod.core.collection.SpectrumPair(signal, noise, binned_signal, binned_noise, snr, resolution_floor, band=None, meta=<factory>)[source]#

Bases: object

A signal spectrum and the noise it is judged against.

Build with compare(), which runs the rescale, the interpolation, the binning and the band search in the order they depend on each other.

Parameters:

signal (Spectrum)
noise (Spectrum)
binned_signal (BinnedSpectrum)
binned_noise (BinnedSpectrum)
snr (ndarray[tuple[Any, ...], dtype[float64]])
resolution_floor (float)
band (tuple[float, float] | None)
meta (Mapping[str, Any])

SETTINGS_KEY: ClassVar[str] = 'compare_settings'#: Where compare() records its own arguments inside meta.

property passes: bool#: Whether a usable band survived.

for_fitting(id='')[source]#

This pair as the flat view a fitter reads. See FittableView.

Parameters:: id (str)
Return type:: FittableView

classmethod compare(signal, noise, *, threshold=3.0, f_min=0.001, f_max=200.0, n_bins=101, scale_parseval=True, resolution_floor=True, rotate_noise=True, noise_model='boost', bandwidth='peak', rotation_inc=0.05, rotation_space=(0.001, 1.001), meta=None)[source]#

Pair the two and select the band.

The order matters and is not arbitrary. The noise is rescaled and moved onto the signal’s frequency axis before binning, which is what makes the two binned arrays share bin edges — the element-wise ratio below is only meaningful because of it, and it holds for every estimator including those whose native axes differ in length.

The floor is captured from the two spectra before the interpolation, because afterwards the noise carries the signal’s axis and its own lowest resolvable frequency is unrecoverable.

Parameters:

signal (Spectrum)
noise (Spectrum)
threshold (float)
f_min (float)
f_max (float)
n_bins (int)
scale_parseval (bool)
resolution_floor (bool)
rotate_noise (bool)
noise_model (str | NoiseModel)
bandwidth (str)
rotation_inc (float)
rotation_space (tuple[float, float])
meta (Mapping[str, Any] | None)

Return type:

SpectrumPair

to_motion(motion)[source]#

This pair in another ground-motion domain, re-binned and re-banded.

Replaces spectral.SNP.integrate/differentiate, which mutated in place and had no way to express “the same event, as displacement” other than destroying the velocity one. This returns a new pair.

The noise is not lifted again. self.noise already carries the lift from the comparison that built this pair, and applying it a second time would compound on every conversion — narrowing the band each time. The pre-refactor code guarded this with a ROTATED flag; here it falls out of the settings being replayed with rotate_noise=False.

The band can move, and not for the reason it first appears. The unbinned signal-to-noise ratio is invariant under a domain change — both spectra are multiplied by the same power of 2*pi*f. The binned ratio is not, because a bin holds the geometric mean of log10(amp) and averaging log10(a/f) over a bin is not log10(a) averaged minus log10(f_centre) unless the centre is the geometric mean of the frequencies in it. Measured on the 28 PNR windows: the binned ratio moves by up to 16%, and 3 of the 28 bands with it.

Parameters:: motion (Motion | str)
Return type:: SpectrumPair

class specmod.core.collection.SpectrumSet(pairs, event='', meta=<factory>)[source]#

Bases: object

The pairs for one event, keyed by trace id.

Replaces spectral.Spectra. A mapping rather than a class with a group attribute, so the obvious operations — iterate, filter, count — are the ones that work.

Parameters:

pairs (Mapping[str, SpectrumPair])
event (str)
meta (Mapping[str, Any])

passing()[source]#

Only the pairs that yielded a usable band.

Return type:: SpectrumSet

to_motion(motion)[source]#

The whole event in another ground-motion domain.

Replaces spectral.Spectra.inte/diff. See SpectrumPair.to_motion() for what is recomputed and what is not.

Parameters:: motion (Motion | str)
Return type:: SpectrumSet

specmod.core.collection.find_bandwidth(freq, snr, threshold, *, method='peak')[source]#

Select the usable band with a named strategy.

A thin front for specmod.core.bandwidth.BANDWIDTH_SELECTORS. The default is "peak", which is what the shipped configuration has always used — the legacy BW_METHOD = 2. See that module for what the strategies assume and why the choice matters.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
snr (ndarray[tuple[Any, ...], dtype[float64]])
threshold (float)
method (str)

Return type:

tuple[float, float] | None

specmod.core.collection.interpolate_onto(target_freq, freq, amp)[source]#

Resample amp onto target_freq.

Warning

np.interp does not extrapolate — it repeats the edge value. Below freq.min() the result is therefore a flat continuation rather than a measurement, and a signal-to-noise ratio computed there has an invented denominator. SpectrumPair.resolution_floor() is what keeps the selected band out of that region; this function does not, and must not be used without it.

Parameters:

target_freq (ndarray[tuple[Any, ...], dtype[float64]])
freq (ndarray[tuple[Any, ...], dtype[float64]])
amp (ndarray[tuple[Any, ...], dtype[float64]])

Return type:

ndarray[tuple[Any, …], dtype[float64]]

specmod.core.collection.log_bin(freq, amp, *, f_min=0.001, f_max=200.0, n_bins=101)[source]#

Average amp into n_bins log-spaced bins, dropping empty ones.

The requested range is clamped to the record’s own, which is what makes the requested bin count the count you get. Unclamped, the shipped defaults (0.001 Hz to 200 Hz) sit far outside any real record — on the PNR data roughly a third of the bins fall below the lowest frequency present and a third above the highest, all of them empty — which is why the surviving axis was always far shorter than n_bins.

The average is geometric (the mean of log10(amp)), matching the log scale the bins themselves are spaced on. Empty bins are expected rather than exceptional — log bins over a linear grid are inevitably sparse at the low end — so they are dropped silently rather than warned about per bin.

Membership is computed, not tested. The bin index comes from the position of log10(f) along the range, which puts every sample in exactly one bin. The previous version tested f >= left and f <= right against each edge in turn: both ends closed, so a sample landing on an interior edge belonged to two bins, and which of the two comparisons succeeded depended on the last bit of np.logspace. That is one of the three places where a last-bit difference changed a result — it moved the surviving bin count by one, and with it the length of bsnr. Computing the index removes the double membership and the edge comparison together.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
amp (ndarray[tuple[Any, ...], dtype[float64]])
f_min (float)
f_max (float)
n_bins (int)

Return type:

BinnedSpectrum

specmod.core.collection.parseval_scale(n_signal, n_noise)[source]#

Factor putting a noise spectrum on the signal’s energy footing.

The two windows are rarely the same length — 1.2 to 1.6 s of noise against 1.8 to 3.5 s of signal on the PNR data — and a shorter record spreads the same power over fewer bins. Comparing them without this compares spectra computed over different durations.

Parameters:

n_signal (int)
n_noise (int)

Return type:

float

Physical units and domains, as types rather than conventions.

The pre-refactor code tracked none of this. A Spectrum did not know whether it held power or amplitude, nor whether it was in displacement, velocity or acceleration; that lived in Models.MOTION, a module global read at import time which the user had to keep in sync by hand with however many times they had called .inte() or .diff(). Getting it wrong returned a wrong seismic moment with no error anywhere.

Making both a typed attribute turns those silent factor errors into exceptions.

Conventions#

The canonical amplitude kind is the folded one-sided Fourier amplitude spectrum (AmplitudeKind.FAS), in units of [signal] * s. “Folded” means the negative-frequency half has been added in, so FAS = 2|X| where X is the Fourier transform. That is what makes energy recoverable by integrating over non-negative frequencies alone, and it is why every estimator here can be held to one Parseval check.

It is not the quantity the source model is written in. Omega, the long-period spectral level, is the plateau of |X| — at zero frequency |X(0)| = |integral u dt|, which is what M0 is proportional to. Reading FAS as Omega puts M0 out by two, which is 0.2 magnitude units. Use AmplitudeKind.MAGNITUDE for that, and let Spectrum.to_kind() apply the factor rather than doing it by hand.

Relationships between the kinds, for a record of duration T:

Kind	Units	From FAS `A`
`FAS`	`[x] * s`	–
`MAGNITUDE`	`[x] * s`	`\|X\| = A / 2`
`PSD`	`[x]^2 / Hz`	`P = A^2 / (2 T)`
`ASD`	`[x] / sqrt(Hz)`	`D = A / sqrt(2 T)`

Parseval takes a different form in each amplitude convention, which is the whole reason both are named here rather than left to the caller:

E = integral A**2 / 2 df        (FAS, folded)
E = 2 * integral |X|**2 df      (MAGNITUDE, unfolded)

T is the physical record duration, n_samples * dt. It is never inferred from the length of the frequency axis: zero-padding changes that length while leaving the duration alone, which is precisely how the old psd_to_amp acquired a padding-dependent error.

class specmod.core.units.AmplitudeKind(value)[source]#

Bases: StrEnum

What the amplitude axis of a spectrum represents.

FAS = 'fas'#: Folded one-sided Fourier amplitude spectrum, 2|X|, in [x] * s. Energy is integral(FAS**2 / 2) df. The default, because it is the convention in which one Parseval check covers every estimator.

MAGNITUDE = 'magnitude'#: Unfolded Fourier transform magnitude, |X| = |rfft(x)| * dt, in [x] * s. This is the one Omega is defined in, and the one to read a long-period spectral level off. Energy is 2 * integral(|X|**2) df.

PSD = 'psd'#: One-sided power spectral density, [x]^2 / Hz.

ASD = 'asd'#: One-sided amplitude spectral density, [x] / sqrt(Hz).

property is_amplitude: bool#

Whether this kind scales linearly with the record.

The distinction that matters for Spectrum.to_motion(): applying a 2*pi*f factor to a squared quantity is wrong by 2*pi*f again.

unit(motion)[source]#

Full unit string for this kind in a given motion domain.

Parameters:: motion (Motion)
Return type:: str

class specmod.core.units.Motion(value)[source]#

Bases: StrEnum

Ground-motion domain of a time series or spectrum.

property derivative_order: int#

Order of time differentiation relative to displacement.

Converting between domains multiplies the spectrum by (2*pi*f) per order, so the difference of two orders gives the exponent directly.

property unit: str#: SI unit of the time-domain signal.

The time-frequency surface a CWT produces, and the QC it makes possible.

A Scalogram is deliberately not an amplitude spectrum. Its power is |W(a,b)|**2 in the L2-Morlet convention, which carries units of [signal]**2 * time — documented, but not comparable to a Fourier amplitude spectrum and not something to fit a source model to.

The conversion happens in exactly one place, Scalogram.time_average(), which applies the C_delta and dj*dt bridge and returns an ordinary Spectrum. One normalisation path, one test. A second “already normalised” surface would be a second thing to get wrong.

References

Torrence, C. and Compo, G.P. (1998). A practical guide to wavelet analysis. Bulletin of the American Meteorological Society 79(1), 61-78.

class specmod.core.scalogram.Scalogram(time, freq, power, scales, coi, c_delta, dj, dt, motion, meta=<factory>)[source]#

Bases: object

Full time-frequency surface from a continuous wavelet transform.

Parameters:

time (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Sample times, shape (n_times,).
freq (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Fourier-equivalent frequencies, shape (n_scales,). These are true Fourier frequencies via the analytic Morlet relation, not scales, so the axis means the same thing as an FFT’s.
power (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – |W(a,b)|**2, shape (n_scales, n_times).
scales (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Wavelet scales in seconds, shape (n_scales,). Needed by the normalisation bridge, which divides by scale.
coi (numpy.ndarray[tuple[Any, ...], numpy.dtype[numpy.float64]]) – Longest resolvable period at each time, shape (n_times,). A frequency is inside the cone of influence where 1/freq > coi.
c_delta (float) – Reconstruction constant for the wavelet actually used. Computed rather than tabulated, so a non-default omega0 stays correct.
dj (float) – Spacing of the log-scale grid, in octaves.
dt (float)
motion (Motion)
meta (MappingProxyType)

coi_mask()[source]#

True where a coefficient is free of edge effects.

Return type:: ndarray[tuple[Any, …], dtype[bool]]

coi_coverage()[source]#

Fraction of the window free of edge effects, per frequency.

This is the number that says whether a window is long enough to constrain the low-frequency plateau, which is what sets Omega.

Return type:: ndarray[tuple[Any, …], dtype[float64]]

time_average(*, mask_coi=True)[source]#

Collapse to an ordinary amplitude spectrum.

Applies the Torrence & Compo normalisation so that the result satisfies the same Parseval contract as every other estimator: summing the wavelet power over scales, weighted by dj*dt/C_delta and divided by scale, returns the record’s energy.

Parameters:: mask_coi (bool) – Exclude coefficients inside the cone of influence and rescale by the surviving fraction. Without this a short window reads low at low frequency — precisely the band that constrains Omega.
Return type:: Spectrum

qc()[source]#

Compute the §4.4.2 checks.

Return type:: ScalogramQC

class specmod.core.scalogram.ScalogramQC(lowest_resolved_frequency, median_coi_coverage, temporal_concentration, half_window_ratio)[source]#

Bases: object

Quality checks a time-frequency surface makes possible.

An amplitude-only signal-to-noise test cannot see any of these: it collapses the time axis before looking. Every field is computed and recorded rather than acted on — a trace is never silently dropped, the numbers travel with the result so they can be filtered downstream.

Parameters:

lowest_resolved_frequency (float)
median_coi_coverage (float)
temporal_concentration (float)
half_window_ratio (float)

lowest_resolved_frequency: float#: Lowest frequency with usable coverage outside the cone of influence. Window length imposes this limit and nothing else in the pipeline enforces it, so a short window can otherwise report usable bandwidth where the transform has no support.

median_coi_coverage: float#: Fraction of the window free of edge effects, per frequency, summarised as the median across the band.

temporal_concentration: float#: Normalised Gini coefficient of energy over time, in [0, 1]. Near 0 is stationary; near 1 means essentially all the energy is in a handful of samples, which is a glitch rather than an arrival.

half_window_ratio: float#: Ratio of spectral energy in the first half of the window to the second. Far from 1 suggests coda contamination, a second arrival, or a window that started late.

to_dict()[source]#

Flat mapping, for landing in a results table as columns.

Return type:: dict[str, float]

Transforms and smoothing#

Pluggable time-to-frequency estimators.

All backends satisfy the contract in specmod.transforms.base, so one test suite pins every one of them.

specmod.transforms.ESTIMATORS: dict[str, type[SpectralEstimator]] = {'cwt': <class 'specmod.transforms.cwt.CWTEstimator'>, 'fft': <class 'specmod.transforms.fft.FFTEstimator'>, 'multitaper': <class 'specmod.transforms.multitaper.MultitaperEstimator'>, 'prieto': <class 'specmod.transforms.prieto.PrietoMultitaperEstimator'>, 'quadratic': <class 'specmod.transforms.quadratic.QuadraticMultitaperEstimator'>, 'welch': <class 'specmod.transforms.fft.WelchEstimator'>}#: Name -> estimator, for resolving TransformConfig.estimator.

class specmod.transforms.CWTEstimator(omega0=6.0, dj=0.125, mask_coi=True, f_min=None, f_max=None, name='cwt')[source]#

Bases: object

Continuous wavelet transform, time-averaged to an amplitude spectrum.

Produces both outputs from one transform: scalogram() returns the full time-frequency surface, and estimate() returns its time average as an ordinary Spectrum, so the fitting pipeline does not know the difference. The surface is what you look at when a fit comes out wrong.

Parameters:

omega0 (float) – Morlet central frequency, dimensionless. Trades time resolution against frequency resolution; 6 is the conventional choice and the value for which the analytic scale-frequency relation is usually quoted.
dj (float) – Scale spacing in octaves. Smaller resolves the frequency axis more finely at proportionally more cost. The normalisation carries an explicit dj factor, so the recovered energy does not depend on it.
mask_coi (bool) – Exclude the cone of influence from the time average. On by default: without it a short window reads low at low frequency, which is the band that constrains Omega.
f_min (float | None) – Frequency range to cover. Defaults span 1/T — the longest period the record can represent — up to the Nyquist frequency.
f_max (float | None) – Frequency range to cover. Defaults span 1/T — the longest period the record can represent — up to the Nyquist frequency.
name (str)

scalogram(data, dt, *, motion=Motion.VELOCITY, meta=None)[source]#

The full time-frequency surface.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)
motion (Motion | str)
meta (dict[str, Any] | None)

Return type:

Scalogram

estimate(data, dt, *, motion=Motion.VELOCITY, meta=None)[source]#

Time-averaged amplitude spectrum, interchangeable with the others.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)
motion (Motion | str)
meta (dict[str, Any] | None)

Return type:

Spectrum

class specmod.transforms.FFTEstimator(taper='tukey', taper_alpha=0.05, taper_correction='energy', n_fft=None, drop_dc=True, name='fft')[source]#

Bases: object

One-sided Fourier amplitude spectrum via numpy.fft.rfft().

Parameters:

taper (str) – Window applied before transforming. tukey with a small alpha suppresses edge discontinuities while leaving the body of the record untouched.
taper_alpha (float) – Window applied before transforming. tukey with a small alpha suppresses edge discontinuities while leaving the body of the record untouched.
taper_correction (Literal['energy', 'amplitude']) – "energy" preserves Parseval, "amplitude" preserves the peak of a coherent sinusoid. See specmod.transforms.base.
n_fft (int | str | None) –
Transform length: None for no padding, an integer, or a strategy — "fast" for the next efficiently-factorised length, "pow2" for the next power of two. See resolve_n_fft().

Padding refines the frequency grid without changing amplitude — the property the pre-refactor normalisation got wrong by keying off len(freq). So it buys two things and neither is leakage suppression, which is the taper’s job: it removes scalloping loss (36% worst case on a line falling between bins, unpadded), and it avoids the slow path for an awkward record length.

"fast" is the one to reach for. Cut windows are not round numbers — of the 28 PNR S-windows, 17 are odd and several are prime — and a prime length costs 1.77x across those. "pow2" is offered because it is what people expect, but it overshoots: numpy’s pocketfft handles 5-smooth lengths, so padding 65537 to 131072 does twice the work of padding it to 65610.
drop_dc (bool) – Discard the zero-frequency bin.
name (str)

class specmod.transforms.MultitaperEstimator(time_bandwidth=3.0, n_tapers=5, adaptive=True, center=False, center_edge_tolerance=0.05, normalize_to_variance=False, drop_dc=True, name='multitaper')[source]#

Bases: object

Multitaper spectrum estimate.

Parameters:

time_bandwidth (float) – The time-bandwidth product NW. Larger values reduce variance and leakage at the cost of frequency resolution. In the pre-refactor code this was the literal 3 passed positionally to mtspec, with no way to configure it.
n_tapers (int) – Number of DPSS tapers. Must not exceed 2*NW - 1, beyond which the tapers are poorly concentrated and add leakage rather than reducing variance; exceeding it raises rather than silently degrading.
adaptive (bool) –
Apply Thomson’s adaptive weighting. When False, tapers are averaged with equal weight.

On by default, because leakage suppression is the reason to reach for multitaper at all and flat weighting does not provide it. Measured on a 2 Hz line 10^6 times stronger than the background — a mild version of what a seismic spectrum does across its band — the recovered noise floor between 20 and 49 Hz sits 287x above the truth with flat weighting and 1.1x with adaptive. A Brune fit reads t* and f_c off exactly that high-frequency decay, so an inflated floor is not a cosmetic problem.

The cost is resolution: adaptive weighting downweights the higher-order tapers wherever leakage would dominate, so it uses fewer effective degrees of freedom and gives a noisier estimate in bands where the signal is strong. Turn it off for a well-conditioned record with little dynamic range, where the extra averaging is worth more than the leakage rejection.
center (bool) –
Circularly shift the record so its energy centroid sits mid-window before estimating.

This removes the position dependence entirely rather than reducing it: measured across start positions from 2% to 78%, the recovered energy ratio is identical to three decimal places once centred, and the spectrum matches the naturally-centred case exactly. It is legitimate because |FFT| is invariant under a circular shift, so the quantity being estimated does not change.

What remains after centring is the taper concentration itself — a compact centred transient still reads about 1.16x high with flat weighting. That bias is consistent rather than position-dependent, which matters: a consistent multiplicative bias cancels in any ratio (signal-to-noise, spectral ratios, relative amplitudes between stations) and can be calibrated, where a position-dependent one cannot.

Off by default because a circular shift wraps. It is safe when the window edges are quiet, and refuses when they are not — see center_edge_tolerance.
center_edge_tolerance (float) – Maximum amplitude at the wrap point, as a fraction of the record’s peak, before centring raises rather than introducing a discontinuity. A window whose coda is still strong at the end cannot be safely rolled.
normalize_to_variance (bool) –
Rescale the whole spectrum so it integrates to the record’s variance, as Prieto’s multitaper package does (mtspec.py: sscal = xvar / (sum(spec)*df)). mtspec wrapped the same lineage, so enable this when reproducing pre-refactor results.

Off by default because it costs you a diagnostic. With it on, Spectrum.energy() returns the input energy whatever else is wrong, so it can no longer tell you whether the normalisation is sound — the one check that would catch a mis-scaled spectrum always passes. Left off, energy() is a live measurement of your own record.

It is a calibration, not a derivation: a single scalar chosen to force the integral, which is why it cannot change spectral shape. Total power comes out right and the shape absorbs whatever error remains, so it does not make the long-period level position-independent — see the warning above, and center for something that does.
drop_dc (bool)
name (str)

class specmod.transforms.PrietoMultitaperEstimator(time_bandwidth=3.0, n_tapers=5, weighting='adaptive', n_fft=None, drop_dc=True, name='prieto')[source]#

Bases: object

Multitaper estimate via Prieto’s multitaper package.

Parameters:

time_bandwidth (float) – As for MultitaperEstimator.
n_tapers (int) – As for MultitaperEstimator.
weighting (Literal['adaptive', 'constant', 'eigenvalue']) – adaptive is Prieto’s default. constant weights the tapers equally; eigenvalue weights by the concentration ratios.
n_fft (int | None) – Transform length. None leaves MTSpec’s own default, which pads. Padding refines the frequency grid; duration stays the physical record length either way.
drop_dc (bool)
name (str)

confidence_interval(data, dt, *, motion=Motion.VELOCITY)[source]#

Jackknife 95% confidence bounds, as two FAS spectra.

The main reason to reach for this backend rather than the native one. The interval is log-symmetric about the estimate — measured at roughly 2.5x either way for nw=3, kspec=5.

Note

Upstream inverts the two bounds for a small fraction of bins (~2% in testing), so low <= high does not hold everywhere. Sort the pair per-bin if you need that guarantee.

Raises:

NotImplementedError – For weighting="constant". This is an upstream bug, not a limitation: in multitaper.utils.jackspec the degrees-of-freedom array se keeps shape (nfft, 1) under constant weighting, so t.ppf(...) * sqrt(var[:, 0]) broadcasts to (nfft, nfft) instead of elementwise and raises. adaptive and eigenvalue are unaffected.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)
motion (Motion | str)

Return type:

tuple[Spectrum, Spectrum]

f_test(data, dt)[source]#

Thomson’s F-test for periodic (line) components.

Returns (freq, F, p) on the one-sided axis. Useful for spotting instrumental or cultural tones that would otherwise be mistaken for source structure.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)

Return type:

tuple[ndarray[tuple[Any, …], dtype[float64]], ndarray[tuple[Any, …], dtype[float64]], ndarray[tuple[Any, …], dtype[float64]]]

class specmod.transforms.QuadraticMultitaperEstimator(time_bandwidth=3.0, n_tapers=5, adaptive=True, center=False, center_edge_tolerance=0.05, normalize_to_variance=False, drop_dc=True, name='quadratic')[source]#

Bases: object

Curvature-corrected multitaper spectrum estimate.

Parameters:

time_bandwidth (float) – As MultitaperEstimator. The correction scales with W = NW/N, so a larger time-bandwidth product means more smoothing to undo and a correspondingly larger correction.
n_tapers (int) – As MultitaperEstimator. The correction scales with W = NW/N, so a larger time-bandwidth product means more smoothing to undo and a correspondingly larger correction.
adaptive (bool) – Whether the eigencoefficients fed to the curvature fit are weighted by Thomson’s adaptive scheme. Applies to the input weights only; the quadratic step itself is unweighted.
normalize_to_variance (bool) –
Rescale so the spectrum integrates to the record variance, as mtspec and Prieto’s package do. Off by default, for the reason given in MultitaperEstimator.

Note that Prieto’s MTSpec.qiinv() applies this unconditionally, and to an already-renormalised input, so reproducing that path needs it on.
center (bool) – Circularly shift the record so its energy centroid sits mid-window. Same rationale and same wrap check as the ordinary estimator.
center_edge_tolerance (float)
drop_dc (bool)
name (str)

class specmod.transforms.SpectralEstimator(*args, **kwargs)[source]#

Bases: Protocol

Anything that turns a real record into a Spectrum.

property name: str#

Short identifier, recorded in the spectrum’s metadata.

Declared read-only so the frozen dataclasses implementing this protocol satisfy it; a plain mutable attribute would not.

estimate(data, dt, *, motion=Motion.VELOCITY, meta=None)[source]#

Estimate the spectrum of data sampled at interval dt.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)
motion (Motion | str)
meta (dict[str, Any] | None)

Return type:

Spectrum

class specmod.transforms.WelchEstimator(segment_length=None, overlap=0.5, taper='hann', drop_dc=True, name='welch')[source]#

Bases: object

Segment-averaged PSD via scipy.signal.welch(), returned as FAS.

Averaging reduces variance at the cost of frequency resolution, which is the right trade for a noise window. Note that duration remains the full record length: it is the physical property of the record, not of the segments, and every kind conversion depends on it.

Parameters:

segment_length (int | None)
overlap (float)
taper (str)
drop_dc (bool)
name (str)

specmod.transforms.get_estimator(name, **kwargs)[source]#

Construct an estimator by name.

Parameters:

name (str)
kwargs (object)

Return type:

SpectralEstimator

The estimator interface and the one place normalisation is defined.

Every backend — FFT, Welch, multitaper, and later the CWT — returns a Spectrum obeying the same contract, so a single test suite pins all of them. That is the point of the abstraction: the pre-refactor code called mtspec directly from Spectrum.__init__ and threaded backend keyword arguments through three layers of public API, so there was nowhere for a shared contract to live.

The contract#

Given a real record x of N samples at interval dt:

The returned spectrum is one-sided, spanning (0, f_Nyquist].
Its default kind is FAS, in [x] * s.
duration is N * dt, the physical record length. Normalisation is keyed off it and never off len(freq), so zero-padding changes resolution without changing amplitude.
Energy is preserved: spectrum.energy() recovers sum(x^2) * dt to within the accuracy of the estimator.

Taper correction#

A taper attenuates the record, and the correction depends on what you are measuring. The two are not interchangeable:

energy (default): Divide by sqrt(mean(w^2)), preserving total power. Parseval then holds exactly, which is what the shared test asserts, and it is the right choice for a transient — which is what a seismic arrival is.
amplitude: Divide by mean(w), preserving the peak of a coherent sinusoid so that it reads A0 * T. Right when measuring a monochromatic line.

For a Tukey taper with alpha=0.05 the two differ by well under a percent, but the choice is explicit rather than implied.

class specmod.transforms.base.SpectralEstimator(*args, **kwargs)[source]#

Bases: Protocol

Anything that turns a real record into a Spectrum.

property name: str#

Short identifier, recorded in the spectrum’s metadata.

Declared read-only so the frozen dataclasses implementing this protocol satisfy it; a plain mutable attribute would not.

estimate(data, dt, *, motion=Motion.VELOCITY, meta=None)[source]#

Estimate the spectrum of data sampled at interval dt.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)
motion (Motion | str)
meta (dict[str, Any] | None)

Return type:

Spectrum

specmod.transforms.base.make_window(kind, n, alpha=0.05)[source]#

Build a taper of length n.

Parameters:

kind (str)
n (int)
alpha (float)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

specmod.transforms.base.prepare_record(data, dt)[source]#

Validate and demean a record, returning (x, n_samples, duration).

Demeaning is unconditional: a non-zero mean puts all of its energy in the DC bin, which is then discarded, so leaving it in loses energy the Parseval check would report as a failure.

Parameters:

data (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
dt (float)

Return type:

tuple[ndarray[tuple[Any, …], dtype[float64]], int, float]

specmod.transforms.base.window_correction(window, correction)[source]#

Scale factor that undoes a taper’s attenuation.

See the module docstring for why there are two.

Parameters:

window (ndarray[tuple[Any, ...], dtype[float64]])
correction (Literal['energy', 'amplitude'])

Return type:

float

Spectral smoothing.

Separate from estimation, and lossy by design — see specmod.smoothing.base for why a smoothed spectrum no longer satisfies the Parseval contract that specmod.transforms guarantees.

specmod.smoothing.SMOOTHERS: dict[str, type[Smoother]] = {'konno_ohmachi': <class 'specmod.smoothing.konno_ohmachi.KonnoOhmachi'>, 'log_bins': <class 'specmod.smoothing.log_bins.LogBinner'>}#: Name -> smoother, for resolving SmoothingConfig.method.

class specmod.smoothing.KonnoOhmachi(bandwidth=40.0, count=1, normalize=True, name='konno_ohmachi')[source]#

Bases: object

Konno-Ohmachi smoothing at fixed log-frequency bandwidth.

Parameters:

bandwidth (float) – The b coefficient. Smaller values smooth more aggressively.
count (int) – Number of times to apply the window. More than one is occasionally used for very noisy spectra.
normalize (bool) – Normalise each window to unit area. Without this the smoothed amplitudes are biased low where the window is truncated at the edges of the frequency axis.
name (str)

class specmod.smoothing.LogBinner(f_min=None, f_max=None, n_bins=151, statistic='geometric', min_count=1, drop_empty=True, name='log_bins')[source]#

Bases: object

Average a spectrum into log-spaced frequency bins.

Parameters:

f_min (float | None) – Bin edges in Hz. None derives them from the record: 1/T for f_min and Nyquist for f_max. Deriving is the sensible default because it is exactly the band the record can represent.
f_max (float | None) – Bin edges in Hz. None derives them from the record: 1/T for f_min and Nyquist for f_max. Deriving is the sensible default because it is exactly the band the record can represent.
n_bins (int) – Number of bins.
statistic (Literal['geometric', 'mean', 'median']) – How to combine samples within a bin. geometric is the mean of log10(amp), which is what the pre-refactor code computed and the right choice for a quantity spanning orders of magnitude — an arithmetic mean over a decade of amplitudes is dominated by its largest member.
min_count (int) – Minimum samples for a bin to be kept when drop_empty is set.
drop_empty (bool) –
Drop bins holding fewer than min_count samples. Log bins over a linearly-spaced frequency axis are inevitably sparse at the low end — a bin is only reliably populated above roughly 1 / (2.3 * dlog10f * T) — so dropping is usually what you want for plotting or fitting.

Set False to keep a fixed-length axis with nan in the empty bins. Combined with explicit f_min/f_max that guarantees two spectra bin onto identical axes, which is what an element-wise signal-to-noise ratio requires.
name (str)

edges_for(spectrum)[source]#

Bin edges for a given spectrum.

Explicit f_min/f_max are honoured exactly and never clamped to the spectrum’s own range. That matters more than it looks: signal and noise windows have different durations, so clamping would bin them onto different axes, and the SNR ratio compares them element-wise. Pinning both edges is how a caller guarantees a shared axis.

Only derived bounds are clamped, since deriving already means “whatever this record supports”.

Parameters:: spectrum (Spectrum)
Return type:: ndarray[tuple[Any, …], dtype[float64]]

class specmod.smoothing.Smoother(*args, **kwargs)[source]#

Bases: Protocol

Anything that maps a spectrum to a smoothed spectrum.

property name: str#: Short identifier, recorded in the spectrum’s metadata.

smooth(spectrum)[source]#

Return a smoothed copy of spectrum.

Parameters:: spectrum (Spectrum)
Return type:: Spectrum

specmod.smoothing.get_smoother(name, **kwargs)[source]#

Construct a smoother by name.

Parameters:

name (str)
kwargs (object)

Return type:

Smoother

specmod.smoothing.is_smoothed(spectrum)[source]#

Whether any smoother has been applied.

Worth checking before comparing a spectrum’s energy against a time-domain expectation: smoothing does not preserve it.

Parameters:: spectrum (Spectrum)
Return type:: bool

Source models#

Source, attenuation and motion models, and the registry that names them.

The spectrum being fitted is a product of three things, so this package keeps them separate and composes them:

log10 A(f) = log10 S(f)  +  log10 D(f)  +  log10 G(f)
             source         attenuation    motion

Splitting them matters because they are chosen independently and confused easily. A Brune source with frequency-dependent attenuation is a legitimate combination; so is a Boatwright source with constant t*. The legacy code bound the pair together at import time through MODEL = which_model(...), which is why a Brune and a Boatwright could not be fitted in one session.

What a source model has to carry#

Not just a spectral shape. The relation between corner frequency and source radius belongs to the model, and forgetting that is a specific trap:

Madariaga is omega-squared like Brune and sits at the same (gamma, n), so adding it as a shape alone changes no fitted parameter — the fit is identical — while the source radius it implies is quite different. Stress drop goes as r**-3, so that is roughly an order of magnitude on identical data, arriving silently. SourceModel.corner_frequency_coefficient exists so the difference cannot be lost that way, and so that a model which changes nothing about the fit still changes what is derived from it.

class specmod.sources.AttenuationModel(*args, **kwargs)[source]#

Bases: Protocol

Anything that attenuates a source spectrum along the path.

property parameters: tuple[str, ...]#: Free parameter names, in the order the fitter should take them.

log10_decay(freq, *values)[source]#

log10 D(f), given this model’s parameters in order.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
values (float)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

class specmod.sources.BoatwrightSource(gamma=2.0, n=2.0)[source]#

Bases: _GeneralisedSource

Boatwright (1980), omega-squared with a sharper corner than Brune.

gamma = 2 narrows the transition; the high-frequency falloff is the same f**-2. The coefficients here are Brune’s, because Boatwright’s formulation shares the kinematic radius relation — recorded explicitly so that it is a stated choice rather than an omission.

Parameters:

gamma (float)
n (float)

class specmod.sources.BruneSource(gamma=1.0, n=2.0)[source]#

Bases: _GeneralisedSource

Brune (1970), omega-squared with a smooth corner.

The corner-frequency coefficients are Brune’s own kinematic values. They are not interchangeable with Madariaga’s — see the note in specmod.sources.

Parameters:

gamma (float)
n (float)

class specmod.sources.ConstantQ[source]#

Bases: object

Frequency-independent t*.

\[\log_{10} D(f) = -\pi f t^{*} / \ln 10\]

The ln 10 is the conversion into base-10 logs, which is the space the fit is performed in. Getting it wrong scales t* by 2.3 and leaves the spectrum looking plausible.

class specmod.sources.FrequencyDependentQ[source]#

Bases: object

t* with a power-law frequency dependence.

\[\log_{10} D(f) = -\pi f^{1-a} t^{*} / \ln 10\]

a = 0 recovers ConstantQ. The two are separate models rather than one with a switch because they have different free parameters, and a fitter needs to know that before it starts.

class specmod.sources.SourceModel(*args, **kwargs)[source]#

Bases: Protocol

A source spectral shape, plus what it implies about the source.

The second half is the part that is easy to leave out. Two models can share a spectral shape exactly and still disagree by an order of magnitude on stress drop, because they disagree about what a given corner frequency says about the rupture dimension.

property name: str#: Short identifier, as configuration refers to it.

property corner_frequency_coefficient: tuple[float, float]#

(k_P, k_S) in f_c = k * beta / r.

The bridge from a fitted corner frequency to a source radius, and therefore to stress drop. It belongs to the model rather than to whatever computes stress drop, because models that share a spectral shape do not share this.

log10_shape(freq, log10_omega, f_c)[source]#

log10 S(f) for a plateau log10_omega and corner f_c.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
log10_omega (float)
f_c (float)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

class specmod.sources.SpectralModel(source, attenuation, motion)[source]#

Bases: object

A source, an attenuation law and a motion, evaluated together.

\[\log_{10} A(f) = \log_{10} S(f) + \log_{10} D(f) + \log_{10} G(f)\]

Immutable and self-describing: it knows its own parameter names, so a fitter does not have to be told them separately and cannot be told them wrongly. That is what the legacy could not do — FitSpectra.set_model took a bare function, and the parameter names came from introspecting its signature, so the model and its parameters could disagree.

Parameters:

source (SourceModel)
attenuation (AttenuationModel)
motion (Motion)

property parameters: tuple[str, ...]#: Free parameters, in the order evaluate() takes them.

evaluate(freq, *values)[source]#

log10 A(f) for the given parameter values, in order.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
values (float)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

describe()[source]#

One line naming every choice this model makes.

Return type:: str

as_callable()[source]#

A plain function of (f, llpsp, fc, ...), for lmfit.Model.

lmfit discovers parameter names with inspect.signature() (checked, rather than assumed — it does not read co_varnames), so attaching a __signature__ is enough and no code has to be generated. A wrapper taking *args would give the fit a single positional parameter and nothing to vary.

Return type:: Callable[[…], ndarray[tuple[Any, …], dtype[float64]]]

specmod.sources.build_model(*, source='brune', motion=Motion.VELOCITY, frequency_dependent_attenuation=False)[source]#

Assemble a SpectralModel from configuration names.

This is the join that did not exist. ModelConfig.source was a Literal["brune", "boatwright"] that nothing read: FitSpectra took the model function as an argument, so the caller supplied Brune or Boatwright by hand and the configured value was decorative. Anything reading configuration now comes through here.

Parameters:

source (str)
motion (Motion | str)
frequency_dependent_attenuation (bool)

Return type:

SpectralModel

specmod.sources.from_config()[source]#

The model the current configuration asks for.

One call, so that [model] source = "boatwright" in a study file is what decides the shape being fitted — which it was not before.

Return type:: SpectralModel

specmod.sources.get_attenuation_model(name)[source]#

Resolve a registered attenuation model by name.

Parameters:: name (str)
Return type:: AttenuationModel

specmod.sources.get_source_model(name)[source]#

Resolve a registered source model by name, with its defaults.

Parameters:: name (str)
Return type:: SourceModel

specmod.sources.motion_scaling(freq, motion)[source]#

log10 G(f): the factor taking displacement to motion.

Source models are written for displacement, because that is where Omega and hence M0 are defined. Records are usually velocity. Each order of differentiation multiplies by 2 pi f, so in log space it is that many copies of log10(2 pi f).

The order comes from specmod.core.units.Motion rather than a string comparison, so an unrecognised motion fails at the enum instead of silently returning zero — which would fit a displacement model to velocity data and put Omega out by a factor of 2 pi f.

Parameters:

freq (ndarray[tuple[Any, ...], dtype[float64]])
motion (Motion | str)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

Fitting#

Fitting a source model to one spectrum, and to a whole event.

fittable_signal() decides what to fit, initial_guess() where to start, FitSpectrum fits one station and FitSpectra an event.

class specmod.fitting.FitSpectra(spectra, model=None, guess=None, fit_bins=None)[source]#

Bases: object

Fit every passing station in an event.

Parameters:

spectra (Any)
model (Any)
guess (Mapping[str, Mapping[str, float]] | None)
fit_bins (bool | None)

spectra: Any#: Declarations, as on FitSpectrum. models = {} at class level was one dictionary shared by every FitSpectra ever built; __init__ rebinds it, so nothing reached the shared copy, but nothing prevented it either. guess = {} was never assigned anywhere at all — a class attribute recording a constructor argument that is not kept.

fit_spectra(weight_method=None, **kwargs)[source]#

Fit every station, with the configured minimiser unless told otherwise.

method and weight_method both come from [fitting] when not given. Neither used to: fit_spectra() fell through to lmfit’s default minimiser, so a study file saying method = "powell" was ignored and the caller had to remember fit_spectra(method="powell") — which the tutorial does and nothing enforced.

It matters. On the 28 PNR windows lmfit’s default returns a negative corner frequency on one station where Powell does not; a corner frequency below zero is not a degraded measurement but a meaningless one, and nothing downstream rejects it.

Parameters:

weight_method (str | None)
kwargs (Any)

Return type:

None

init_fitting(model, guess, fit_bins)[source]#

Build a fit per passing station.

model=None resolves through the configuration once per station, which is cheap and keeps every fit in a run agreeing on what it is fitting.

Parameters:

model (Any)
guess (Mapping[str, Mapping[str, float]])
fit_bins (bool)

Return type:

None

reset(name='all')[source]#

Unbind every parameter, on one station or all of them.

The lookup tested name.upper() for membership and then indexed with name, so any id not already upper-case passed the check and raised KeyError on the next line. Station ids are upper-case in practice, which is why it never fired.

Parameters:: name (str)
Return type:: None

static write_flatfile(path, fits)[source]#

Write the group fit table, in the format path’s suffix names.

.parquet is typed, compressed and queryable without loading; .csv is what journal supplements want. See specmod.tables.

The previous implementation was os.makedirs(os.path.join( *path.split("/")[:-1])), which raised TypeError: join() missing 1 required positional argument for any path without a directory component — write_flatfile("out.csv", fits) could not work. It also split on / literally, so it did nothing useful on Windows.

Parameters:

path (str | Path)
fits (FitSpectra)

Return type:

Path

static read_flatfile(path)[source]#

Read a fit table back. Format follows the suffix.

Parameters:: path (str | Path)
Return type:: pd.DataFrame

class specmod.fitting.FitSpectrum(signal, model=None, fit_bins=False, **params)[source]#

Bases: object

Fit a source model to one spectrum with lmfit.

Takes anything carrying REQUIRED_SPECTRUM_ATTRIBUTES — in practice a FittableView from fittable_signal().

Parameters:

signal (Spectrumish)
model (Any)
fit_bins (bool)
params (float)

sig: Spectrumish#: Declarations, not defaults. These were class attributes carrying None and {}, which meant two things at once: every read had to cope with a None that __init__ had in fact replaced, and meta = {} was one dictionary shared by every instance ever constructed. __init__ assigns all of them, so the type is what it is after construction — and the shared-mutable-default hazard is gone rather than merely unreached.

result: lm.ModelResult | None#: None until fit_mod() runs. This one really is optional, and callers test it — see specmod.plotting.plot_pair().

spectral_model: sources.SpectralModel | None#: The specmod.sources.SpectralModel behind the fit, when there is one. None if a bare callable was supplied.

fit_mod(**kwargs)[source]#

Fit, judge the result, then record it — in that order.

The judgement used to be made after the recording, so the pass_fitting column of every flat file held the value from before the fit ran — True, the class default, on a fresh FitSpectrum. The attribute and the table disagreed, and the table is what gets written out and regressed on.

Parameters:: kwargs (Any)
Return type:: None

set_model(model=None, **params)[source]#

Set the model to fit.

Accepts a specmod.sources.SpectralModel, a bare callable, or None — in which case the model is whatever [model] in the configuration asks for. That default is the point: before it existed, config.model.source was read by nothing and the caller had to pass the right function by hand, so a study file saying source = "boatwright" silently got Brune.

A bare callable still works, because fitting an ad-hoc shape is a legitimate thing to want. It simply carries no provenance: spectral_model is None and nothing can report what was fitted.

Parameters:

model (Any)
params (float)

Return type:

None

property fitted: lm.ModelResult#

The fit result, or a message saying it has not been fitted.

Every private reader below went straight through self.result, which is None until fit_mod() runs — so calling quick_vis() on an unfitted spectrum raised AttributeError: 'NoneType' object has no attribute 'best_fit', from a line that names neither the station nor the missing step.

describe_model()[source]#

What is being fitted, or None for a bare callable.

Return type:: str | None

specmod.fitting.SpectraLike#: alias of Any

specmod.fitting.Spectrumish#: alias of Any

specmod.fitting.fittable_signal(pair, id='')[source]#

The signal to fit from a paired spectrum, or None to skip it.

Skipping is a decision the container should not have to spell out at every call site: a pair is unfittable when the signal-to-noise gate rejected it.

What comes back for a SpectrumPair is its FittableView, not its signal. The pair keeps the unbinned and binned spectra as separate objects, which is right for the comparison and wrong for a fitter that wants freq, amp, bfreq and bamp side by side; the view is what puts them there. id names the station on it, since a frozen pair does not carry one of its own.

The getattr fallback below is what a spectrum-like object that is not a pair takes — a bare view, or anything else presenting the same attributes. It is not a legacy shim; it is what lets the fitter be given something constructed by hand.

Parameters:

pair (Any)
id (str)

Return type:

Any | None

specmod.fitting.initial_guess(spectra, model=None)[source]#

Starting parameters for every fittable spectrum in spectra.

Replaces model_guess.create_simple_guess and its _fdep twin, which were two near-identical functions differing only in whether they added an a for frequency-dependent Q — so adding a third model meant writing a third guess function, and picking the wrong one gave lmfit a parameter the model did not take.

Which parameters are needed is asked of the model, not assumed. The fitted callable declares them in its signature, so a model gets exactly the guesses it takes and nothing else. Values that cannot be read off the spectrum come from [fitting] in the configuration.

The two that are read off the spectrum:

llpsp: log10 of the largest amplitude inside the selected band — the long-period plateau, which is what Omega is.
fc: the frequency at which that maximum falls.

Both assume a velocity spectrum, which is where a fit belongs anyway: the model carries a motion factor, so llpsp is the displacement plateau whichever domain is fitted, but converting first is not a neutral change of view — integrating implicitly low-passes and differentiating amplifies high-frequency noise, so the record to fit is the one the sensor recorded.

In velocity the peak is not merely near the corner, it is the corner, for any omega-squared source: the stationary point of f * [1 + (f/fc)**(gamma*n)]**(-1/gamma) sits at f = fc whenever n == 2, whatever the corner sharpness. In displacement and acceleration the spectrum is monotonic across the band, so the peak is whichever band edge it was handed and the guess is meaningless. Handed one of those, this warns rather than proceeding quietly.

Stations with no band are omitted rather than given None guesses. The old version emitted {"llpsp": None, "fc": None, "ts": None} on IndexError, which lmfit cannot use — the failure simply moved to the fit call.

Parameters:

spectra (Any)
model (Any)

Return type:

dict[str, dict[str, float]]

specmod.fitting.plot_columns()[source]#

How many columns a multi-panel figure uses, from [viz].

A function rather than a constant, and that is the whole point. This was PLOT_COLUMNS = cfg.load_config().config.viz.plot_columns evaluated at import time, so importing specmod.fitting resolved configuration against whatever directory the process happened to start in and froze the answer for the life of the interpreter. Measured: importing from a project whose specmod.toml says 5, then moving to one that resolves to 3, left the constant at 5 — a worker serving two projects would use the first one’s layout for both.

One home for the setting either way: it used to be defined in both the SPECTRAL and FITTING dicts, and the two copies could disagree.

Return type:: int

specmod.fitting.selected_band(spectrum)[source]#

The band to fit over, or None to fit everything available.

None rather than an empty array, because “no band survived” and “a band from 0 to 0” are different claims and the legacy spelling — an empty ubfreqs — could be read as either.

Parameters:: spectrum (Any)
Return type:: tuple[float, float] | None

The two-stage event fit, and the channel selection that feeds it.

Fitting a source model to one spectrum is not a unique inversion. The source corner and the path attenuation trade off against each other on the falling limb, and two minimisers can reach the same reduced chi-squared at corner frequencies differing by tens of percent — a factor of several in stress drop, which scales as fc**3. On the 28 PNR windows, Powell and leastsq land at 21.26 Hz and 14.75 Hz on one station at redchi 0.0259 against 0.0254.

That is not resolvable from one spectrum by any minimiser, and the published workflow does not try. f_c belongs to the source: every station sees the same rupture, so there is one value of it for the event. t* belongs to the path, and every station has a different one. So a station whose t* came out too high returns a corner that is too high, and the next station’s error does not point the same way. Averaging over the ensemble is not cosmetic smoothing — it uses the fact that the quantity being averaged is common to all of them while the contaminating one is not.

The two stages#

Omega, f_c and t* free at every station independently. The output is not the answer; it is N noisy estimates of one number plus N estimates of N different numbers.
The event f_c — a weighted mean over the stations that survived selection — is held fixed, and every station refits Omega and t* against a corner it can no longer trade against.

Measured on the same 28 windows, with the same 28 channels contributing to both: the two minimisers differ by a factor 1.44 in f_c at the worst station and by 0.4% on the event value. After stage two they agree to 0.23% on t* and 1.7e-3 log10 units on Omega.

“With the same channels contributing” is load-bearing, and is the trap in this module. See ChannelSelection.require_pass — comparing two minimisers under the default selection compares two different ensembles, and gives 125% rather than 0.4%.

Be clear about what that last number is not. Fixing f_c removes the parameter the minimisers were disagreeing about, so of course they then agree. What it shows is that the residual two-parameter problem is well conditioned: once the corner is pinned, Omega and t* are determined by the spectrum rather than negotiable. The judgement is concentrated into one number for the whole event, and that number came from the ensemble.

Choosing which channels contribute#

Selection is the part that cannot be automated away, because it is where quality control enters. A station with a bad instrument response, a clipped record or a pick on the wrong phase produces a corner frequency that is confidently wrong, and averaging it in moves the event value for every other station.

So the ensemble is chosen by ChannelSelection, which reads [fitting] and can be overridden per call. The order is fixed and each step is recorded with a reason, so StagedFit.excluded says why any given channel is not contributing:

anything not matching include (when include is non-empty),
anything matching exclude,
anything whose stage-1 fit failed its bounds, when require_pass,
anything left with no fit at all.

Patterns are shell globs, and they match at whichever level you write them. A trace id is NET.STA.LOC.CHA, and a pattern is tried against the whole id, against NET.STA, and against each component on its own — so all of these do what they look like they do:

"AQ07"              every channel of that station
"UR"                every station of that network
"UR.AQ07"           that station, spelled unambiguously
"HHE"               every east component, at every station
"HH?"               every high-gain broadband channel
"UR.AQ07.00.HHE"    exactly that channel
"LV.L00[123]..HH?"  a glob over the full id

Writing the station code alone is the common case after quality control — a clipped record or a bad response is a property of the instrument, not of one component — and needing "UR.AQ07.*" for it is the kind of detail that gets mistyped as "AQ07" and silently matches nothing. Station and channel codes do not collide in practice; where a pattern could match at two levels it matches, and StagedFit.excluded records which level it hit.

specmod.staged.WEIGHT_MODELS: dict[str, Any] = {'inverse_distance': <class 'specmod.staged.InverseDistance'>, 'inverse_epicentral_distance': <function <lambda>>, 'inverse_hypocentral_distance': <function <lambda>>, 'inverse_variance': <class 'specmod.staged.InverseVariance'>, 'uniform': <class 'specmod.staged.Uniform'>}#: Registered weightings, resolved by name from [fitting] event_weighting.

class specmod.staged.ChannelSelection(include=(), exclude=(), require_pass=True)[source]#

Bases: object

Which channels contribute to the event value.

Defaults come from [fitting]; pass one of these to override per call. The point of it being a value rather than four arguments is that a selection can be written down, compared and stored — a run that dropped three stations should be able to say so afterwards.

Parameters:

include (tuple[str, ...])
exclude (tuple[str, ...])
require_pass (bool)

require_pass: bool#

Drop a station whose stage-1 fit ended with a parameter against a bound.

This interacts with the minimiser, and not symmetrically. pass_fitting asks whether value +/- stderr reaches a bound. Powell — the shipped default — estimates no covariance matrix, so stderr is None, the spread is zero, and the test degenerates to “is the value exactly on the bound”, which essentially never fires. On these 28 windows it drops 0 stations under Powell and 6 under leastsq.

So changing the minimiser changes which stations vote, not just how each one is fitted. Comparing two minimisers under the default selection therefore compares two different ensembles: 12.081 Hz from 28 channels against 5.376 Hz from 22, a 125% difference that says almost nothing about the minimisers. Holding the ensemble fixed with require_pass=False gives 12.081 against 12.029 — 0.4%, which is the honest comparison and the number quoted in the module docstring.

Left True because the flag is doing the right thing when it fires: a pinned parameter reports its bound rather than a measurement, and averaging that in averages in a constant. But a study comparing minimisers, or reporting a corner frequency alongside one fitted another way, has to set this False or account for it.

choose(fit)[source]#

Split fit’s stations into contributors and reasoned exclusions.

Parameters:: fit (FitSpectra)
Return type:: tuple[list[str], dict[str, str]]

class specmod.staged.StagedFit(stage1, stage2, parameter, value, weighting, contributing, excluded=<factory>)[source]#

Bases: object

The result of fit_event(): both stages and how they were joined.

Parameters:

stage1 (FitSpectra)
stage2 (FitSpectra | None)
parameter (str)
value (float)
weighting (str)
contributing (tuple[str, ...])
excluded (Mapping[str, str])

stage1: FitSpectra#: Every station fitted independently. Kept because the spread across it is the evidence for how well constrained the event value is, and discarding it would leave only a number with no error on it.

stage2: FitSpectra | None#: Every station refitted with the event parameter held fixed. None when no station survived selection, in which case there is nothing to fix it to and the honest result is stage one alone.

parameter: str#: The aggregated value and what produced it.

property table: Any#: The stage-2 table where there is one, else stage one’s.

spread()[source]#

How much the contributing stations disagreed, before aggregation.

The number worth reporting beside the event value. A 2% spread and a 60% spread give the same weighted mean and mean very different things.

Return type:: dict[str, float]

describe()[source]#

One paragraph a caller can print or paste into a log.

Return type:: str

class specmod.staged.WeightModel(*args, **kwargs)[source]#

Bases: Protocol

How much each station’s stage-1 estimate counts toward the event value.

Given the stage-1 table and the spectra it was fitted to, return one weight per row. Weights need not sum to anything; they are normalised on use.

specmod.staged.fit_event(spectra, *, model=None, guess=None, fit_bins=None, parameter=None, weighting=None, selection=None, quiet=True, **fit_kwargs)[source]#

Fit an event in two stages, with the ensemble deciding the source term.

Every argument has a configured default, so fit_event(spectra) is the published workflow and the rest is there for when the default is wrong.

Parameters:

spectra (Any) – A SpectrumSet, or anything mapping trace ids to paired spectra.
model (Any) – Passed to FitSpectra for both stages, so the two are fitting the same thing.
guess (Mapping[str, Mapping[str, float]] | None) – Passed to FitSpectra for both stages, so the two are fitting the same thing.
fit_bins (bool | None) – Passed to FitSpectra for both stages, so the two are fitting the same thing.
parameter (str | None) – Which parameter the ensemble determines and stage two holds fixed. [fitting] event_parameter, default "fc" — the corner frequency is the term that belongs to the source. "ts" is the meaningful alternative, for a study with an independent handle on Q.
weighting (str | WeightModel | None) – A registered name or a WeightModel. [fitting] event_weighting, default inverse hypocentral distance.
selection (ChannelSelection | None) – Which channels contribute. Defaults to ChannelSelection.from_config().
quiet (bool) – Suppress the fitter’s per-station chatter, which is two full passes of it here. The failures are on StagedFit.excluded either way.
**fit_kwargs (Any) – Passed to fit_spectra() for both stages — method= most usefully.

Return type:

StagedFit

Notes

Stage two is skipped, with stage2=None, when selection leaves nothing. Fixing the parameter to a mean of no stations would be inventing the very number the second stage exists to constrain, and a caller who asked for two stages and got one should be able to see that rather than read a value.

specmod.staged.get_weight_model(name)[source]#

Resolve a registered weighting by name, with its defaults.

Parameters:: name (str)
Return type:: WeightModel

Magnitude#

Seismic moment and moment magnitude from a fitted long-period plateau.

The final step of the Edwards et al. (2010) spectral method. The fit supplies Omega, the long-period displacement plateau in m s — llpsp is its base-10 logarithm — and each channel supplies a distance.

M0 = 4 pi rho beta**3 R_0 (Omega / G(R)) / (Theta F) [N m] Mw = (2/3) (log10 M0 - 9.1)

Units#

rho in kg/m^3, beta in m/s, R_0 in metres, Omega in m s give M0 in newton-metres. G(R) is specmod.spreading, whose distance is in kilometres. These are two different distances and both units are correct; they cancel only at spreading exponent 1.

Constants#

Defaults are the S-wave values of Holt (2019) Ch. 1 §1.4 and Ch. 2 §2.2, all taken at the source:

`density`	2700	kg/m^3
`velocity`	3500	m/s, cubed here
`radiation_pattern`	0.55	average SH pattern over the focal sphere
`free_surface`	2	vertically incident SH
`reference_dist...`	1000	metres

0.55 is the SH average and pairs with free_surface = 2; the textbook 0.63 is the RMS over total S, a different quantity. No partition factor is applied.

Use seismic_moment() and moment_magnitude() for arrays, station_moments() to append m0 and mw to a fit table, and event_magnitude() for an event value aggregated over stations.

The absolute calibration is unverified — see §4.7 of docs/REFACTOR_PLAN.md.

class specmod.magnitude.MediumConstants(density=2700.0, velocity=3500.0, radiation_pattern=0.55, free_surface=2.0, reference_distance_m=1000.0)[source]#

Bases: object

Medium properties at the source, and the geometry of the measurement.

Values are those at the rupture, not at the station and not a crustal average. velocity enters cubed.

Parameters:

density (float)
velocity (float)
radiation_pattern (float)
free_surface (float)
reference_distance_m (float)

density: float = 2700.0#: kg/m^3.

velocity: float = 3500.0#: m/s, shear-wave velocity at the source.

radiation_pattern: float = 0.55#: Average SH radiation pattern over the focal sphere (Boatwright, 1978).

free_surface: float = 2.0#: Free-surface factor; 2 for vertically incident SH.

reference_distance_m: float = 1000.0#: Metres. The distance at which the source spectrum is defined.

class specmod.magnitude.MomentMagnitude(value, m0, stations, station_magnitudes, excluded=<factory>, constants=<factory>, spreading_model='power_law', distance_measure='rhyp', unit='Mw')[source]#

Bases: object

An event magnitude and the station estimates behind it.

Per-station values are kept so spread() can report how well constrained the event value is.

Parameters:

value (float)
m0 (float)
stations (tuple[str, ...])
station_magnitudes (tuple[float, ...])
excluded (dict[str, str])
constants (MediumConstants)
spreading_model (str)
distance_measure (str)
unit (str)

unit: str = 'Mw'#: Always "Mw". Carried so the number cannot be read without its scale.

spread()[source]#

How much the stations disagreed, in magnitude units.

Return type:: dict[str, float]

specmod.magnitude.event_magnitude(staged, *, constants=None, spreading=None, distance_measure='rhyp', outlier_sigma=2.5, min_stations=3)[source]#

Event Mw: the mean of station estimates, after rejecting outliers.

Follows Holt (2019) §2.2 — the sample mean of station magnitudes, excluding anything beyond outlier_sigma standard deviations, and no event value at all on fewer than min_stations.

Averaging in magnitude rather than in M0 makes this a geometric mean of moments. Rejection is a single pass, not iterated to convergence.

Parameters:

staged (Any)
constants (MediumConstants | None)
spreading (SpreadingModel | None)
distance_measure (str)
outlier_sigma (float)
min_stations (int)

Return type:

MomentMagnitude

specmod.magnitude.moment_magnitude(m0)[source]#

Hanks and Kanamori (1979), for m0 in N m.

Parameters:: m0 (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
Return type:: ndarray[tuple[Any, …], dtype[float64]]

specmod.magnitude.seismic_moment(omega, distance_km, *, constants=None, spreading=None)[source]#

Scalar seismic moment in N m, from plateaus and their distances.

omega is the long-period displacement plateau in m s — the fit’s llpsp is its base-10 logarithm. distance_km is source-to-site, in kilometres, one per plateau.

Parameters:

omega (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
distance_km (_Buffer | _SupportsArray[dtype[Any]] | _NestedSequence[_SupportsArray[dtype[Any]]] | complex | bytes | str | _NestedSequence[complex | bytes | str])
constants (MediumConstants | None)
spreading (SpreadingModel | None)

Return type:

ndarray[tuple[Any, …], dtype[float64]]

specmod.magnitude.station_moments(table, *, constants=None, spreading=None, distance_measure='rhyp', plateau_column='llpsp')[source]#

Per-station M0 and Mw from a fit table.

Takes a table rather than a StagedFit, so either stage or a flatfile read from disk all work. Returns a copy with m0 and mw appended; the input is not modified.

Parameters:

table (Any)
constants (MediumConstants | None)
spreading (SpreadingModel | None)
distance_measure (str)
plateau_column (str)

Return type:

Any

Geometric spreading models, as a registry.

Each model returns the dimensionless amplitude ratio between the reference distance and the site, so an observed plateau is corrected to the source by dividing by it. Distances are in kilometres, as published spreading tables are written.

PowerLaw is the default at exponent=1, the theoretical body-wave value. Piecewise takes the contiguous segments regional models are published as, Tabulated interpolates a supplied curve in log-log space, and HOLT_2019_UTAH is the refined Utah model of Holt (2019) Table 2.1.

Register new models in SPREADING_MODELS; build one by name with get_spreading_model(). Why the default is not fitted, and why a single event cannot measure an exponent, is §4.7 of docs/REFACTOR_PLAN.md.

specmod.spreading.SPREADING_MODELS: dict[str, Any] = {'piecewise': <class 'specmod.spreading.Piecewise'>, 'power_law': <class 'specmod.spreading.PowerLaw'>, 'tabulated': <class 'specmod.spreading.Tabulated'>}#: Registered models, by the name configuration refers to them by.

class specmod.spreading.Piecewise(segments, reference_km=1.0, name='piecewise')[source]#

Bases: object

A contiguous piecewise power law, the shape regional models are published in.

segments is ((exponent, upper_km), ...) in increasing distance. Each segment starts where the previous one ended, so the decay accumulates as a product and the curve is continuous across every hinge — which is what the published tables mean.

Beyond the last boundary the final exponent continues rather than raising, so one far station cannot drop an event. The segment table records where the fitted evidence stopped.

Parameters:

segments (tuple[tuple[float, float], ...])
reference_km (float)
name (str)

class specmod.spreading.PowerLaw(exponent=1.0, reference_km=1.0, name='power_law')[source]#

Bases: object

(R_0 / R) ** exponent — a single power law.

exponent=1 is body-wave amplitude decay in a homogeneous whole space and is the default; 0.5 is the surface-wave value, for distance ranges where an Lg phase dominates the window. Note that 1/R**2 is how energy decays, and a moment expression corrects an amplitude.

Parameters:

exponent (float)
reference_km (float)
name (str)

class specmod.spreading.SpreadingModel(*args, **kwargs)[source]#

Bases: Protocol

Amplitude decay from the reference distance to the site.

name and reference_km are read-only so frozen dataclasses satisfy the protocol.

class specmod.spreading.Tabulated(distances_km, values, reference_km=1.0, name='tabulated')[source]#

Bases: object

A spreading curve supplied as data, interpolated in log-log space.

For spreading that has no functional form — a non-parametric G(R) inversion, for instance. Interpolation is linear in log10 of both axes, the space such curves are inverted in.

Parameters:

distances_km (tuple[float, ...])
values (tuple[float, ...])
reference_km (float)
name (str)

specmod.spreading.get_spreading_model(name, **kwargs)[source]#

Build a registered spreading model by name.

Parameters:

name (str)
kwargs (Any)

Return type:

SpreadingModel

HOLT_2019_UTAH is the piecewise model fitted in Holt (2019), pre-built: Piecewise(segments=((0.90, 43.0), (2.57, 76.0), (0.44, 136.0), (1.54, 400.0))). It is excluded above because autodoc cannot format a signature for a callable dataclass instance documented as module data — inspect.signature handles it fine, autodoc’s own formatter raises. Excluding one constant is cheaper than working around that.

Output#

Writing an event’s spectra to disk and reading them back.

Replaces spectral.Spectra.write_spectra/read_spectra, which pickled the container. Pickle stores the import path of every class it holds, so a stored result stops loading the moment a class is renamed or moved — which is exactly what happened to the shipped Tutorial/Spectra/*.spec: unreadable since the Spectral.py → spectral.py rename, years before the classes were deleted. A format that breaks when you refactor is a cache, not a format.

HDF5, per REFACTOR_PLAN §4.6. The four rules there are lessons from how pickle failed rather than general good practice, and each is visible in the layout below:

Never store class identity.: Plain arrays and a documented layout. Nothing here names a Python type, so nothing here can be broken by renaming one.
Every file carries a format version.: FORMAT_VERSION, on the root group. A reader checks it and fails naming both versions. The absence of this is the whole problem with what came before.
Self-describing units.: motion, kind, duration and sampling_rate are stored attributes, not conventions the reader has to know. This is the typing of Spectrum expressed on disk, and it is what stops a file being silently misread as displacement.
One file per event, one group per channel.: Which matches how the science is done and how it is re-examined, and sidesteps HDF5’s single-writer limitation if the workflow is ever parallelised across events.

Nothing here can pickle. h5py has no mechanism to; trace metadata is stored as a JSON string attribute rather than as an object array, which is how numpy would otherwise smuggle pickling back in.

The tables — fit results — go to Parquet instead, in specmod.tables. Arrays and tables are used in genuinely different ways (random access into one event, versus a columnar scan over every event ever fitted), and one format for both would be worse at each.

specmod.io.FORMAT_VERSION = 1#: Bumped when the layout changes in a way an older reader cannot handle.

specmod.io.load(path)[source]#

Read back what save() wrote.

The arrays come back read-only, as they went in: a spectrum loaded from disk gives the same immutability guarantee as one just computed, which is what makes reload-and-refit safe to do in a loop.

Parameters:: path (str | Path)
Return type:: SpectrumSet

specmod.io.save(path, spectra)[source]#

Write an event’s spectra to path.

One group per trace id, keyed by the id itself — HDF5 group names may contain dots, so LV.L001..HHE needs no mangling and the file browses with the same names the pipeline uses.

The parent directory is created if absent. The legacy version raised FileNotFoundError from inside open, naming the file rather than the directory that did not exist.

Parameters:

path (str | Path)
spectra (SpectrumSet)

Return type:

Path

Fit results as a table, on disk.

The other half of REFACTOR_PLAN §4.6. Spectra are arrays asked about one event at a time and live in HDF5 (specmod.io); fit results are a columnar scan over every event ever fitted — “give me f_c and Omega for all 635 events and regress them” — and live here. The published Magna run produced 11,226 rows of exactly that, and a multi-event catalogue is larger.

Parquet as the primary, CSV as an export. This is not fashion. CSV loses dtypes, so a column of floats comes back as floats only if pandas guesses right and as strings if one cell says None; it round-trips every float through decimal text; and it has to be read in full to read any of it. Parquet is typed, compressed, and queryable with DuckDB or polars without loading the file. CSV stays because journal supplements want it, and because a human with a text editor is a legitimate reader.

The format follows from the suffix, so the choice is visible at the call site rather than buried in a keyword.

specmod.tables.read_table(path)[source]#

Read a fit table back, choosing the format from path’s suffix.

Parameters:: path (str | Path)
Return type:: DataFrame

specmod.tables.write_table(path, table, *, meta=None)[source]#

Write a fit table, choosing the format from path’s suffix.

meta is stored in the file’s own metadata for Parquet — where it survives as key-value pairs a reader can get at without parsing the data — and dropped for CSV, which has nowhere to put it. That asymmetry is stated rather than hidden: a CSV export is lossy about provenance, and pretending otherwise is how a run stops being reproducible.

The parent directory is created if absent.

Parameters:

path (str | Path)
table (pd.DataFrame)
meta (Mapping[str, Any] | None)

Return type:

Path

Looking at a spectrum, a pair, or a whole event.

Replaces spectral.SNP.quick_vis and spectral.Spectra.quick_vis, which were methods on the mutable containers and went with them. A spectral package where you cannot look at a spectrum is not usable, so this is not optional furniture.

Functions over methods. The legacy version was a method, so plotting a spectrum required owning the container it lived in — which is why the fitter grew its own near-duplicate rather than reusing it. These take a SpectrumPair and draw it; anything that has one can call them, and a caller supplying their own Axes gets full control of the figure.

Nothing here mutates its argument, and nothing calls plt.show(). Returning the axes is what lets a caller compose these into a larger figure, annotate them, or save without a window ever opening — which is also what makes them usable from a script and from a notebook without behaving differently.

specmod.plotting.plot_pair(pair, ax=None, *, id='', fit=None, show_binned=False)[source]#

Draw one signal-and-noise pair, with the band it selected.

Parameters:

pair (SpectrumPair) – What to draw.
ax (Axes | None) – Draw here; a new figure is made if omitted.
id (str) – Station label. A frozen pair does not carry one, so the caller that knows the key supplies it — pair.signal.meta["id"] is used when it is there and this is not.
fit (Any) –
A FitSpectrum whose model to overlay, if it has been fitted — or a mapping of label to fit, to draw several at once. Passed in rather than read off the spectrum: the containers are frozen precisely so a result cannot write itself back into its own input.

Several matters more than it sounds. Fitting a source model is not a unique inversion, and two minimisers can reach the same goodness of fit at corner frequencies differing by tens of percent — which is a factor of several in stress drop, since it scales as fc**3. Drawing them together is how that stops being invisible.
show_binned (bool) – Also draw the log-binned spectra the signal-to-noise ratio is actually computed on. Off by default because it doubles the lines, on when the question is why a band came out where it did.

Return type:

Axes

specmod.plotting.plot_set(spectra, *, fits=None, columns=None, passing_only=False, **kwargs)[source]#

Draw every pair in an event on one grid.

fits may be a FitSpectra, or any mapping from trace id to a fit; each pair gets its own model overlaid where one exists.

columns defaults to viz.plot_columns in the configuration, which is the one place that number lives now — it used to be defined in both the SPECTRAL and FITTING dicts, where the two copies could disagree.

Parameters:

spectra (SpectrumSet)
fits (Any)
columns (int | None)
passing_only (bool)
kwargs (Any)

Return type:

Figure

Configuration#

Configuration: semantic sections, layered overrides, recorded provenance.

See docs/REFACTOR_PLAN.md §4.7.

Defaults reproduce the behaviour shipped before the refactor. A study pins its own values in a committed TOML file; personal experimentation goes in specmod.local.toml, which is gitignored, and is promoted deliberately with specmod config freeze.

class specmod.config.AcquireConfig(client='IRIS', event_id=None, origin_time=None, latitude=None, longitude=None, depth_km=None, magnitude=None, networks=('*',), stations=('*',), locations=('*',), channels=('HH?', 'BH?', 'EN?'), max_radius_km=400.0, seconds_before=60.0, seconds_after=300.0, remove_response=False)[source]#

Bases: object

Waveform acquisition. Consumed by specmod.acquire.

Parameters:

client (str)
event_id (str | None)
origin_time (str | None)
latitude (float | None)
longitude (float | None)
depth_km (float | None)
magnitude (float | None)
networks (tuple[str, ...])
stations (tuple[str, ...])
locations (tuple[str, ...])
channels (tuple[str, ...])
max_radius_km (float)
seconds_before (float)
seconds_after (float)
remove_response (bool)

client: str#: FDSN data centre. Different centres serve different holdings for the same event, so this is part of the provenance record, not a detail.

origin_time: str | None#

explicit origin.

Type:: Fallback when event_id is not given

seconds_before: float#: Seconds either side of origin to request.

remove_response: bool#: Store raw counts plus the response rather than a deconvolved trace, so response removal stays under test and no ObsPy version is baked in.

class specmod.config.Config(acquire=<factory>, windows=<factory>, geometry=<factory>, transform=<factory>, smoothing=<factory>, snr=<factory>, model=<factory>, fitting=<factory>, viz=<factory>)[source]#

Bases: object

The whole resolved configuration.

Parameters:

acquire (AcquireConfig)
windows (WindowsConfig)
geometry (GeometryConfig)
transform (TransformConfig)
smoothing (SmoothingConfig)
snr (SnrConfig)
model (ModelConfig)
fitting (FittingConfig)
viz (VizConfig)

to_dict()[source]#

Return a plain nested dict, suitable for TOML or JSON.

Return type:: dict[str, Any]

classmethod from_dict(data)[source]#

Build a Config from a nested mapping, validating section names.

Parameters:: data (dict[str, Any])
Return type:: Self

class specmod.config.FittingConfig(method='powell', fit_bins=False, weight_method='none', initial_t_star=0.01, initial_alpha=1e-05, t_star_min=0.0001, corner_frequency_min=0.0, event_parameter='fc', event_weighting='inverse_distance', include=(), exclude=(), require_pass=True)[source]#

Bases: object

Minimisation.

Parameters:

method (str)
fit_bins (bool)
weight_method (Literal['none', 'log'])
initial_t_star (float)
initial_alpha (float)
t_star_min (float)
corner_frequency_min (float)
event_parameter (str)
event_weighting (str)
include (tuple[str, ...])
exclude (tuple[str, ...])
require_pass (bool)

initial_t_star: float#: Initial guesses that were hardcoded in ModelGuess.

t_star_min: float#

Lower bounds on the fitted parameters that have one physically.

A negative t* says the wave gained energy travelling; a corner frequency at or below zero is not a poor measurement but a meaningless one. Neither is prevented by the misfit surface, and lmfit will return either if the surface leans that way — with the shipped multitaper default it returned fc = -4.45 Hz on one PNR station.

Zero rather than a small positive number, deliberately: a parameter that lands on its bound is flagged by pass_fitting, so the fit is rejected rather than reported as a corner frequency of nothing.

event_parameter: str#

The two-stage event fit; see specmod.staged.

One spectrum cannot separate the source corner from the path attenuation — they trade off on the falling limb — so the corner is determined by the ensemble and then held fixed while each station refits the rest. event_parameter is what the ensemble decides. "fc" because that is the term belonging to the source; "ts" is the meaningful alternative for a study with an independent handle on Q.

event_weighting: str#: How stations are weighted into the event value. The published choice is inverse hypocentral distance: the nearer station has less path, so less of its falloff can be attenuation. See specmod.staged.WEIGHT_MODELS.

include: tuple[str, ...]#

Which channels contribute to the event value, as shell globs matched against the trace id and each of its SEED components — so "AQ07" means the station, "HHE" means the component, "UR" means the network. Empty include means “everything not excluded”.

These exist to be edited after looking at a first pass. Quality control is a judgement — a clipped record, a bad response, a pick on the wrong phase — and a station that is confidently wrong moves the event value for every other station. Putting the decision in the study file is what makes it part of the record rather than something done in a notebook and forgotten.

require_pass: bool#: Drop a station whose stage-1 fit ended with a parameter pinned against one of its bounds. The value reported there is the bound rather than a measurement, so averaging it in is averaging in a constant.

class specmod.config.GeometryConfig(distance_measure='repi')[source]#

Bases: object

Source-to-site geometry.

Its own section because more than one stage needs it. Distance feeds the ensemble weighting of the two-stage fit (specmod.staged) and the geometric spreading a moment calculation corrects for, and a setting two consumers each keep their own copy of is how the two come to disagree.

It lived in [windows] until there was a second reader, which was the wrong home even then: cutting a window does not depend on how distance is measured.

Parameters:: distance_measure (str)

distance_measure: str#

Which distance, resolved through specmod.distance.DISTANCE_MEASURES.

repi is the default and is the honest one wherever sensor depths are not known. rhyp is built from the source depth and the station elevation, so it assumes every sensor sits at the surface — for a borehole deployment that is wrong by the burial depth, and nothing in the metadata says so.

The choice is not a detail at short range: on the PNR data the nearest station is 1.02 km epicentral against 2.30 km hypocentral, a factor of 2.24, while the farthest agree to 1.004 — so anything weighted by inverse distance is most sensitive to it exactly where it matters most.

rrup and rjb are registered and raise: both need a rupture surface, and for a point source they degenerate to rhyp and repi.

class specmod.config.ModelConfig(source='brune', motion='velocity', frequency_dependent_attenuation=False)[source]#

Bases: object

Source and attenuation model.

Parameters:

source (Literal['brune', 'boatwright'])
motion (Literal['displacement', 'velocity', 'acceleration'])
frequency_dependent_attenuation (bool)

motion: Literal['displacement', 'velocity', 'acceleration']#: The motion the model is expressed in. Once Spectrum carries its own motion (§4.2) this is a default rather than a global that must be kept in sync by hand.

class specmod.config.Provenance(specmod_version, config, config_hash, created_at, sources)[source]#

Bases: object

The record attached to every output.

Parameters:

specmod_version (str)
config (dict[str, Any])
config_hash (str)
created_at (str)
sources (dict[str, str])

class specmod.config.ResolvedConfig(config, sources)[source]#

Bases: object

A Config plus the layer each value came from.

The provenance is what makes specmod config show able to answer “why did this run differ”, which is otherwise guesswork once local overrides and environment variables are in play.

Parameters:

config (Config)
sources (dict[str, str])

sources#: Maps "section.key" to the name of the layer that set it.

explain()[source]#

Render the resolved config with the origin of every value.

Return type:: str

class specmod.config.SmoothingConfig(method='log_bins', f_min=0.001, f_max=200.0, n_bins=151, konno_ohmachi_bandwidth=40.0)[source]#

Bases: object

Spectral smoothing and log-space binning.

Parameters:

method (Literal['log_bins', 'konno_ohmachi', 'none'])
f_min (float | None)
f_max (float | None)
n_bins (int)
konno_ohmachi_bandwidth (float)

f_min: float | None#

fmin from 1/T, fmax from Nyquist. The old code hardcoded 0.001-200 Hz regardless.

Type:: Log bin edges. None derives them from the record

konno_ohmachi_bandwidth: float#: Konno-Ohmachi bandwidth b. Smaller smooths harder.

class specmod.config.SnrConfig(tolerance=3.0, min_points=10, assert_bandwidths=False, bands=((2.0, 4.0), (4.0, 6.0), (6.0, 8.0)), scale_parseval=True, interpolate_noise=True, bandwidth_method='peak', resolution_floor=True, rotate_noise=True, rotation_method='boost', rotation_increment=0.05, rotation_space=(0.001, 1.001))[source]#

Bases: object

Signal-to-noise assessment and usable bandwidth selection.

Parameters:

tolerance (float)
min_points (int)
assert_bandwidths (bool)
bands (tuple[tuple[float, float], ...])
scale_parseval (bool)
interpolate_noise (bool)
bandwidth_method (Literal['widest', 'peak'])
resolution_floor (bool)
rotate_noise (bool)
rotation_method (Literal['rotate', 'boost'])
rotation_increment (float)
rotation_space (tuple[float, float])

assert_bandwidths: bool#: Require SNR above tolerance in every band. The published Magna run used this as its selection criterion; it ships disabled.

scale_parseval: bool#: Scale noise amplitude by sqrt(len(signal)/len(noise)) when the noise window is shorter than the signal window.

bandwidth_method: Literal['widest', 'peak']#: Names come from specmod.core.bandwidth.BANDWIDTH_SELECTORS. This said "integral" while the registry said "widest", from the period when the selector was still a percentile of a sign integral — a config value that named nothing the code would accept.

resolution_floor: bool#: Impose a low-frequency floor from the window length (~1/T), or the cone of influence when the spectrum came from a CWT. Nothing enforced this before, so a short window could report bandwidth it could not resolve.

class specmod.config.TransformConfig(estimator='multitaper', time_bandwidth=3.0, n_tapers=5, adaptive=True, normalize_to_variance=False, taper='tukey', taper_alpha=0.05, n_fft=None, welch_segment_length=None, wavelet='morlet', omega0=6.0, dj=0.125, mask_coi=True, drop_dc=True)[source]#

Bases: object

Time-to-frequency conversion. Consumed by specmod.transforms.

Parameters:

estimator (Literal['multitaper', 'fft', 'welch', 'cwt', 'prieto', 'quadratic', 'mtspec'])
time_bandwidth (float)
n_tapers (int)
adaptive (bool)
normalize_to_variance (bool)
taper (Literal['hann', 'tukey', 'boxcar'])
taper_alpha (float)
n_fft (int | str | None)
welch_segment_length (int | None)
wavelet (Literal['morlet'])
omega0 (float)
dj (float)
mask_coi (bool)
drop_dc (bool)

time_bandwidth: float#: Multitaper. time_bandwidth was previously the literal 3 passed positionally to mtspec, with no way to configure it.

adaptive: bool#

leakage suppression is the point of multitaper, and flat weighting leaves the high-frequency floor ~287x high under a strong low-frequency peak. See specmod.transforms.multitaper.

Type:: On by default

normalize_to_variance: bool#: Rescale the spectrum to integrate to the record variance, as mtspec and Prieto’s multitaper do. Needed to reproduce pre-refactor results; off by default because it makes the Parseval check circular.

taper: Literal['hann', 'tukey', 'boxcar']#: FFT / Welch.

n_fft: int | str | None#: None for no padding, an integer, or “fast”/”pow2”. Padding is a pure interpolation here – the normalisation is keyed off duration, not len(freq), which is what the old psd_to_amp got wrong. Use “fast” to avoid the slow FFT path on prime-length cut windows.

wavelet: Literal['morlet']#: Continuous wavelet transform.

dj: float#

voices per octave.

Type:: Scale resolution

drop_dc: bool#: Drop the DC bin. The old code did this unconditionally and before using len(freq) for normalisation, biasing amplitudes slightly.

class specmod.config.VizConfig(plot_columns=3)[source]#

Bases: object

Plotting.

PLOT_COLUMNS was previously defined in both the SPECTRAL and FITTING dicts, and the two copies could disagree. One home makes that impossible.

Parameters:: plot_columns (int)

class specmod.config.WindowsConfig(p_velocity=5.9, s_velocity=2.9, emergency_ratio=1.7, s_start_ratio=0.8, s_length=20.0, s_length_mode='absolute_time', p_before=0.0, p_length=0.8, p_length_mode='relative_time', refine=True, refine_percentiles=(1.0, 99.0), noise_shift=0.2, noise_length=1.0, pad_seconds=0.0, pad_value=0.0, station_shifts=<factory>)[source]#

Bases: object

Phase arrivals and signal/noise window construction.

Parameters:

p_velocity (float)
s_velocity (float)
emergency_ratio (float)
s_start_ratio (float)
s_length (float)
s_length_mode (Literal['absolute_time', 'relative_ps'])
p_before (float)
p_length (float)
p_length_mode (Literal['absolute_time', 'relative_time'])
refine (bool)
refine_percentiles (tuple[float, float])
noise_shift (float)
noise_length (float)
pad_seconds (float)
pad_value (float)
station_shifts (dict[str, float])

p_velocity: float#: Group velocities (km/s) for theoretical arrivals. The published Magna run used s=3.4; 2.9 is the shipped default and is kept as such.

emergency_ratio: float#

s_time = p_time + emergency_ratio * (p - o).

Type:: Used when an S pick is missing

s_start_ratio: float#

opens at s_start_ratio of the P-S time, runs s_length.

Type:: S window

p_before: float#: P window.

refine: bool#: Refine windows to percentiles of the cumulative squared-amplitude integral. This is step 5 of the published Magna workflow.

noise_shift: float#: Noise window ends this many seconds before the P arrival. The published run used 0.5; 0.2 is the shipped default.

specmod.config.config_hash(config, *, length=12)[source]#

Short, stable digest of a configuration.

Comparing two runs starts here: same hash means same settings, so any difference is in the data or the code, not the configuration.

Parameters:

config (Config)
length (int)

Return type:

str

specmod.config.load_config(start=None, *, project_file=None, use_local=True, use_env=True, **overrides)[source]#

Resolve configuration through all layers.

Parameters:

start (Path | str | None) – Directory to search for config files. Defaults to the current directory. The search does not walk upwards — an implicit parent search makes it unclear which file a run actually used.
project_file (Path | str | None) – Explicit path to a committed config, bypassing the search. This is how a study config (studies/magna_2020_paper.toml) is applied, and how tests pin an explicit configuration rather than inheriting defaults.
use_local (bool) – Disable the local-file and environment layers. Tests set both to False so a developer’s machine cannot influence a result.
use_env (bool) – Disable the local-file and environment layers. Tests set both to False so a developer’s machine cannot influence a result.
**overrides (dict[str, Any]) – Section-keyed dicts, e.g. snr={"tolerance": 4}. Highest precedence.

Return type:

ResolvedConfig

Reading picks and catalogues, and looking at waveforms.

What is left after the migration. The rotation machinery that used to live here — find_rotation_angle, find_rotation_angle_v2, rotate, rotate_noise_full, get_centroid_freq and non_lin_boost_noise_func — is now specmod.core.noise, where the two methods are registered models rather than an integer flag. The SAC-discovery scaffolding that used to sit between them (DataSet and its channel-ranking helpers, cps, path_to_utc) had no callers anywhere and was hardcoded to another study’s station codes; it is replaced by the acquisition layer in REFACTOR_PLAN §5.2.

specmod.utils.logger = <Logger specmod.utils (WARNING)>#

diagnostics go here, and nothing in this package configures logging on a caller’s behalf.

Type:: See the note in specmod.fitting.event

specmod.utils.read_pyrocko_picks(path)[source]#

Snuffler marker file to a PickSet.

Thin wrapper over specmod.picks.SnufflerReader, which is where the format lives.

Parameters:: path (str | PathLike[str])
Return type:: PickSet

specmod.utils.read_pyrocko(path)[source]#

Snuffler marker file to {"NET.STA.LOC": {"P": UTCDateTime, ...}}.

Keyed per sensor, not per channel. An arrival is an observation of one sensor, and the component it was picked on is incidental — on the shipped PNR file every station has P on HHZ and S on HHN, so keying by full SEED id would leave the horizontals with no S at all.

The location code stays in the key, because that is what distinguishes two sensors at one site. A surface and a borehole instrument share a station name and see genuinely different arrivals; collapsing them would silently give one the other’s picks.

A flat mapping cannot express a partial identity or a second event; prefer read_pyrocko_picks() with specmod.picks.resolve().

Parameters:: path (str | PathLike[str])
Return type:: dict[str, dict[str, UTCDateTime]]

specmod.utils.read_quakeml_picks(source)[source]#

QuakeML picks to {"NET.STA.LOC": {"P": UTCDateTime, ...}}.

The same shape read_pyrocko() returns, so the two are interchangeable at the call site. source is a path or an obspy.core.event.Catalog.

Phase hints are normalised to P or S on their first letter, which folds Pg/Pn/Pb and their S counterparts together. That matches what this pipeline does with them — it wants a direct arrival, and the branch distinction is not something the windowing uses.

Picks with an evaluation_status of rejected are dropped. Anything else — reviewed, preliminary, unset — is kept, matching how the Snuffler reader treats weights.

A multi-event source raises rather than merging. A flat mapping cannot express a partial identity or a second event; prefer specmod.picks.read() with specmod.picks.resolve().

Parameters:: source (str | PathLike[str] | Catalog)
Return type:: dict[str, dict[str, UTCDateTime]]

specmod.utils.picks_to_quakeml(picks, *, event_id=None)[source]#

The inverse: a pick mapping as an obspy Catalog.

Used to convert a Snuffler marker file to the standard format once, rather than teaching every consumer both. The location code -- is written back as an empty string, which is what StationXML and miniSEED use.

Parameters:

picks (dict[str, dict[str, UTCDateTime]])
event_id (str | None)

Return type:

Catalog

specmod.utils.stream_distance_sort(st, dist_met='repi')[source]#

A copy of st ordered by distance. Not in place.

Returns the stream unsorted, with a warning, when the traces carry no distance — which is why the copy is taken on the way out rather than only on the sorted path.

Parameters:

st (Stream)
dist_met (str)

Return type:

Stream

specmod.utils.cat2kstyle(row)[source]#

Catalogue date and time as dot-separated whole-number fields.

Sub-second precision is dropped, because keith2utc() parses every field with int(). It used to be dropped with a fixed [:-3] slice, which silently assumed exactly two decimal places on the seconds:

"13:09:31.00" worked,
"13:09:31.000" left a trailing separator and made keith2utc raise invalid literal for int(),
"13:09:31" — no decimals at all — lost the seconds, quietly returning 13:09 as though it were 13:09:31.

The last is the one worth fixing: it is a wrong answer rather than an error, and a catalogue written without fractional seconds is not unusual. Splitting on the decimal point handles all three.

Parameters:: row (Any)
Return type:: str