Source code for openquake.hmtk.seismicity.smoothing.utils

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

#
# LICENSE
#
# Copyright (C) 2010-2023 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.

"""
Module :mod:`openquake.hmtk.seismicity.smoothing.utils` implements
utility functions for smoothed seismicity analysis
"""

import numpy as np


[docs]def hermann_adjustment_factors(bval, min_mag, mag_inc):
    """
    Returns the adjustment factors (fval, fival) proposed by Hermann (1978)

    :param float bval:
        Gutenberg & Richter (1944) b-value

    :param np.ndarray min_mag:
        Minimum magnitude of completeness table

    :param non-negative float mag_inc:
        Magnitude increment of the completeness table
    """

    fval = 10.0 ** (bval * min_mag)
    fival = 10.0 ** (bval * (mag_inc / 2.0)) - 10.0 ** (
        -bval * (mag_inc / 2.0)
    )
    return fval, fival


[docs]def incremental_a_value(bval, min_mag, mag_inc):
    """
    Incremental a-value from cumulative - using the version of the
    Hermann (1979) formula described in Wesson et al. (2003)

    :param float bval:
        Gutenberg & Richter (1944) b-value

    :param np.ndarray min_mag:
        Minimum magnitude of completeness table

    :param float mag_inc:
        Magnitude increment of the completeness table
    """
    a_cum = 10.0 ** (bval * min_mag)
    a_inc = a_cum + np.log10(
        (10.0 ** (bval * mag_inc)) - (10.0 ** (-bval * mag_inc))
    )

    return a_inc


[docs]def get_weichert_factor(beta, cmag, cyear, end_year):
    """
    Gets the Weichert adjustment factor for each the magnitude bins

    :param float beta:
        Beta value of Gutenberg & Richter parameter (b * log(10.))

    :param np.ndarray cmag:
        Magnitude values of the completeness table

    :param np.ndarray cyear:
        Year values of the completeness table

    :param float end_year:
        Last year for consideration in the catalogue

    :returns:
        Weichert adjustment factor (float)
    """
    if len(cmag) > 1:
        # cval corresponds to the mid-point of the completeness bins
        # In the original code it requires that the magnitude bins be
        # equal sizedclass IsotropicGaussian(BaseSmoothingKernel):
        dmag = (cmag[1:] + cmag[:-1]) / 2.0
        cval = np.hstack([dmag, cmag[-1] + (dmag[-1] - cmag[-2])])
    else:
        # Single completeness value so Weichert factor is unity
        return 1.0 / (end_year - cyear[0] + 1), None

    t_f = sum(np.exp(-beta * cval)) / sum(
        (end_year - cyear + 1) * np.exp(-beta * cval)
    )
    return t_f, cval


[docs]def check_completeness_table(completeness_table, catalogue):
    """
    Check to ensure completeness table is in the correct format
    `completeness_table = np.array([[year_, mag_i]]) for i in number of bins`

    :param np.ndarray completeness_table:
        Completeness table in format [[year, mag]]

    :param catalogue:
        Instance of openquake.hmtk.seismicity.catalogue.Catalogue class

    :returns:
        Correct completeness table

    """
    if isinstance(completeness_table, np.ndarray):
        assert np.shape(completeness_table)[1] == 2
        return completeness_table
    elif isinstance(completeness_table, list):
        # Assuming list has only two elements
        assert len(completeness_table) == 2
        return np.array([[completeness_table[0], completeness_table[1]]])
    else:
        # Accepts the minimum magnitude and earliest year of the catalogue
        return np.array(
            [
                [
                    np.min(catalogue.data["year"]),
                    np.min(catalogue.data["magnitude"]),
                ]
            ]
        )


[docs]def get_even_magnitude_completeness(completeness_table, catalogue=None):
    """
    To make the magnitudes evenly spaced, render to a constant 0.1
    magnitude unit

    :param np.ndarray completeness_table:
        Completeness table in format [[year, mag]]

    :param catalogue:
        Instance of openquake.hmtk.seismicity.catalogue.Catalogue class

    :returns:
        Correct completeness table

    """
    mmax = np.floor(10.0 * np.max(catalogue.data["magnitude"])) / 10.0
    check_completeness_table(completeness_table, catalogue)
    cmag = np.hstack([completeness_table[:, 1], mmax + 0.1])
    cyear = np.hstack([completeness_table[:, 0], completeness_table[-1, 0]])
    if np.shape(completeness_table)[0] == 1:
        # Simple single-valued table
        return completeness_table, 0.1

    for iloc in range(0, len(cmag) - 1):
        mrange = np.arange(
            np.floor(10.0 * cmag[iloc]) / 10.0,
            (np.ceil(10.0 * cmag[iloc + 1]) / 10.0),
            0.1,
        )
        temp_table = np.column_stack(
            [cyear[iloc] * np.ones(len(mrange), dtype=float), mrange]
        )
        if iloc == 0:
            completeness_table = np.copy(temp_table)
        else:
            completeness_table = np.vstack([completeness_table, temp_table])
    # completeness_table = np.vstack([completeness_table,
    #    np.array([[cyear[-1], cmag[-1]]])])
    return completeness_table, 0.1