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

# -*- 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 warnings
import numpy as np
from openquake.hmtk.seismicity.occurrence.base import (
    SeismicityOccurrence, OCCURRENCE_METHODS)
from openquake.hmtk.seismicity.occurrence.utils import recurrence_table, input_checks


[docs]@OCCURRENCE_METHODS.add('calculate', completeness=True) class AkiMaxLikelihood(SeismicityOccurrence):
[docs] def calculate(self, catalogue, config=None, completeness=None): """ Calculation of b-value and its uncertainty for a given catalogue, using the maximum likelihood method of Aki (1965), with a correction for discrete bin width (Bender, 1983). :param catalogue: See :class:`openquake.hmtk.seismicity.occurrence.base.py` for further explanation :param config: The configuration in this case do not contains specific information :keyword float completeness: Completeness magnitude :return float bval: b-value of the Gutenberg-Richter relationship :return float sigma_b: Standard deviation of the GR b-value """ # Input checks _cmag, _ctime, _ref_mag, dmag, config = input_checks(catalogue, config, completeness) rt = recurrence_table( catalogue.data['magnitude'], dmag, catalogue.data['year']) bval, sigma_b = self._aki_ml(rt[:, 0], rt[:, 1]) return bval, sigma_b
def _aki_ml(self, mval, number_obs, dmag=0.1, m_c=0.0): """ :param numpy.ndarray mval: array of reference magnitudes (column 0 from recurrence table) :param numpy.ndarray number_obs: number of observations in magnitude bin (column 1 from recurrence table) :keyword float dmag: magnitude interval :keyword float m_c: completeness magnitude :return float bval: b-value of the Gutenberg-Richter relationship :return float sigma_b: Standard deviation of the GR b-value """ # Exclude data below Mc id0 = mval >= m_c mval = mval[id0] number_obs = number_obs[id0] # Get Number of events, minimum magnitude and mean magnitude neq = np.sum(number_obs) if neq <= 1: # Cannot determine b-value (too few event) return NaNs warnings.warn('Too few events (<= 1) to calculate b-value') return np.nan, np.nan m_min = np.min(mval) m_ave = np.sum(mval * number_obs) / neq # Calculate b-value bval = np.log10(np.exp(1.0)) / (m_ave - m_min + (dmag / 2.)) # Calculate sigma b from Bender estimator sigma_b = np.sum(number_obs * ((mval - m_ave) ** 2.0)) /\ (neq * (neq - 1)) sigma_b = np.log(10.) * (bval ** 2.0) * np.sqrt(sigma_b) return bval, sigma_b