pysmo.tools.iccs

Iterative Cross-Correlation and Stack (ICCS).

Warning

This module is being developed alongside a complete rewrite of AIMBAT. Expect major changes until the rewrite is complete.

The ICCS¹ method is an iterative algorithm to rapidly determine the best fitting delay times between an arbitrary number of seismograms. In each iteration, individual seismograms are cross-correlated with the stack of all input seismograms to improve the previous (or initial) phase arrival picks.

The results of ICCS are are similar to those of the the MCCC² method, while also requiring fewer cross-correlations to be computed (each individual seismogram is only cross-correlated with the stack, whereas in MCCC all seismograms are cross-correlated with each other). ICCS is therefore particularly useful to perepare data for a successful MCCC run (e.g. if the initial picks are calculated rather than hand picked).

The iccs module requires seismograms containing extra attributes specific to the ICCS method. Hence it provides a protocol class (ICCSSeismogram) and corresponding Mini class (MiniICCSSeismogram) for that purpose. In addition to the common attributes of a Seismogram in pysmo, the following parameters are required:

Attribute	Description
t0	Initial pick (typically computed). Serves as input only when t1 is not set.
t1	Improved pick. Serves as both input (if not None) and output (always) when running the ICCS algorithm.
select	Determines if a seismogram should be used for the stack. Can be automatically set when autoselectis True during a run.
flip	Determines if the seismogram data should be flipped (i.e. multiplied with -1) before using it in the stack and cross-correlation. Can be automatically set when autoflipis True during a run.

The diagram below shows execution flow, and how the above parameters are used when the ICCS algorithm is executed (see here for parameters and default values)':

flowchart TD
Start(["ICCSSeismograms with initial parameters."])
Stack0["Generate windowed seismograms and create stack from them."]
C["Cross-correlate windowed seismograms with stack to obtain updated picks and normalised correlation coefficients."]
FlipQ{"Is **autoflip**
True?"}
Flip["Toggle **flip** attribute of seismograms with negative correlation coefficients."]
QualQ{"Is **autoselect**
True?"}
Qual1["Set **select** attribute of seismograms to False for poor quality seismograms."]
Stack1["Recompute windowed seismograms and stack with updated parameters."]
H{"Convergence
criteria met?"}
I{"Maximum
iterations
reached?"}
End(["ICCSSeismograms with updated **t1**, **flip**, and **select** parameters."])
Start --> Stack0 --> C --> FlipQ -->|No| QualQ -->|No| Stack1 --> H -->|No| I -->|No| C
FlipQ -->|Yes| Flip --> QualQ
QualQ -->|Yes| Qual1 -->  Stack1
H -->|Yes| End
I -->|Yes| End

---

1. Lou, X., et al. “AIMBAT: A Python/Matplotlib Tool for Measuring Teleseismic Arrival Times.” Seismological Research Letters, vol. 84, no. 1, Jan. 2013, pp. 85–93, https://doi.org/10.1785/0220120033. ↩

2. VanDecar, J. C., and R. S. Crosson. “Determination of Teleseismic Relative Phase Arrival Times Using Multi-Channel Cross-Correlation and Least Squares.” Bulletin of the Seismological Society of America, vol. 80, no. 1, Feb. 1990, pp. 150–69, https://doi.org/10.1785/BSSA0800010150. ↩

Classes:

• ICCS –

  Class to store a list of ICCSSeismograms and run the ICCS algorithm.

• ICCSSeismogram –

  Protocol class to define the ICCSSeismogram type.

• MiniICCSSeismogram –

  Minimal implementation of the ICCSSeismogram type.

Functions:

• plotstack –

  Plot the ICCS stack.

• stack_pick –

  Manually pick t1 in the stack and apply it to all seismograms.

• stack_tw_pick –

  Pick new time window limits in the stack.

• update_pick –

  Update t1 in all seismograms by the same amount.

ICCS

Class to store a list of ICCSSeismograms and run the ICCS algorithm.

The ICCS class serves as a container to store a list of seismograms (typically recordings of the same event at different stations), and to then run the ICCS algorithm when an instance of this class is called. Processing parameters that are common to all seismograms are stored as attributes (e.g. time window limits).

Examples:

We begin with a set of SAC files of the same event, recorded at different stations. All files have a preliminary phase arrival estimate saved in the T0 SAC header, so we can use these files to create instances of the MiniICCSSeismogram class for use with the ICCS class:

>>> from pysmo.classes import SAC
>>> from pysmo.functions import clone_to_mini
>>> from pysmo.tools.iccs import MiniICCSSeismogram
>>> from pathlib import Path
>>>
>>> sacfiles = sorted(Path("iccs-example/").glob("*.bhz"))
>>>
>>> seismograms = []
>>> for sacfile in sacfiles:
...     sac = SAC.from_file(sacfile)
...     update = {"t0": sac.timestamps.t0}
...     iccs_seismogram = clone_to_mini(MiniICCSSeismogram, sac.seismogram, update=update)
...     seismograms.append(iccs_seismogram)
...
>>>

To better illustrate the different modes of running the ICCS algorithm, we modify the data and picks in the seismograms to make them worse than they actually are:

>>> from datetime import timedelta
>>> from copy import deepcopy
>>> import numpy as np
>>>
>>> # change the sign of the data in the first seismogram
>>> seismograms[0].data *= -1
>>>
>>> # move the initial pick 2 seconds earlier in second seismogram
>>> seismograms[1].t0 += timedelta(seconds=-2)
>>>
>>> # move the initial pick 2 seconds later in third seismogram
>>> seismograms[2].t0 += timedelta(seconds=2)
>>>
>>> # create a seismogram with completely random data
>>> iccs_random: MiniICCSSeismogram = deepcopy(seismograms[-1])
>>> np.random.seed(42)  # set this for consistent results during testing
>>> iccs_random.data = np.random.rand(len(iccs_random))
>>> seismograms.append(iccs_random)
>>>

We can now create an instance of the ICCS class and plot the initial stack together with the input seismograms:

>>> from pysmo.tools.iccs import ICCS, plotstack
>>> iccs = ICCS(seismograms)
>>> plotstack(iccs, padded=False)
>>>

The phase emergance is not visible in the stack, and the (absolute) correlation coefficents of the the seismograms are not very high. This shows the initial picks are not very good and/or that the data are of low quality. To run the ICCS algorithm, we simply call (execute) the ICCS instance:

>>> convergence_list = iccs()  # this runs the ICCS algorithm and returns
>>>                            # a list of the convergence value after each
>>>                            # iteration.
>>> plotstack(iccs, padded=False)
>>>

Despite of the random noise seismogram, the phase arrival is now visible in the stack. Seismograms with low correlation coefficients can automatically be deselected from the calculations by running ICCS again with autoselect=True:

>>> _ = iccs(autoselect=True)
>>> plotstack(iccs, padded=False)
>>>

Seismograms that fit better with their polarity reversed can be flipped automatically by setting autoflip=True:

>>> _ = iccs(autoflip=True)
>>> plotstack(iccs, padded=False)
>>>

To further improve results, you can interactively update the picks and the time window using stack_pick and stack_tw_pick, respectively, and then run the ICCS algorithm again.

Methods:

• __call__ –

  Run the iccs algorithm.

Attributes:

• plot_padding (timedelta) –

  Padding to apply before and after the time window.

• seismograms (Sequence[ICCSSeismogram]) –

  Input seismograms.

• seismograms_ccnorm (list[float]) –

  Returns a list of the normalised cross-correlation coefficients.

• seismograms_for_plotting (list[MiniSeismogram]) –

  Returns the windowed, detrended, normalised, and optionally flipped seismograms.

• seismograms_prepared (list[MiniSeismogram]) –

  Returns the windowed, detrended, normalised, tapered, and optionally flipped seismograms.

• stack (MiniSeismogram) –

  Returns the stacked seismograms (seismograms_prepared).

• stack_for_plotting (MiniSeismogram) –

  Returns the stacked seismograms (seismograms_for_plotting).

• taper_width (timedelta | float) –

  Taper width.

• window_post (timedelta) –

  End of the time window relative to the pick.

• window_pre (timedelta) –

  Begining of the time window relative to the pick.

Source code in pysmo/tools/iccs/_iccs.py

@define
class ICCS:
    """Class to store a list of [`ICCSSeismograms`][pysmo.tools.iccs.ICCSSeismogram] and run the ICCS algorithm.

    The [`ICCS`][pysmo.tools.iccs.ICCS] class serves as a container to store a
    list of seismograms (typically recordings of the same event at different
    stations), and to then run the ICCS algorithm when an instance of this
    class is called. Processing parameters that are common to all seismograms
    are stored as attributes (e.g. time window limits).

    Examples:
        We begin with a set of SAC files of the same event, recorded at different
        stations. All files have a preliminary phase arrival estimate saved in the
        `T0` SAC header, so we can use these files to create instances of the
        [`MiniICCSSeismogram`][pysmo.tools.iccs.MiniICCSSeismogram] class for use
        with the [`ICCS`][pysmo.tools.iccs.ICCS] class:

        ```python
        >>> from pysmo.classes import SAC
        >>> from pysmo.functions import clone_to_mini
        >>> from pysmo.tools.iccs import MiniICCSSeismogram
        >>> from pathlib import Path
        >>>
        >>> sacfiles = sorted(Path("iccs-example/").glob("*.bhz"))
        >>>
        >>> seismograms = []
        >>> for sacfile in sacfiles:
        ...     sac = SAC.from_file(sacfile)
        ...     update = {"t0": sac.timestamps.t0}
        ...     iccs_seismogram = clone_to_mini(MiniICCSSeismogram, sac.seismogram, update=update)
        ...     seismograms.append(iccs_seismogram)
        ...
        >>>
        ```

        To better illustrate the different modes of running the ICCS algorithm, we
        modify the data and picks in the seismograms to make them **worse** than
        they actually are:

        ```python
        >>> from datetime import timedelta
        >>> from copy import deepcopy
        >>> import numpy as np
        >>>
        >>> # change the sign of the data in the first seismogram
        >>> seismograms[0].data *= -1
        >>>
        >>> # move the initial pick 2 seconds earlier in second seismogram
        >>> seismograms[1].t0 += timedelta(seconds=-2)
        >>>
        >>> # move the initial pick 2 seconds later in third seismogram
        >>> seismograms[2].t0 += timedelta(seconds=2)
        >>>
        >>> # create a seismogram with completely random data
        >>> iccs_random: MiniICCSSeismogram = deepcopy(seismograms[-1])
        >>> np.random.seed(42)  # set this for consistent results during testing
        >>> iccs_random.data = np.random.rand(len(iccs_random))
        >>> seismograms.append(iccs_random)
        >>>
        ```

        We can now create an instance of the [`ICCS`][pysmo.tools.iccs.ICCS] class
        and plot the initial stack together with the input seismograms:

        ```python
        >>> from pysmo.tools.iccs import ICCS, plotstack
        >>> iccs = ICCS(seismograms)
        >>> plotstack(iccs, padded=False)
        >>>
        ```

        ![initial stack](../../../examples/tools/iccs/stack_initial.png#only-light){ loading=lazy }
        ![initial stack](../../../examples/tools/iccs/stack_initial_dark.png#only-dark){ loading=lazy }

        The phase emergance is not visible in the stack, and the (absolute)
        correlation coefficents of the the seismograms are not very high. This
        shows the initial picks are not very good and/or that the data are of low
        quality. To run the ICCS algorithm, we simply call (execute) the ICCS
        instance:

        ```python
        >>> convergence_list = iccs()  # this runs the ICCS algorithm and returns
        >>>                            # a list of the convergence value after each
        >>>                            # iteration.
        >>> plotstack(iccs, padded=False)
        >>>
        ```
        ![initial stack](../../../examples/tools/iccs/stack_first_run.png#only-light){ loading=lazy }
        ![initial stack](../../../examples/tools/iccs/stack_first_run_dark.png#only-dark){ loading=lazy }

        Despite of the random noise seismogram, the phase arrival is now visible in
        the stack. Seismograms with low correlation coefficients can automatically
        be deselected from the calculations by running ICCS again with
        `autoselect=True`:


        ```python
        >>> _ = iccs(autoselect=True)
        >>> plotstack(iccs, padded=False)
        >>>
        ```

        ![initial stack](../../../examples/tools/iccs/stack_autoselect.png#only-light){ loading=lazy }
        ![initial stack](../../../examples/tools/iccs/stack_autoselect_dark.png#only-dark){ loading=lazy }

        Seismograms that fit better with their polarity reversed can be flipped
        automatically by setting `autoflip=True`:

        ```python
        >>> _ = iccs(autoflip=True)
        >>> plotstack(iccs, padded=False)
        >>>
        ```

        ![initial stack](../../../examples/tools/iccs/stack_autoflip.png#only-light){ loading=lazy }
        ![initial stack](../../../examples/tools/iccs/stack_autoflip_dark.png#only-dark){ loading=lazy }

        To further improve results, you can interactively update the picks and
        the time window using [`stack_pick`][pysmo.tools.iccs.stack_pick] and
        [`stack_tw_pick`][pysmo.tools.iccs.stack_tw_pick], respectively, and
        then run the ICCS algorithm again.
    """

    seismograms: Sequence[ICCSSeismogram] = field(
        factory=lambda: list[ICCSSeismogram](), validator=_clear_cache_on_update
    )
    """Input seismograms.

    These seismograms provide the input data for ICCS. They are used to store
    processing parameters and create shorter seismograms (based on pick and
    time window) that are then used for cross-correlation. The shorter
    seismograms are created on the fly and then cached within an [`ICCS`]
    [pysmo.tools.iccs.ICCS] instance.
    """

    window_pre: timedelta = field(
        default=timedelta(seconds=-10),
        validator=[validators.lt(timedelta(seconds=0)), _clear_cache_on_update],
    )
    """Begining of the time window relative to the pick."""

    window_post: timedelta = field(
        default=timedelta(seconds=10),
        validator=[validators.gt(timedelta(seconds=0)), _clear_cache_on_update],
    )
    """End of the time window relative to the pick."""

    plot_padding: timedelta = field(
        default=timedelta(seconds=10),
        validator=[validators.gt(timedelta(seconds=0)), _clear_cache_on_update],
    )
    """Padding to apply before and after the time window.

    This padding is *not* used for the cross-correlation."""

    taper_width: timedelta | float = field(
        default=0.0, validator=_clear_cache_on_update
    )
    """Taper width.

    Can be either a timedelta or a float - see [`pysmo.functions.taper()`][pysmo.functions.taper]
    for details.
    """

    # Prepared seismograms and stack are cached to prevent unnecessary
    # processing. Setting the caches to None will force a new calculation.
    _seismograms_prepared: list[MiniSeismogram] | None = field(init=False)
    _seismograms_for_plotting: list[MiniSeismogram] | None = field(init=False)
    _seismograms_ccnorm: list[float] | None = field(init=False)
    _stack: MiniSeismogram | None = field(init=False)
    _stack_for_plotting: MiniSeismogram | None = field(init=False)

    def __attrs_post_init__(self) -> None:
        self._clear_caches()

    def _clear_caches(self) -> None:
        """Clear all cached attributes."""
        self._seismograms_prepared = None
        self._seismograms_for_plotting = None
        self._seismograms_ccnorm = None
        self._stack = None
        self._stack_for_plotting = None

    @property
    def seismograms_prepared(self) -> list[MiniSeismogram]:
        """Returns the windowed, detrended, normalised, tapered, and optionally flipped seismograms."""

        if self._seismograms_prepared is None:
            self._seismograms_prepared = _prepare_seismograms(
                self.seismograms,
                self.window_pre,
                self.window_post,
                self.taper_width,
            )

        return self._seismograms_prepared

    @property
    def seismograms_for_plotting(self) -> list[MiniSeismogram]:
        """Returns the windowed, detrended, normalised, and optionally flipped seismograms."""

        if self._seismograms_for_plotting is None:
            self._seismograms_for_plotting = _prepare_seismograms(
                self.seismograms,
                self.window_pre,
                self.window_post,
                self.taper_width,
                True,
                self.plot_padding,
            )

        return self._seismograms_for_plotting

    @property
    def seismograms_ccnorm(self) -> list[float]:
        """Returns a list of the normalised cross-correlation coefficients."""

        if self._seismograms_ccnorm is None:
            self._seismograms_ccnorm = _calc_ccnorms(
                self.seismograms_prepared, self.stack
            )

        return self._seismograms_ccnorm

    @property
    def stack(self) -> MiniSeismogram:
        """Returns the stacked seismograms ([`seismograms_prepared`][pysmo.tools.iccs.ICCS.seismograms_prepared]).

        The stack is calculated as the average of all seismograms with the
        attribute [`select`][pysmo.tools.iccs.ICCSSeismogram.select] set to
        [`True`][True]. The [`begin_time`][pysmo.MiniSeismogram.begin_time] of
        the returned stack is the average of the [`begin_time`]
        [pysmo.tools.iccs.ICCSSeismogram.begin_time] of the input seismograms.

        Returns:
            Stacked input seismograms.
        """
        if self._stack is not None:
            return self._stack

        self._stack = _create_stack(self.seismograms_prepared, self.seismograms)
        return self._stack

    @property
    def stack_for_plotting(self) -> MiniSeismogram:
        """Returns the stacked seismograms ([`seismograms_for_plotting`][pysmo.tools.iccs.ICCS.seismograms_for_plotting]).

        Returns:
            Stacked input seismograms.
        """
        if self._stack_for_plotting is not None:
            return self._stack_for_plotting

        self._stack_for_plotting = _create_stack(
            self.seismograms_for_plotting, self.seismograms
        )
        return self._stack_for_plotting

    def __call__(
        self,
        autoflip: bool = False,
        autoselect: bool = False,
        convergence_limit: float = 10e-6,
        convergence_method: CONVERGENCE_METHOD = "corrcoef",
        max_iter: int = 10,
        max_shift: timedelta | None = None,
        min_ccnorm: float = 0.5,
        parallel: bool = False,
    ) -> npt.NDArray:
        """Run the iccs algorithm.

        Parameters:
            autoflip: Automatically toggle [`flip`][pysmo.tools.iccs.ICCSSeismogram.flip] attribute of seismograms.
            autoselect: Automatically set `select` attribute to `False` for poor quality seismograms.
            convergence_limit: Convergence limit at which the algorithm stops.
            convergence_method: Method to calculate convergence criterion.
            max_iter: Maximum number of iterations.
            max_shift: Maximum shift in seconds (see [`delay()`][pysmo.tools.signal.delay]).
            min_ccnorm: Minimum normalised cross-correlation coefficient.
                When `autoselect` is `True`, seismograms with correlation
                coefficients below this value are set to `select=False`.
            parallel: Whether to use parallel processing. Setting this to `True`
                will perform the cross-correlation calculations in parallel
                using [`ProcessPoolExecutor`][concurrent.futures.ProcessPoolExecutor].
                In most cases this will *not* be faster.

        Returns:
            convergence: Array of convergence criterion values.
        """
        convergence_list = []

        for _ in range(max_iter):
            prev_stack = clone_to_mini(MiniSeismogram, self.stack)

            if parallel:
                with ProcessPoolExecutor() as executor:
                    for prepared_seismogram, seismogram in zip(
                        self.seismograms_prepared, self.seismograms
                    ):
                        future = executor.submit(
                            delay,
                            self.stack,
                            prepared_seismogram,
                            max_shift=max_shift,
                            abs_max=autoflip,
                        )
                        future.add_done_callback(
                            partial(
                                _update_seismogram_fn,
                                seismogram,
                                autoflip,
                                autoselect,
                                min_ccnorm,
                            )
                        )
            else:
                for prepared_seismogram, seismogram in zip(
                    self.seismograms_prepared, self.seismograms
                ):
                    _delay, _ccnorm = delay(
                        self.stack,
                        prepared_seismogram,
                        max_shift=max_shift,
                        abs_max=autoflip,
                    )

                    _update_seismogram(
                        _delay, _ccnorm, seismogram, autoflip, autoselect, min_ccnorm
                    )

            self._clear_caches()

            convergence = _calc_convergence(self.stack, prev_stack, convergence_method)
            convergence_list.append(convergence)
            if convergence <= convergence_limit:
                break

        return np.array(convergence_list)

plot_padding class-attribute instance-attribute

plot_padding: timedelta = field(
    default=timedelta(seconds=10),
    validator=[
        gt(timedelta(seconds=0)),
        _clear_cache_on_update,
    ],
)

Padding to apply before and after the time window.

This padding is not used for the cross-correlation.

seismograms class-attribute instance-attribute

seismograms: Sequence[ICCSSeismogram] = field(
    factory=lambda: list[ICCSSeismogram](),
    validator=_clear_cache_on_update,
)

Input seismograms.

These seismograms provide the input data for ICCS. They are used to store processing parameters and create shorter seismograms (based on pick and time window) that are then used for cross-correlation. The shorter seismograms are created on the fly and then cached within an ICCS instance.

seismograms_ccnorm property

seismograms_ccnorm: list[float]

Returns a list of the normalised cross-correlation coefficients.

seismograms_for_plotting property

seismograms_for_plotting: list[MiniSeismogram]

Returns the windowed, detrended, normalised, and optionally flipped seismograms.

seismograms_prepared property

seismograms_prepared: list[MiniSeismogram]

Returns the windowed, detrended, normalised, tapered, and optionally flipped seismograms.

stack property

stack: MiniSeismogram

Returns the stacked seismograms (seismograms_prepared).

The stack is calculated as the average of all seismograms with the attribute select set to True. The begin_time of the returned stack is the average of the begin_time of the input seismograms.

Returns:

• MiniSeismogram –

  Stacked input seismograms.

stack_for_plotting property

stack_for_plotting: MiniSeismogram

Returns the stacked seismograms (seismograms_for_plotting).

Returns:

• MiniSeismogram –

  Stacked input seismograms.

taper_width class-attribute instance-attribute

taper_width: timedelta | float = field(
    default=0.0, validator=_clear_cache_on_update
)

Taper width.

Can be either a timedelta or a float - see pysmo.functions.taper() for details.

window_post class-attribute instance-attribute

window_post: timedelta = field(
    default=timedelta(seconds=10),
    validator=[
        gt(timedelta(seconds=0)),
        _clear_cache_on_update,
    ],
)

End of the time window relative to the pick.

window_pre class-attribute instance-attribute

window_pre: timedelta = field(
    default=timedelta(seconds=-10),
    validator=[
        lt(timedelta(seconds=0)),
        _clear_cache_on_update,
    ],
)

Begining of the time window relative to the pick.

__call__

__call__(
    autoflip: bool = False,
    autoselect: bool = False,
    convergence_limit: float = 1e-05,
    convergence_method: CONVERGENCE_METHOD = "corrcoef",
    max_iter: int = 10,
    max_shift: timedelta | None = None,
    min_ccnorm: float = 0.5,
    parallel: bool = False,
) -> NDArray

Run the iccs algorithm.

Parameters:

• autoflip (bool, default: False ) –

  Automatically toggle flip attribute of seismograms.

• autoselect (bool, default: False ) –

  Automatically set select attribute to False for poor quality seismograms.

• convergence_limit (float, default: 1e-05 ) –

  Convergence limit at which the algorithm stops.

• convergence_method (CONVERGENCE_METHOD, default: 'corrcoef' ) –

  Method to calculate convergence criterion.

• max_iter (int, default: 10 ) –

  Maximum number of iterations.

• max_shift (timedelta | None, default: None ) –

  Maximum shift in seconds (see delay()).

• min_ccnorm (float, default: 0.5 ) –

  Minimum normalised cross-correlation coefficient. When autoselect is True, seismograms with correlation coefficients below this value are set to select=False.

• parallel (bool, default: False ) –

  Whether to use parallel processing. Setting this to True will perform the cross-correlation calculations in parallel using ProcessPoolExecutor. In most cases this will not be faster.

Returns:

• convergence ( NDArray ) –

  Array of convergence criterion values.

Source code in pysmo/tools/iccs/_iccs.py

def __call__(
    self,
    autoflip: bool = False,
    autoselect: bool = False,
    convergence_limit: float = 10e-6,
    convergence_method: CONVERGENCE_METHOD = "corrcoef",
    max_iter: int = 10,
    max_shift: timedelta | None = None,
    min_ccnorm: float = 0.5,
    parallel: bool = False,
) -> npt.NDArray:
    """Run the iccs algorithm.

    Parameters:
        autoflip: Automatically toggle [`flip`][pysmo.tools.iccs.ICCSSeismogram.flip] attribute of seismograms.
        autoselect: Automatically set `select` attribute to `False` for poor quality seismograms.
        convergence_limit: Convergence limit at which the algorithm stops.
        convergence_method: Method to calculate convergence criterion.
        max_iter: Maximum number of iterations.
        max_shift: Maximum shift in seconds (see [`delay()`][pysmo.tools.signal.delay]).
        min_ccnorm: Minimum normalised cross-correlation coefficient.
            When `autoselect` is `True`, seismograms with correlation
            coefficients below this value are set to `select=False`.
        parallel: Whether to use parallel processing. Setting this to `True`
            will perform the cross-correlation calculations in parallel
            using [`ProcessPoolExecutor`][concurrent.futures.ProcessPoolExecutor].
            In most cases this will *not* be faster.

    Returns:
        convergence: Array of convergence criterion values.
    """
    convergence_list = []

    for _ in range(max_iter):
        prev_stack = clone_to_mini(MiniSeismogram, self.stack)

        if parallel:
            with ProcessPoolExecutor() as executor:
                for prepared_seismogram, seismogram in zip(
                    self.seismograms_prepared, self.seismograms
                ):
                    future = executor.submit(
                        delay,
                        self.stack,
                        prepared_seismogram,
                        max_shift=max_shift,
                        abs_max=autoflip,
                    )
                    future.add_done_callback(
                        partial(
                            _update_seismogram_fn,
                            seismogram,
                            autoflip,
                            autoselect,
                            min_ccnorm,
                        )
                    )
        else:
            for prepared_seismogram, seismogram in zip(
                self.seismograms_prepared, self.seismograms
            ):
                _delay, _ccnorm = delay(
                    self.stack,
                    prepared_seismogram,
                    max_shift=max_shift,
                    abs_max=autoflip,
                )

                _update_seismogram(
                    _delay, _ccnorm, seismogram, autoflip, autoselect, min_ccnorm
                )

        self._clear_caches()

        convergence = _calc_convergence(self.stack, prev_stack, convergence_method)
        convergence_list.append(convergence)
        if convergence <= convergence_limit:
            break

    return np.array(convergence_list)

ICCSSeismogram

Bases: Seismogram, Protocol

Protocol class to define the ICCSSeismogram type.

The ICCSSeismogram type extends the Seismogram type with the addition of parameters that are required for ICCS.

Methods:

• __len__ –

  The length of the Seismogram.

Attributes:

• begin_time (datetime) –

  Seismogram begin time.

• data (NDArray) –

  Seismogram data.

• delta (timedelta) –

  The sampling interval.

• end_time (datetime) –

  Seismogram end time.

• flip (bool) –

  Data in seismogram should be flipped for ICCS.

• select (bool) –

  Use seismogram to create stack.

• t0 (datetime) –

  Initial pick.

• t1 (datetime | None) –

  Updated pick.

Source code in pysmo/tools/iccs/_types.py

@runtime_checkable
class ICCSSeismogram(Seismogram, Protocol):
    """Protocol class to define the `ICCSSeismogram` type.

    The `ICCSSeismogram` type extends the [`Seismogram`][pysmo.Seismogram] type
    with the addition of parameters that are required for ICCS.
    """

    @property
    def flip(self) -> bool:
        """Data in seismogram should be flipped for ICCS."""
        ...

    @flip.setter
    def flip(self, value: bool) -> None: ...

    @property
    def select(self) -> bool:
        """Use seismogram to create stack."""
        ...

    @select.setter
    def select(self, value: bool) -> None: ...

    @property
    def t0(self) -> datetime:
        """Initial pick."""
        ...

    @t0.setter
    def t0(self, value: datetime) -> None: ...

    @property
    def t1(self) -> datetime | None:
        """Updated pick."""
        ...

    @t1.setter
    def t1(self, value: datetime | None) -> None: ...

begin_time property writable

begin_time: datetime

Seismogram begin time.

data property writable

data: NDArray

Seismogram data.

delta property writable

delta: timedelta

The sampling interval.

end_time property

end_time: datetime

Seismogram end time.

flip property writable

flip: bool

Data in seismogram should be flipped for ICCS.

select property writable

select: bool

Use seismogram to create stack.

t0 property writable

t0: datetime

Initial pick.

t1 property writable

t1: datetime | None

Updated pick.

__len__

__len__() -> int

The length of the Seismogram.

Returns:

• int –

  Number of samples in the data array.

Source code in pysmo/_types/_seismogram.py

def __len__(self) -> int:
    """The length of the Seismogram.

    Returns:
        Number of samples in the data array.
    """
    ...

MiniICCSSeismogram

Minimal implementation of the ICCSSeismogram type.

The MiniICCSSeismogram class provides a minimal implementation of class that is compatible with the ICCSSeismogram.

Examples:

Because ICCSSeismogram inherits from Seismogram, we can easily create MiniICCSSeismogram instances from exisiting seismograms using the clone_to_mini() function, whereby the update parameter is used to provide the extra information needed:

>>> from pysmo.classes import SAC
>>> from pysmo.functions import clone_to_mini
>>> from pysmo.tools.iccs import MiniICCSSeismogram
>>> sac = SAC.from_file("example.sac")
>>> sac_seis = sac.seismogram
>>> update = {"t0": sac.timestamps.t0}
>>> mini_iccs_seis = clone_to_mini(MiniICCSSeismogram, sac_seis, update=update)
>>>

Methods:

• __len__ –

  The length of the Seismogram.

Attributes:

• begin_time (datetime) –

  Seismogram begin time.

• data (NDArray) –

  Seismogram data.

• delta (timedelta) –

  Seismogram sampling interval.

• end_time (datetime) –

  Seismogram end time.

• flip (bool) –

  Data in seismogram should be flipped for ICCS.

• select (bool) –

  Use seismogram to create stack.

• t0 (datetime) –

  Initial pick.

• t1 (datetime | None) –

  Updated pick.

Source code in pysmo/tools/iccs/_types.py

@define(kw_only=True, slots=True)
class MiniICCSSeismogram:
    """Minimal implementation of the [`ICCSSeismogram`][pysmo.tools.iccs.ICCSSeismogram] type.

    The [`MiniICCSSeismogram`][pysmo.tools.iccs.ICCSSeismogram] class provides
    a minimal implementation of class that is compatible with the
    [`ICCSSeismogram`][pysmo.tools.iccs.ICCSSeismogram].

    Examples:
        Because [`ICCSSeismogram`][pysmo.tools.iccs.ICCSSeismogram] inherits
        from [`Seismogram`][pysmo.Seismogram], we can easily create
        [`MiniICCSSeismogram`][pysmo.tools.iccs.MiniICCSSeismogram] instances
        from exisiting seismograms using the
        [`clone_to_mini()`][pysmo.functions.clone_to_mini] function, whereby
        the `update` parameter is used to provide the extra information needed:

        ```python
        >>> from pysmo.classes import SAC
        >>> from pysmo.functions import clone_to_mini
        >>> from pysmo.tools.iccs import MiniICCSSeismogram
        >>> sac = SAC.from_file("example.sac")
        >>> sac_seis = sac.seismogram
        >>> update = {"t0": sac.timestamps.t0}
        >>> mini_iccs_seis = clone_to_mini(MiniICCSSeismogram, sac_seis, update=update)
        >>>
        ```
    """

    begin_time: datetime = field(
        default=SEISMOGRAM_DEFAULTS.begin_time.value, validator=datetime_is_utc
    )
    """Seismogram begin time."""

    delta: timedelta = SEISMOGRAM_DEFAULTS.delta.value
    """Seismogram sampling interval."""

    data: npt.NDArray = field(factory=lambda: np.array([]))
    """Seismogram data."""

    t0: datetime = field(validator=datetime_is_utc)
    """Initial pick."""

    t1: datetime | None = field(
        default=None, validator=validators.optional(datetime_is_utc)
    )
    """Updated pick."""

    flip: bool = False
    """Data in seismogram should be flipped for ICCS."""

    select: bool = True
    """Use seismogram to create stack."""

    def __len__(self) -> int:
        """The length of the Seismogram.

        Returns:
            Number of samples in the data array.
        """
        return np.size(self.data)

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

begin_time class-attribute instance-attribute

begin_time: datetime = field(
    default=begin_time.value, validator=datetime_is_utc
)

Seismogram begin time.

data class-attribute instance-attribute

data: NDArray = field(factory=lambda: array([]))

Seismogram data.

delta class-attribute instance-attribute

delta: timedelta = delta.value

Seismogram sampling interval.

end_time property

end_time: datetime

Seismogram end time.

flip class-attribute instance-attribute

flip: bool = False

Data in seismogram should be flipped for ICCS.

select class-attribute instance-attribute

select: bool = True

Use seismogram to create stack.

t0 class-attribute instance-attribute

t0: datetime = field(validator=datetime_is_utc)

Initial pick.

t1 class-attribute instance-attribute

t1: datetime | None = field(
    default=None, validator=optional(datetime_is_utc)
)

Updated pick.

__len__

__len__() -> int

The length of the Seismogram.

Returns:

• int –

  Number of samples in the data array.

Source code in pysmo/tools/iccs/_types.py

def __len__(self) -> int:
    """The length of the Seismogram.

    Returns:
        Number of samples in the data array.
    """
    return np.size(self.data)

plotstack

plotstack(
    iccs: ICCS,
    padded: bool = True,
    return_fig: bool = False,
) -> tuple[Figure, Axes] | None

Plot the ICCS stack.

Parameters:

• iccs (ICCS) –

  Instance of the ICCS class.

• padded (bool, default: True ) –

  If True, the plot is padded on both sides of the time window by the amount defined in ICCS.plot_padding.

• return_fig (bool, default: False ) –

  If True, the Figure and Axes objects are returned instead of shown.

Returns:

• tuple[Figure, Axes] | None –

  Figure of the stack with the seismograms if return_fig is True.

Examples:

The default plotting mode is to pad the stack beyond the time window used for the cross-correlations (highlighted in light green). This is useful particularly useful for narrow time windows. Note that because of the padding, the displayed stack isn't exactly what is used for the cross-correlations.

>>> from pysmo.tools.iccs import ICCS, plotstack
>>> iccs = ICCS(iccs_seismograms)
>>> _ = iccs(autoselect=True, autoflip=True)
>>>
>>> plotstack(iccs)
>>>

To view the stack exactly as it is used in the cross-correlations, set the padded argument to False:

>>> plotstack(iccs, padded=False)
>>>

Source code in pysmo/tools/iccs/_functions.py

def plotstack(
    iccs: ICCS, padded: bool = True, return_fig: bool = False
) -> tuple[Figure, Axes] | None:
    """Plot the ICCS stack.

    Parameters:
        iccs: Instance of the [`ICCS`][pysmo.tools.iccs.ICCS] class.
        padded: If True, the plot is padded on both sides of the
            time window by the amount defined in
            [`ICCS.plot_padding`][pysmo.tools.iccs.ICCS.plot_padding].
        return_fig: If `True`, the [`Figure`][matplotlib.figure.Figure] and
            [`Axes`][matplotlib.axes.Axes] objects are returned instead of
            shown.

    Returns:
        Figure of the stack with the seismograms if `return_fig` is `True`.

    Examples:
        The default plotting mode is to pad the stack beyond the time window
        used for the cross-correlations (highlighted in light green). This is
        useful particularly useful for narrow time windows. Note that because
        of the padding, the displayed stack isn't exactly what is used for the
        cross-correlations.

        ```python
        >>> from pysmo.tools.iccs import ICCS, plotstack
        >>> iccs = ICCS(iccs_seismograms)
        >>> _ = iccs(autoselect=True, autoflip=True)
        >>>
        >>> plotstack(iccs)
        >>>
        ```

        ![plotstack padded](../../../examples/tools/iccs/plotstack_padded.png#only-light){ loading=lazy }
        ![plotstack padded](../../../examples/tools/iccs/plotstack_padded_dark.png#only-dark){ loading=lazy }

        To view the stack exactly as it is used in the cross-correlations, set
        the `padded` argument to `False`:

        ```python
        >>> plotstack(iccs, padded=False)
        >>>
        ```

        ![plotstack](../../../examples/tools/iccs/plotstack.png#only-light){ loading=lazy }
        ![plotstack](../../../examples/tools/iccs/plotstack_dark.png#only-dark){ loading=lazy }
    """

    fig, ax = _plot_common_stack(iccs, padded)
    if return_fig:
        return fig, ax
    plt.show()
    return None

stack_pick

stack_pick(
    iccs: ICCS,
    padded: bool = True,
    return_fig: bool = False,
) -> tuple[Figure, Axes] | None

Manually pick t1 in the stack and apply it to all seismograms.

This function launches an interactive figure to manually pick a new phase arrival in the stack, and then apply it to all seismograms.

Parameters:

• iccs (ICCS) –

  Instance of the ICCS class.

• padded (bool, default: True ) –

  If True, the plot is padded on both sides of the time window by the amount defined in ICCS.plot_padding.

• return_fig (bool, default: False ) –

  If True, the Figure and Axes objects are returned instead of shown.

Returns:

• tuple[Figure, Axes] | None –

  Figure of the stack with the picker if return_fig is True.

Examples:

>>> from pysmo.tools.iccs import ICCS, stack_pick
>>> iccs = ICCS(iccs_seismograms)
>>> _ = iccs(autoselect=True, autoflip=True)
>>>
>>> stack_pick(iccs)
>>>

Source code in pysmo/tools/iccs/_functions.py

def stack_pick(
    iccs: ICCS, padded: bool = True, return_fig: bool = False
) -> tuple[Figure, Axes] | None:
    """Manually pick [`t1`][pysmo.tools.iccs.ICCSSeismogram.t1] in the stack and apply it to all seismograms.

    This function launches an interactive figure to manually pick a new phase
    arrival in the stack, and then apply it to all seismograms.

    Parameters:
        iccs: Instance of the [`ICCS`][pysmo.tools.iccs.ICCS] class.
        padded: If True, the plot is padded on both sides of the
            time window by the amount defined in
            [`ICCS.plot_padding`][pysmo.tools.iccs.ICCS.plot_padding].
        return_fig: If `True`, the [`Figure`][matplotlib.figure.Figure] and
            [`Axes`][matplotlib.axes.Axes] objects are returned instead of
            shown.

    Returns:
        Figure of the stack with the picker if `return_fig` is `True`.

    Examples:
        ```python
        >>> from pysmo.tools.iccs import ICCS, stack_pick
        >>> iccs = ICCS(iccs_seismograms)
        >>> _ = iccs(autoselect=True, autoflip=True)
        >>>
        >>> stack_pick(iccs)
        >>>
        ```

        ![Stack Pick Figure](../../../examples/tools/iccs/stack_pick.png#only-light){ loading=lazy }
        ![Stack Pick Figure](../../../examples/tools/iccs/stack_pick_dark.png#only-dark){ loading=lazy }
    """

    fig, ax = _plot_common_stack(iccs, padded, figsize=(10, 6), constrained=False)

    fig.subplots_adjust(bottom=0.2, left=0.09, right=1.05, top=0.96)

    pick_line = ax.axvline(0, color="g", linewidth=2)

    cursor = Cursor(  # noqa: F841
        ax, useblit=True, color="r", linewidth=2, horizOn=False
    )

    def onclick(event: MouseEvent) -> Any:
        if event.inaxes is ax:
            pick_line.set_xdata(np.array((event.xdata, event.xdata)))
            fig.canvas.draw()
            fig.canvas.flush_events()

    class SaveOrCancel:
        def save(self, _: Event) -> None:
            pickdelta = timedelta(seconds=pick_line.get_xdata(orig=True)[0])  # type: ignore
            plt.close()
            update_pick(iccs, pickdelta)

        def cancel(self, _: Event) -> None:
            plt.close()

    callback = SaveOrCancel()
    ax_save = fig.add_axes((0.7, 0.05, 0.1, 0.075))
    ax_cancel = fig.add_axes((0.81, 0.05, 0.1, 0.075))
    b_save = Button(ax_save, "Save", color="darkgreen", hovercolor="green")
    b_save.on_clicked(callback.save)
    b_abort = Button(ax_cancel, "Cancel", color="darkred", hovercolor="red")
    b_abort.on_clicked(callback.cancel)

    _ = fig.canvas.mpl_connect("button_press_event", onclick)  # type: ignore

    if return_fig:
        return fig, ax
    plt.show()
    return None

stack_tw_pick

stack_tw_pick(
    iccs: ICCS,
    padded: bool = True,
    return_fig: bool = False,
) -> tuple[Figure, Axes] | None

Pick new time window limits in the stack.

This function launches an interactive figure to pick new values for window_pre and window_post.

Parameters:

• iccs (ICCS) –

  Instance of the ICCS class.

• padded (bool, default: True ) –

  If True, the plot is padded on both sides of the time window by the amount defined in ICCS.plot_padding.

• return_fig (bool, default: False ) –

  If True, the Figure and Axes objects are returned instead of shown.

Returns:

• tuple[Figure, Axes] | None –

  Figure of the stack with the picker if return_fig is True.

Info

The new time window may not be chosen such that the pick lies outside the the window. The picker will therefore automatically correct itself for invalid window choices.

Examples:

>>> from pysmo.tools.iccs import ICCS, stack_tw_pick
>>> iccs = ICCS(iccs_seismograms)
>>> _ = iccs(autoselect=True, autoflip=True)
>>>
>>> stack_tw_pick(iccs)
>>>

Source code in pysmo/tools/iccs/_functions.py

def stack_tw_pick(
    iccs: ICCS, padded: bool = True, return_fig: bool = False
) -> tuple[Figure, Axes] | None:
    """Pick new time window limits in the stack.

    This function launches an interactive figure to pick new values for
    [`window_pre`][pysmo.tools.iccs.ICCS.window_pre] and
    [`window_post`][pysmo.tools.iccs.ICCS.window_post].

    Parameters:
        iccs: Instance of the [`ICCS`][pysmo.tools.iccs.ICCS] class.
        padded: If True, the plot is padded on both sides of the
            time window by the amount defined in
            [`ICCS.plot_padding`][pysmo.tools.iccs.ICCS.plot_padding].
        return_fig: If `True`, the [`Figure`][matplotlib.figure.Figure] and
            [`Axes`][matplotlib.axes.Axes] objects are returned instead of
            shown.

    Returns:
        Figure of the stack with the picker if `return_fig` is `True`.

    Info:
        The new time window may not be chosen such that the pick lies
        outside the the window. The picker will therefore automatically correct
        itself for invalid window choices.

    Examples:
        ```python
        >>> from pysmo.tools.iccs import ICCS, stack_tw_pick
        >>> iccs = ICCS(iccs_seismograms)
        >>> _ = iccs(autoselect=True, autoflip=True)
        >>>
        >>> stack_tw_pick(iccs)
        >>>
        ```

        ![Stack Pick Figure](../../../examples/tools/iccs/stack_tw_pick.png#only-light){ loading=lazy }
        ![Stack Pick Figure](../../../examples/tools/iccs/stack_tw_pick_dark.png#only-dark){ loading=lazy }
    """

    fig, ax = _plot_common_stack(iccs, padded, figsize=(10, 6), constrained=False)

    fig.subplots_adjust(bottom=0.2, left=0.09, right=1.05, top=0.96)

    def onselect(xmin: float, xmax: float) -> None:
        # phase arrival pick must be inside time window
        if xmin >= 0:
            xmin = -iccs.stack.delta.total_seconds()
        if xmax <= 0:
            xmax = iccs.stack.delta.total_seconds()

        # time window should not excced the plot
        xlim_min, xlim_max = ax.get_xlim()
        if xmin < xlim_min:
            xmin = xlim_min
        if xmax > xlim_max:
            xmax = xlim_max

        span.extents = (xmin, xmax)

    span = SpanSelector(
        ax,
        onselect,
        "horizontal",
        useblit=True,
        props=dict(alpha=0.5, facecolor="tab:blue"),
        interactive=True,
        drag_from_anywhere=True,
    )

    span.extents = (
        iccs.window_pre.total_seconds(),
        iccs.window_post.total_seconds(),
    )

    class SaveOrCancel:
        def save(self, _: Event) -> None:
            iccs.window_pre = timedelta(seconds=span.extents[0])
            iccs.window_post = timedelta(seconds=span.extents[1])
            plt.close()

        def cancel(self, _: Event) -> None:
            plt.close()

    callback = SaveOrCancel()
    ax_save = fig.add_axes((0.7, 0.05, 0.1, 0.075))
    ax_cancel = fig.add_axes((0.81, 0.05, 0.1, 0.075))
    b_save = Button(ax_save, "Save", color="darkgreen", hovercolor="green")
    b_save.on_clicked(callback.save)
    b_abort = Button(ax_cancel, "Cancel", color="darkred", hovercolor="red")
    b_abort.on_clicked(callback.cancel)

    if return_fig:
        return fig, ax
    plt.show()
    return None

update_pick

update_pick(iccs: ICCS, pickdelta: timedelta) -> None

Update t1 in all seismograms by the same amount.

Parameters:

• iccs (ICCS) –

  Instance of the ICCS class.

• pickdelta (timedelta) –

  delta applied to all picks.

Source code in pysmo/tools/iccs/_functions.py

def update_pick(iccs: ICCS, pickdelta: timedelta) -> None:
    """Update [`t1`][pysmo.tools.iccs.ICCSSeismogram.t1] in all seismograms by the same amount.

    Parameters:
        iccs: Instance of the [`ICCS`][pysmo.tools.iccs.ICCS] class.
        pickdelta: delta applied to all picks.
    """
    for seismogram in iccs.seismograms:
        current_pick = seismogram.t1 or seismogram.t0
        seismogram.t1 = current_pick + pickdelta
    iccs._clear_caches()  # seismograms and stack need to be refreshed