Tutorial

This tutorial uses a simplified ambient noise scenario to show how pysmo fits into real code. It is not a guide to ambient noise processing. Along the way it covers:

• Defining a custom seismogram class for a specific use case.
• Writing functions that operate on it.
• Using pysmo types to make those functions reusable.

Custom seismogram class

The scenario involves ambient noise data. It needs to track whether earthquake signals are present, but has no need for event information. A dataclass fits this well:

from dataclasses import dataclass, field

import numpy as np
import pandas as pd


@dataclass  # (1)!
class NoiseSeismogram:
    begin_time: pd.Timestamp  # (2)!
    delta: pd.Timedelta = pd.Timedelta(seconds=0.01)  # (3)!
    data: np.ndarray = field(default_factory=lambda: np.array([]))  # (4)!

    @property
    def end_time(self) -> pd.Timestamp:  # (5)!
        if len(self.data) == 0:
            return self.begin_time
        return self.begin_time + self.delta * (len(self.data) - 1)

    contains_earthquake: bool = False  # (6)!

1. dataclass is a decorator that generates the methods a data-holding class normally needs, based on the attributes declared in the class body: one to create instances, one to show them as readable text, one to compare them for equality. It saves writing them by hand.
2. Instance attributes are defined simply by declaring them in the class body with type annotations. Note the use of pandas.Timestamp here. It is used throughout pysmo as the standard type for time information.
3. Attributes can have default values too.
4. Mutable default values (like lists or dictionaries) need field(default_factory=...), so that each instance gets its own separate copy.
5. A read-only end_time property computes the end time from the start time, number of samples, and sampling interval.
6. Finally, an attribute records whether the seismogram contains earthquake signals.

A real project would have more attributes, but this is enough to demonstrate the pattern. Creating an instance:

$ uv run python -i noise_seismogram.py
>>> begin_time = Timestamp("2023-01-01", tz="UTC")
>>> data = np.random.randn(1000)  # Simulated noise data
>>> noise_seis = NoiseSeismogram(begin_time=begin_time, data=data)
>>>

Functions that operate on the new class

Two functions handle the processing:

• check_for_earthquakes(): checks whether earthquake signals are present.
• detrend(): detrends the seismogram data.

A first version:

import scipy
from noise_seismogram import NoiseSeismogram


def check_for_earthquakes(seismogram: NoiseSeismogram) -> None:
    if seismogram.contains_earthquake is True:
        print("Seismogram contains an earthquake.")
    elif seismogram.contains_earthquake is False:
        print("Seismogram does not contain an earthquake.")
    else:
        print("Seismogram earthquake status is unknown.")


def detrend(seismogram: NoiseSeismogram) -> None:
    seismogram.data = scipy.signal.detrend(seismogram.data)

The type hints are correct and mypy confirms it:

$ uv run mypy functions_v1.py
Success: no issues found in 1 source file

With type checking in place, mypy can also identify unreachable code. Running with --warn-unreachable:

$ uv run mypy --warn-unreachable functions_v1.py
functions_v1.py:11: error: Statement is unreachable  [unreachable]
Found 1 error in 1 file (checked 1 source file)

The else branch is unreachable because contains_earthquake is non-optional: it can only be True or False. Removing it:

import scipy
from noise_seismogram import NoiseSeismogram


def check_for_earthquakes(seismogram: NoiseSeismogram) -> None:
    if seismogram.contains_earthquake is True:
        print("Seismogram contains an earthquake.")
    else:  # (1)!
        print("Seismogram does not contain an earthquake.")


def detrend(seismogram: NoiseSeismogram) -> None:
    seismogram.data = scipy.signal.detrend(seismogram.data)

1. At this point seismogram.contains_earthquake can only be False, so the elif check is no longer needed.

Mypy reports no errors here either:

$ uv run mypy --warn-unreachable functions_v2.py
Success: no issues found in 1 source file

Reusing functions in other contexts

Comparing the two functions, only check_for_earthquakes() relies on contains_earthquake, the one attribute specific to this project. The remaining attributes form a common baseline, suggesting detrend() should work with other seismogram classes too. To test this, consider a second project that stores the season alongside seismogram data:

from dataclasses import dataclass, field
from enum import StrEnum

import numpy as np
import pandas as pd


class Season(StrEnum):  # (1)!
    SPRING = "spring"
    SUMMER = "summer"
    AUTUMN = "autumn"
    WINTER = "winter"


@dataclass
class SeasonSeismogram:
    begin_time: pd.Timestamp
    delta: pd.Timedelta = pd.Timedelta(seconds=0.01)
    data: np.ndarray = field(default_factory=lambda: np.array([]))

    @property
    def end_time(self) -> pd.Timestamp:
        if len(self.data) == 0:
            return self.begin_time
        return self.begin_time + self.delta * (len(self.data) - 1)

    season: Season = Season.SUMMER  # (2)!

1. StrEnum limits the values a string attribute can take.
2. Much like NoiseSeismogram, this class has just one project-specific attribute (season).

class SeismogramEndtimeMixin:
    """Add a computed `end_time` property.

    Mix into any class that provides `begin_time`, `delta`, and `data`.
    """

    __slots__ = ()

    @property
    def end_time(self: Seismogram) -> pd.Timestamp:
        """Seismogram end time."""
        if len(self.data) == 0:
            return self.begin_time
        return self.begin_time + self.delta * (len(self.data) - 1)

@dataclass
class SeasonSeismogram(SeismogramEndtimeMixin):  # (1)!
    begin_time: Timestamp
    delta: Timedelta = Timedelta(seconds=0.01)
    data: np.ndarray = field(default_factory=lambda: np.array([]))
    season: Season = Season.SUMMER

Next comes a script that pairs this new class with the detrend() function from earlier. The season_detrend_v*.py scripts that follow are identical apart from which functions_v*.py they import detrend from.

import numpy as np
import pandas as pd
from functions_v2 import detrend
from season_seismogram import Season, SeasonSeismogram


def main() -> None:
    # Create a sample SeasonSeismogram instance with random data
    begin_time = pd.Timestamp(2023, 1, 1, 0, 0, 0)
    data = np.random.randn(1000)
    season_seismogram = SeasonSeismogram(
        begin_time=begin_time, data=data, season=Season.WINTER
    )

    # Use the season_seismogram with the detrend function
    detrend(season_seismogram)


if __name__ == "__main__":
    main()

This script runs correctly:

$ uv run season_detrend_v1.py && echo "success!"
success!

But mypy flags a type mismatch. detrend() expects a NoiseSeismogram and is being passed a SeasonSeismogram:

$ uv run mypy season_detrend_v1.py
season_detrend_v1.py:16: error: Argument 1 to "detrend" has incompatible type "SeasonSeismogram"; expected "NoiseSeismogram"  [arg-type]

Type annotations prevent using non-existent attributes, but they do not require using all of them. detrend() only touches data, which both classes happen to share. That was luck, not design.

Fixing this means amending the type annotations of the detrend() function:

import scipy
from noise_seismogram import NoiseSeismogram
from season_seismogram import SeasonSeismogram  # (1)!


def check_for_earthquakes(seismogram: NoiseSeismogram) -> None:
    if seismogram.contains_earthquake is True:
        print("Seismogram contains an earthquake.")
    else:
        print("Seismogram does not contain an earthquake.")


def detrend(seismogram: NoiseSeismogram | SeasonSeismogram) -> None:
    seismogram.data = scipy.signal.detrend(seismogram.data)

1. SeasonSeismogram has to be imported before it can be used in the annotations.

With these changes, mypy reports no errors:

$ uv run mypy season_detrend_v2.py
Success: no issues found in 1 source file

Introducing pysmo

Writing a custom class for each project is fine. Updating every shared function whenever a new class appears is not. Each new class means touching function annotations, and a change to any class risks breaking the functions that depend on it. The standard solution is an interface between functions and classes: functions target the interface, and classes conform to it.

Pysmo provides such an interface for seismogram (and other) classes. These interfaces use Python's Protocol class. Below is the actual implementation of pysmo's Seismogram interface:

class Seismogram(Protocol):
    """Protocol class to define the `Seismogram` type.

    Examples:
        A function annotated with `Seismogram` accepts any compatible class.
        This one returns the begin time in ISO format:

        ```python
        >>> from pysmo import Seismogram
        >>> from pysmo.classes import SAC  # SAC implements the Seismogram protocol
        >>>
        >>> def example_function(seis_in: Seismogram) -> str:
        ...     return seis_in.begin_time.isoformat()
        ...
        >>> sac = SAC.from_file("example.sac")
        >>> seismogram = sac.seismogram
        >>> example_function(seismogram)
        '2010-02-27T06:44:06.069538+00:00'
        >>>
        ```
    """

    begin_time: pd.Timestamp
    """Seismogram begin time."""

    data: npt.NDArray[np.floating]
    """Seismogram data."""

    delta: pd.Timedelta
    """Seismogram sampling interval.

    Must be a positive `pd.Timedelta`.
    """

    @property
    def end_time(self) -> pd.Timestamp:
        """Seismogram end time."""
        ...

Strip away the docstrings and this looks much like the common structure of NoiseSeismogram and SeasonSeismogram. The key difference is that end_time is declared but not implemented. Protocol classes provide type information only and cannot be instantiated.

Through structural subtyping, any class with the matching structure is treated as a subtype of the protocol. Instances of NoiseSeismogram and SeasonSeismogram therefore satisfy Seismogram as well.

Annotating detrend() with the Seismogram type rather than listing every class:

import scipy
from noise_seismogram import NoiseSeismogram

from pysmo import Seismogram  # (1)!


def check_for_earthquakes(seismogram: NoiseSeismogram) -> None:
    if seismogram.contains_earthquake is True:
        print("Seismogram contains an earthquake.")
    else:
        print("Seismogram does not contain an earthquake.")


def detrend(seismogram: Seismogram) -> None:  # (2)!
    seismogram.data = scipy.signal.detrend(seismogram.data)

1. Seismogram replaces the import of SeasonSeismogram.
2. Any class that satisfies the Seismogram structure is now accepted, with no further changes to detrend().

With Seismogram in place, mypy accepts the season script unchanged:

$ uv run mypy season_detrend_v3.py
Success: no issues found in 1 source file

Conclusion

This tutorial introduced the core ideas behind pysmo rather than its API:

• Pysmo is not centred on a single seismogram class. Monolithic classes tend to reflect the use cases their authors had in mind, not the ones users actually have.
• Custom seismogram classes fit specific use cases well, but create friction when writing reusable code.
• Pysmo addresses this with interfaces, the pysmo types, that capture what different classes have in common. Functions target the interface, and any conforming class works.
• Pysmo types are intentionally narrow: few attributes, almost no methods.

The same principles apply to the processing modules pysmo ships with, so they work just as well outside pysmo as within it.

The Usage chapter goes further: how the types are designed, the conventions they follow, and how to adapt an existing class to them.