Source code for openquake.hmtk.seismicity.declusterer.dec_afteran

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4

#
# LICENSE
#
# Copyright (C) 2010-2025 GEM Foundation, G. Weatherill, M. Pagani,
# D. Monelli.
#
# The Hazard Modeller's Toolkit is free software: you can redistribute
# it and/or modify it under the terms of the GNU Affero General Public
# License as published by the Free Software Foundation, either version
# 3 of the License, or (at your option) any later version.
#
# You should have received a copy of the GNU Affero General Public License
# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>
#
# DISCLAIMER
#
# The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein
# is released as a prototype implementation on behalf of
# scientists and engineers working within the GEM Foundation (Global
# Earthquake Model).
#
# It is distributed for the purpose of open collaboration and in the
# hope that it will be useful to the scientific, engineering, disaster
# risk and software design communities.
#
# The software is NOT distributed as part of GEM’s OpenQuake suite
# (https://www.globalquakemodel.org/tools-products) and must be considered as a
# separate entity. The software provided herein is designed and implemented
# by scientific staff. It is not developed to the design standards, nor
# subject to same level of critical review by professional software
# developers, as GEM’s OpenQuake software suite.
#
# Feedback and contribution to the software is welcome, and can be
# directed to the hazard scientific staff of the GEM Model Facility
# (hazard@globalquakemodel.org).
#
# The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT
# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# The GEM Foundation, and the authors of the software, assume no
# liability for use of the software.

import numpy as np

from openquake.hmtk.seismicity.declusterer.base import (
    BaseCatalogueDecluster,
    DECLUSTERER_METHODS,
)
from openquake.hmtk.seismicity.utils import decimal_year, haversine
from openquake.hmtk.seismicity.declusterer.distance_time_windows import (
    TIME_DISTANCE_WINDOW_FUNCTIONS,
)


[docs]@DECLUSTERER_METHODS.add(
    "decluster",
    time_distance_window=TIME_DISTANCE_WINDOW_FUNCTIONS,
    time_window=float,
)
class Afteran(BaseCatalogueDecluster):
    """
    This implements the Afteran algorithm as described in this paper:
    Musson, R. (1999), Probabilistic seismic hazard maps for the North
    Balkan Region, Annali Di Geofisica, 42(6), 1109 - 1124
    """

[docs]    def decluster(self, catalogue, config):
        """
        catalogue_matrix, window_opt=TDW_GARDNERKNOPOFF, time_window=60.):

        :param catalogue: a catalogue object
        :type catalogue: Instance of the openquake.hmtk.seismicity.catalogue.Catalogue()
                         class
        :keyword window_opt: method used in calculating distance and time
            windows
        :type window_opt: string
        :keyword time_window: Length (in days) of moving time window
        :type time_window: positive float
        :returns: **vcl vector** indicating cluster number,
                  **flagvector** indicating which earthquakes belong to a
                  cluster
        :rtype: numpy.ndarray
        """
        # Convert time window from days to decimal years
        time_window = config["time_window"] / 365.0
        # Pre-processing steps are the same as for Gardner & Knopoff
        # Get relevent parameters
        mag = catalogue.data["magnitude"]
        neq = np.shape(mag)[0]  # Number of earthquakes
        # Get decimal year (needed for time windows)
        year_dec = decimal_year(
            catalogue.data["year"],
            catalogue.data["month"],
            catalogue.data["day"],
        )
        # Get space windows corresponding to each event
        sw_space, _ = config["time_distance_window"].calc(
            catalogue.data["magnitude"]
        )

        # Pre-allocate cluster index vectors
        vcl = np.zeros(neq, dtype=int)
        flagvector = np.zeros(neq, dtype=int)
        # Rank magnitudes into descending order
        id0 = np.flipud(np.argsort(mag, kind="heapsort"))

        clust_index = 0
        for imarker in id0:
            # Earthquake not allocated to cluster - perform calculation
            if vcl[imarker] == 0:
                # Perform distance calculation
                mdist = haversine(
                    catalogue.data["longitude"],
                    catalogue.data["latitude"],
                    catalogue.data["longitude"][imarker],
                    catalogue.data["latitude"][imarker],
                ).flatten()

                # Select earthquakes inside distance window, later than
                # mainshock and not already assigned to a cluster
                vsel1 = np.where(
                    np.logical_and(
                        vcl == 0,
                        np.logical_and(
                            mdist <= sw_space[imarker],
                            year_dec > year_dec[imarker],
                        ),
                    )
                )[0]
                has_aftershocks = False
                if len(vsel1) > 0:
                    # Earthquakes after event inside distance window
                    temp_vsel1, has_aftershocks = self._find_aftershocks(
                        vsel1, year_dec, time_window, imarker, neq
                    )
                    if has_aftershocks:
                        flagvector[temp_vsel1] = 1
                        vcl[temp_vsel1] = clust_index + 1

                # Select earthquakes inside distance window, earlier than
                # mainshock and not already assigned to a cluster
                has_foreshocks = False
                vsel2 = np.where(
                    np.logical_and(
                        vcl == 0,
                        np.logical_and(
                            mdist <= sw_space[imarker],
                            year_dec < year_dec[imarker],
                        ),
                    )
                )[0]
                if len(vsel2) > 0:
                    # Earthquakes before event inside distance window
                    temp_vsel2, has_foreshocks = self._find_foreshocks(
                        vsel2, year_dec, time_window, imarker, neq
                    )
                    if has_foreshocks:
                        flagvector[temp_vsel2] = -1
                        vcl[temp_vsel2] = clust_index + 1

                if has_aftershocks or has_foreshocks:
                    # Assign mainshock to cluster
                    vcl[imarker] = clust_index + 1
                    clust_index += 1

        return vcl, flagvector

    def _find_aftershocks(self, vsel, year_dec, time_window, imarker, neq):
        """
        Function to identify aftershocks from a set of potential
        events inside the distance window of an earthquake.
        :param vsel: Pointer vector to the location of the events in distance
                     window
        :type vsel: numpy.ndarray
        :param year_dec: Vector of decimal catalogue event times
        :type year_dec: numpy.ndarray
        :param time_window: Moving time window for selection of time clusters
        :type time_window: float
        :param imarker: Index of the mainshock in the catalogue vector
        :type imarker: Integer
        :param neq: Number of events in distance window of mainshock
        :type neq: Integer
        """
        temp_vsel1 = np.zeros(neq, dtype=bool)
        has_aftershocks = False

        # Finds the time difference between events
        delta_time = np.diff(np.hstack([year_dec[imarker], year_dec[vsel]]))
        for iloc in range(0, len(vsel)):
            # If time difference between event is smaller than
            # time window - is an aftershock -> continue

            if delta_time[iloc] < time_window:
                temp_vsel1[vsel[iloc]] = True
                has_aftershocks = True
            else:
                # Time difference between events is larger than
                # window -> no more aftershocks -> return
                return temp_vsel1, has_aftershocks

        return temp_vsel1, has_aftershocks

    def _find_foreshocks(self, vsel, year_dec, time_window, imarker, neq):
        """
        Finds foreshocks from a set of potential events within
        the distance window of a mainshock.
        :param vsel: Pointer vector to the location of the events in distance
                     window
        :type vsel: numpy.ndarray
        :param year_dec: Vector of decimal catalogue event times
        :type year_dec: numpy.ndarray
        :param time_window: Moving time window for selection of time clusters
        :type time_window: float
        :param imarker: Index of the mainshock in the catalogue vector
        :type imarker: Integer
        :param neq: Number of events in distance window of mainshock
        :type neq: Integer
        """

        temp_vsel2 = np.zeros(neq, dtype=bool)
        has_foreshocks = False

        # The initial time is the time of the mainshock
        initial_time = year_dec[imarker]
        year_dec = year_dec[vsel]
        for jloc in range(len(vsel) - 1, -1, -1):
            # If the time between the mainshock and the preceeding
            # event is smaller than the time_window then event
            # is a foreshock

            if (initial_time - year_dec[jloc]) < time_window:
                temp_vsel2[vsel[jloc]] = True
                has_foreshocks = True
                # Update target time to consider current foreshock
                # Then continue
                initial_time = year_dec[jloc]
            else:
                # No events inside time window
                # end of foreshock sequence - return
                return temp_vsel2, has_foreshocks

        return temp_vsel2, has_foreshocks