Source code for openquake.hmtk.seismicity.occurrence.kijko_smit

# -*- 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.


import numpy as np
from openquake.hmtk.seismicity.occurrence.base import (
    SeismicityOccurrence,
    OCCURRENCE_METHODS,
)
from openquake.hmtk.seismicity.occurrence.utils import (
    input_checks,
    recurrence_table,
)
from openquake.hmtk.seismicity.occurrence.aki_maximum_likelihood import (
    AkiMaxLikelihood,
)


[docs]@OCCURRENCE_METHODS.add( "calculate", completeness=True, reference_magnitude=0.0, magnitude_interval=0.1, ) class KijkoSmit(SeismicityOccurrence): """ Class to Implement the Kijko & Smit (2012) algorithm for estimation of a- and b-value """
[docs] def calculate(self, catalogue, config, completeness=None): """ Main function to calculate the a- and b-value """ # Input checks cmag, ctime, ref_mag, dmag, config = input_checks( catalogue, config, completeness ) ival = 0 tolerance = 1e-7 number_intervals = np.shape(ctime)[0] b_est = np.zeros(number_intervals, dtype=float) neq = np.zeros(number_intervals, dtype=float) nyr = np.zeros(number_intervals, dtype=float) for ival in range(0, number_intervals): id0 = np.abs(ctime - ctime[ival]) < tolerance m_c = np.min(cmag[id0]) if ival == 0: id1 = np.logical_and( catalogue.data["year"] >= (ctime[ival] - tolerance), catalogue.data["magnitude"] >= (m_c - tolerance), ) nyr[ival] = float(catalogue.end_year) - ctime[ival] + 1.0 elif ival == number_intervals - 1: id1 = np.logical_and( catalogue.data["year"] < (ctime[ival - 1] - tolerance), catalogue.data["magnitude"] >= (m_c - tolerance), ) nyr[ival] = ctime[ival - 1] - ctime[ival] else: id1 = np.logical_and( catalogue.data["year"] >= (ctime[ival] - tolerance), catalogue.data["year"] < (ctime[ival - 1] - tolerance), ) id1 = np.logical_and( id1, catalogue.data["magnitude"] > (m_c - tolerance) ) nyr[ival] = ctime[ival - 1] - ctime[ival] neq[ival] = np.sum(id1) # Get a- and b- value for the selected events temp_rec_table = recurrence_table( catalogue.data["magnitude"][id1], dmag, catalogue.data["year"][id1], ) aki_ml = AkiMaxLikelihood() b_est[ival] = aki_ml._aki_ml( temp_rec_table[:, 0], temp_rec_table[:, 1], dmag, m_c )[0] ival += 1 total_neq = float(np.sum(neq)) bval = self._harmonic_mean(b_est, neq) sigma_b = bval / np.sqrt(total_neq) aval = self._calculate_a_value(bval, total_neq, nyr, cmag, ref_mag) sigma_a = self._calculate_a_value( bval + sigma_b, total_neq, nyr, cmag, ref_mag ) if not config["reference_magnitude"]: aval = np.log10(aval) sigma_a = np.log10(sigma_a) - aval else: sigma_a = sigma_a - aval return bval, sigma_b, aval, sigma_a
def _harmonic_mean(self, parameters, neq): """ Calculates the Harmonic mean of a vector of parameters """ weight = neq.astype(float) / np.sum(neq) if np.shape(parameters)[0] != np.shape(weight)[0]: raise ValueError("Parameter vector not same shape as weights") else: average_value = np.zeros(np.shape(parameters)[0], dtype=float) id0 = np.logical_not(np.isnan(parameters)) average_value = 1.0 / np.sum(weight[id0] / parameters[id0]) return average_value def _calculate_a_value(self, bval, nvalue, nyr, cmag, ref_mag): """ Calculates the rate of events >= ref_mag using the b-value estimator and Eq. 10 of Kijko & Smit """ denominator = np.sum(nyr * np.exp(-bval * (cmag - ref_mag))) return nvalue / denominator