import csv
from functools import cached_property
from typing import Self
import numpy as np
from . import signal
from ..io.record import Record
from ..optimization.fit_eval import goodness_of_fit, relative_error
[docs]
class GroundMotion:
"""
Container for ground motion data and related operations.
This class stores acceleration, velocity, and displacement time series
along with derived intensity measures. All transformation methods return
new instances rather than modifying in place.
Parameters
----------
dt : float
Time step in seconds.
ac : ndarray
Acceleration time series in **cm/s²**.
tag : str, optional
Identifier for the ground motion record (default is None).
Notes
-----
Ground motion instances should be treated as immutable. Direct modification
of ac or dt may lead to inconsistent cached properties.
Examples
--------
Load from file:
>>> gm = GroundMotion.from_file("RSN1.AT2", fmt="NGA", unit="g")
>>> gm.pga
0.45
Create from arrays:
>>> import numpy as np
>>> dt = 0.01
>>> ac = np.random.randn(1000)
>>> gm = GroundMotion.from_array(dt, ac, unit="cm/s2")
>>> gm_trimmed = gm.trim_by_energy((0.05, 0.95))
>>> gm_trimmed.npts
900
"""
def __init__(self, dt, ac, tag=None):
self.dt = float(dt)
self.ac = np.asarray(ac, dtype=np.float64)
self.tag = tag
# Class methods ==================================================================
[docs]
@classmethod
def from_record(cls, record: Record, tag: str = None) -> Self:
"""
Create a GroundMotion instance from a Record object.
Parameters
----------
record : Record
Parsed ground motion record.
tag : str, optional
Identifier for the ground motion record.
Returns
-------
GroundMotion
"""
return cls(record.dt, record.ac, tag=tag)
[docs]
@classmethod
def from_file(
cls,
path: str,
fmt: str,
unit: str = "cm/s2",
skiprows: int = 1,
tag: str = None,
) -> Self:
"""
Load a record from a plain-text file.
Parameters
----------
path : str or Path
Path to the record file.
fmt : str
Format key — ``'NGA'``, ``'ESM'``, ``'COL'``, ``'RAW'``,
``'COR'`` (case-insensitive).
unit : str, optional
Unit of the acceleration values *as stored in the file*.
The loaded array is converted to cm/s² automatically.
Defaults to ``'cm/s2'``.
skiprows : int, optional
Header rows to skip — only relevant for the ``'COL'`` format
(default 1).
tag : str, optional
Record identifier.
Returns
-------
GroundMotion
"""
rec = Record.from_file(path, fmt, unit=unit, skiprows=skiprows)
return cls.from_record(rec, tag=tag)
[docs]
@classmethod
def from_zip(
cls,
zip_path: str,
filename: str,
fmt: str,
unit: str = "cm/s2",
skiprows: int = 1,
tag: str = None,
) -> Self:
"""
Load a record from a file inside a zip archive.
Parameters
----------
zip_path : str or Path
Path to the ``.zip`` archive.
filename : str
Name of the target file inside the archive.
fmt : str
Format key — ``'NGA'``, ``'ESM'``, ``'COL'``, ``'RAW'``,
``'COR'`` (case-insensitive).
unit : str, optional
Unit of the acceleration values *as stored in the file*.
The loaded array is converted to cm/s² automatically.
Defaults to ``'cm/s2'``.
skiprows : int, optional
Header rows to skip — only relevant for the ``'COL'`` format
(default 1).
tag : str, optional
Record identifier.
Returns
-------
GroundMotion
"""
rec = Record.from_zip(zip_path, filename, fmt, unit=unit, skiprows=skiprows)
return cls.from_record(rec, tag=tag)
[docs]
@classmethod
def from_array(
cls, dt: float, ac: np.ndarray, unit: str = "cm/s2", tag: str = None
) -> Self:
"""
Wrap an existing array as a ground motion record.
Parameters
----------
dt : float
Time step in seconds.
ac : array-like
1-D acceleration values expressed in ``unit``.
unit : str, optional
Unit of ``ac``. The array is converted to cm/s² automatically.
Defaults to ``'cm/s2'``.
tag : str, optional
Record identifier.
Returns
-------
GroundMotion
"""
rec = Record.from_array(ac, dt, unit=unit)
return cls.from_record(rec, tag=tag)
# Core derived series ========================================================
[docs]
@cached_property
def npts(self) -> int:
"""
Number of time points in the record.
Returns
-------
int
"""
return int(self.ac.shape[-1])
[docs]
@cached_property
def vel(self) -> np.ndarray:
"""
Velocity time series.
Computed by integrating the acceleration time series.
Returns
-------
ndarray
"""
return signal.integrate(self.dt, self.ac)
[docs]
@cached_property
def disp(self) -> np.ndarray:
"""
Displacement time series.
Computed by integrating the velocity time series.
Returns
-------
ndarray
"""
return signal.integrate(self.dt, self.vel)
# Processing Methods =========================================================
[docs]
def trim_by_index(self, start_index: int, end_index: int) -> Self:
"""
Trim ground motion by index range.
Extracts a subset of the time series between specified indices,
creating a new GroundMotion instance with reduced duration.
Parameters
----------
start_index : int
Starting index (inclusive).
end_index : int
Ending index (exclusive).
Returns
-------
GroundMotion
New instance with trimmed time series.
Raises
------
ValueError
If indices are out of bounds (start_index < 0 or end_index > npts).
See Also
--------
trim_by_slice : Trim using Python slice notation.
trim_by_energy : Trim based on cumulative energy range.
"""
if start_index < 0 or end_index > self.npts:
raise ValueError("start_index and end_index must be within current npts")
return type(self)(
dt=self.dt, ac=self.ac[..., start_index:end_index], tag=self.tag
)
[docs]
def trim_by_slice(self, slicer: slice) -> Self:
"""
Trim ground motion using Python slice object.
Provides flexible slicing similar to NumPy array indexing with support
for negative indices and step values.
Parameters
----------
slicer : slice
Python slice object (e.g., slice(100, 500, 2)).
Returns
-------
GroundMotion
New instance with sliced time series.
Raises
------
TypeError
If slicer is not a slice object.
See Also
--------
trim_by_index : Trim with explicit start/end indices.
"""
if not isinstance(slicer, slice):
raise TypeError("Expected a slice object")
return type(self)(dt=self.dt, ac=self.ac[..., slicer], tag=self.tag)
[docs]
def trim_by_energy(self, energy_range: tuple[float, float]) -> Self:
"""
Trim ground motion to retain specified cumulative energy range.
Identifies time window containing the target energy range (e.g., 5%-95%)
based on cumulative energy of acceleration. Useful for focusing on
significant motion and removing weak pre/post-event portions.
Parameters
----------
energy_range : tuple of float
(start_fraction, end_fraction) where fractions are in [0, 1].
Example: (0.05, 0.95) retains central 90% of energy.
Returns
-------
GroundMotion
New instance trimmed to energy range.
Raises
------
ValueError
If fractions are not in [0, 1] or start >= end.
"""
slicer = signal.slice_energy(self.ce, energy_range)
return type(self)(dt=self.dt, ac=self.ac[..., slicer], tag=self.tag)
[docs]
def trim_by_amplitude(self, threshold: float) -> Self:
"""
Trim ground motion based on acceleration amplitude threshold.
Identifies the time window where acceleration exceeds the specified
threshold, removing weak motion at start and end.
Parameters
----------
threshold : float
Amplitude threshold in same units as acceleration (typically g).
Returns
-------
GroundMotion
New instance trimmed to significant motion window.
"""
slicer = signal.slice_amplitude(self.ac, threshold)
return type(self)(dt=self.dt, ac=self.ac[..., slicer], tag=self.tag)
[docs]
def taper(self, alpha: float = 0.05) -> Self:
"""
Apply Tukey window tapering to ground motion.
Smoothly tapers the beginning and end of the time series to zero using
a Tukey (tapered cosine) window. Reduces spectral leakage in frequency
domain analysis and prevents edge effects in filtering.
Parameters
----------
alpha : float, optional
Taper fraction in [0, 1]. Fraction of window inside cosine tapered region.
- 0: Rectangular window (no tapering)
- 1: Hann window (full taper)
- 0.05: Default, tapers 5% at each end (default is 0.05).
Returns
-------
GroundMotion
New instance with tapered acceleration.
"""
new_ac = signal.taper(self.ac, alpha)
return type(self)(dt=self.dt, ac=new_ac, tag=self.tag)
[docs]
def butterworth_filter(
self, bandpass_freqs: tuple[float, float], order: int = 4
) -> Self:
"""
Apply Butterworth bandpass filter using second-order sections (SOS).
Zero-phase Butterworth filter for frequency content selection without
introducing phase distortion. Uses SOS format for improved numerical
stability compared to transfer function representation.
Parameters
----------
bandpass_freqs : tuple of float
(low_freq, high_freq) in Hz. Defines passband range.
order : int, optional
Filter order controlling steepness of rolloff (default is 4).
Higher orders give sharper cutoffs but may introduce instability.
Returns
-------
GroundMotion
New instance with filtered acceleration.
"""
new_ac = signal.butterworth_filter(self.dt, self.ac, *bandpass_freqs, order)
return type(self)(dt=self.dt, ac=new_ac, tag=self.tag)
[docs]
def baseline_correction(self, degree: int = 1) -> Self:
"""
Apply polynomial baseline correction.
Removes long-period drift by fitting and subtracting a polynomial trend
from acceleration. Common preprocessing step for integrating to velocity
and displacement.
Parameters
----------
degree : int, optional
Polynomial degree for trend fitting (default is 1).
- 0: Remove mean (DC offset)
- 1: Remove linear trend
- 2: Remove quadratic trend
Returns
-------
GroundMotion
New instance with corrected acceleration.
"""
new_ac = signal.baseline_correction(self.ac, degree)
return type(self)(dt=self.dt, ac=new_ac, tag=self.tag)
[docs]
def resample(self, dt: float) -> Self:
"""
Resample to new time step using Fourier method.
Changes time step by resampling in frequency domain, preserving
frequency content up to new Nyquist frequency.
Parameters
----------
dt : float
New time step in seconds.
Returns
-------
GroundMotion
New instance with resampled time step.
"""
_, dt_new, ac_new = signal.resample(self.dt, dt, self.ac)
return type(self)(dt=dt_new, ac=ac_new, tag=self.tag)
[docs]
def response_spectra(self, periods: np.ndarray, damping: float = 0.05):
"""
Calculate response spectra for given periods and damping.
Computes spectral displacement (Sd), velocity (Sv), and acceleration (Sa)
for elastic single-degree-of-freedom oscillators.
Parameters
----------
periods : ndarray
Array of natural periods in seconds.
damping : float, optional
Damping ratio (default is 0.05 for 5% damping).
Returns
-------
sd : ndarray
Spectral displacement in cm.
sv : ndarray
Spectral velocity in cm/s.
sa : ndarray
Spectral acceleration in g.
"""
return signal.response_spectra(self.dt, self.ac, period=periods, zeta=damping)
[docs]
def compute_intensity_measures(
self, ims: list[str], periods: np.ndarray = None
) -> dict:
"""
Compute selected intensity measures.
Batch computation of multiple IMs with optimized calculation of
spectral quantities.
Parameters
----------
ims : list of str
IM names to compute (e.g., ['pga', 'sa', 'cav']).
periods : ndarray, optional
Periods for spectral IMs (sa, sv, sd). Required if any spectral IM requested.
Returns
-------
dict
Dictionary with IM names as keys. Spectral IMs have keys like 'sa_0.200'.
Values are floats (single record) or arrays (multiple records).
"""
periods = np.asarray(periods) if periods is not None else None
results = {}
# Pre-compute spectra if requested (optimization)
spectral_ims = [im for im in ims if im.lower() in ("sa", "sv", "sd")]
spectral_data = {}
if spectral_ims:
if periods is None:
raise ValueError(
"Periods must be provided to compute spectral quantities (sa, sv, sd)."
)
# Compute once for all spectral types
sd, sv, sa = self.response_spectra(periods)
spectral_data["sd"] = sd
spectral_data["sv"] = sv
spectral_data["sa"] = sa
# Iterate and collect data
for im in ims:
im_l = im.lower()
# Case A: Spectral IMs
if im_l in spectral_data:
data_matrix = spectral_data[im_l]
for idx, period in enumerate(periods):
key = f"{im_l}_{period:.3f}"
results[key] = data_matrix[..., idx]
# Case B: Fourier Amplitude Spectra
elif im_l == "fas":
for idx, freq in enumerate(self.freq):
key = f"fas_{freq:.3f}"
results[key] = self.fas[..., idx]
# Case C: Scalar IMs
else:
attr = getattr(self, im_l)
results[im_l] = attr
return results
[docs]
def to_csv(self, filename: str, ims: list[str], periods: np.ndarray = None):
"""
Export intensity measures to CSV file.
Writes computed IMs to comma-separated values format suitable for
further analysis or database storage.
Parameters
----------
filename : str
Output file path.
ims : list of str
IM names to export.
periods : ndarray, optional
Periods for spectral IMs.
"""
data = self.compute_intensity_measures(ims, periods)
if not data:
return
fieldnames = list(data.keys())
columns = list(data.values())
# Determine if we have a single row (scalars) or multiple rows (arrays)
if np.isscalar(columns[0]):
rows = [columns]
else:
# zip(*) perfectly transposes columns into rows dynamically
rows = zip(*columns)
with open(filename, "w", newline="", encoding="utf-8") as f:
writer = csv.writer(f)
writer.writerow(fieldnames)
writer.writerows(rows)
[docs]
def compare(
self,
other: Self,
ims: list[str],
periods: np.ndarray = None,
method: str = "gof",
) -> dict:
"""
Compare with another ground motion using goodness-of-fit metrics.
Quantifies similarity between this ground motion and a target/model
using specified intensity measures.
Parameters
----------
other : GroundMotion
Target ground motion for comparison.
ims : list of str
IM names to compare.
periods : ndarray, optional
Periods for spectral IMs.
method : str, optional
Comparison metric: 'gof' (goodness of fit) or 're' (relative error).
Default is 'gof'.
Returns
-------
dict
Comparison scores for each IM. Lower is better for 'gof',
closer to 0 is better for 're'.
"""
criterion_map = {"gof": goodness_of_fit, "re": relative_error}
if method.lower() not in criterion_map:
raise ValueError(
f"Unknown method: {method}. Supported: {list(criterion_map.keys())}"
)
func = criterion_map[method.lower()]
my_data = self.compute_intensity_measures(ims, periods)
other_data = other.compute_intensity_measures(ims, periods)
scores = {}
for key in my_data:
if key in other_data:
scores[key] = func(my_data[key], other_data[key])
return scores
# Properties ========================================================
@property
def vsi(self):
"""
Velocity spectrum intensity (0.1-2.5s range).
"""
return self.spectrum_intensity[1]
@property
def asi(self):
"""
Acceleration spectrum intensity (0.1-2.5s range).
"""
return self.spectrum_intensity[2]
@property
def dsi(self):
"""
Displacement spectrum intensity (0.1-2.5s range).
"""
return self.spectrum_intensity[0]
[docs]
@cached_property
def t(self):
"""Time array corresponding to recorded points."""
return signal.time(self.npts, self.dt)
[docs]
@cached_property
def freq(self):
"""Frequency array for Fourier transform."""
return signal.frequency(self.npts, self.dt)
[docs]
@cached_property
def fas(self):
"""Fourier amplitude spectrum of acceleration."""
return signal.fas(self.dt, self.ac)
[docs]
@cached_property
def fas_vel(self):
"""Fourier amplitude spectrum of velocity."""
return signal.fas(self.dt, self.vel)
[docs]
@cached_property
def fas_disp(self):
"""Fourier amplitude spectrum of displacement."""
return signal.fas(self.dt, self.disp)
[docs]
@cached_property
def fps(self):
"""Fourier phase spectrum of acceleration (unwrapped)."""
return signal.fps(self.ac)
[docs]
@cached_property
def ce(self):
"""Cumulative energy of acceleration time series."""
return signal.ce(self.dt, self.ac)
[docs]
@cached_property
def pga(self):
"""Peak ground acceleration."""
return signal.peak_abs_value(self.ac)
[docs]
@cached_property
def pgv(self):
"""Peak ground velocity."""
return signal.peak_abs_value(self.vel)
[docs]
@cached_property
def pgd(self):
"""Peak ground displacement."""
return signal.peak_abs_value(self.disp)
[docs]
@cached_property
def cav(self):
"""Cumulative absolute velocity."""
return signal.cav(self.dt, self.ac)
[docs]
@cached_property
def spectrum_intensity(self):
"""Spectrum intensities (Sd, Sv, Sa) over 0.1-2.5s period range."""
vsi_tp = np.arange(0.1, 2.5, 0.05)
sd, sv, sa = signal.response_spectra(self.dt, self.ac, period=vsi_tp, zeta=0.05)
dsi = np.sum(sd, axis=-1) * 0.05
vsi = np.sum(sv, axis=-1) * 0.05
asi = np.sum(sa, axis=-1) * 0.05
return dsi, vsi, asi
[docs]
@cached_property
def zc_ac(self):
"""Zero-crossing of acceleration."""
return signal.zc(self.ac)
[docs]
@cached_property
def zc_vel(self):
"""Zero-crossing of velocity."""
return signal.zc(self.vel)
[docs]
@cached_property
def zc_disp(self):
"""Zero-crossing of displacement."""
return signal.zc(self.disp)
[docs]
@cached_property
def pmnm_ac(self):
"""Positive-minima to negative-maxima of acceleration."""
return signal.pmnm(self.ac)
[docs]
@cached_property
def pmnm_vel(self):
"""Positive-minima to negative-maxima ratio of velocity."""
return signal.pmnm(self.vel)
[docs]
@cached_property
def pmnm_disp(self):
"""Positive-minima to negative-maxima ratio of displacement."""
return signal.pmnm(self.disp)
[docs]
@cached_property
def le_ac(self):
"""Mean local extrema of acceleration."""
return signal.le(self.ac)
[docs]
@cached_property
def le_vel(self):
"""Mean local extrema of velocity."""
return signal.le(self.vel)
[docs]
@cached_property
def le_disp(self):
"""Mean local extrema of displacement."""
return signal.le(self.disp)
def __repr__(self) -> str:
tag_str = f", tag={self.tag!r}" if self.tag else ""
return f"GroundMotion(npts={self.npts}, dt={self.dt}{tag_str})"
def __len__(self) -> int:
return self.npts