Source code for openquake.hmtk.seismicity.occurrence.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.

import numpy as np


[docs]def recurrence_table(mag, dmag, year, time_interval=None): """ Table of recurrence statistics for each magnitude [Magnitude, Number of Observations, Cumulative Number of Observations >= M, Number of Observations (normalised to annual value), Cumulative Number of Observations (normalised to annual value)] Counts number and cumulative number of occurrences of each magnitude in catalogue :param numpy.ndarray mag: Catalog matrix magnitude column :param numpy.ndarray dmag: Magnitude interval :param numpy.ndarray year: Catalog matrix year column :returns numpy.ndarray recurrence table: Recurrence table """ # Define magnitude vectors if time_interval is None: num_year = np.max(year) - np.min(year) + 1.0 else: num_year = time_interval upper_m = np.max(np.ceil(10.0 * mag) / 10.0) lower_m = np.min(np.floor(10.0 * mag) / 10.0) mag_range = np.arange(lower_m, upper_m + (1.5 * dmag), dmag) mval = mag_range[:-1] + (dmag / 2.0) # Find number of earthquakes inside range number_obs = np.histogram(mag, mag_range)[0] number_rows = np.shape(number_obs)[0] # Cumulative number of events n_c = np.zeros((number_rows, 1)) i = 0 while i < number_rows: n_c[i] = np.sum(number_obs[i:], axis=0) i += 1 # Normalise to Annual Rate number_obs_annual = number_obs / num_year n_c_annual = n_c / num_year rec_table = np.column_stack( [mval, number_obs, n_c, number_obs_annual, n_c_annual] ) return rec_table
[docs]def input_checks(catalogue, config, completeness): """Performs a basic set of input checks on the data""" if isinstance(completeness, np.ndarray): # completeness table is a numpy array (i.e. [year, magnitude]) if np.shape(completeness)[1] != 2: raise ValueError("Completeness Table incorrectly configured") else: cmag = completeness[:, 1] ctime = completeness[:, 0] elif isinstance(completeness, float): # Completeness corresponds to a single magnitude (i.e. applies to # the entire catalogue) cmag = np.array(completeness) ctime = np.array(np.min(catalogue.data["year"])) else: # Everything is valid - i.e. no completeness magnitude cmag = np.array(np.min(catalogue.data["magnitude"])) ctime = np.array(np.min(catalogue.data["year"])) # Set the config parameters if missing config = {} if config is None else config key = "reference_magnitude" if key not in config or config[key] is None: config["reference_magnitude"] = 0.0 key = "magnitude_interval" if key not in config or config[key] is None: config["magnitude_interval"] = 0.1 ref_mag = config["reference_magnitude"] dmag = config["magnitude_interval"] return cmag, ctime, ref_mag, dmag, config
[docs]def generate_trunc_gr_magnitudes(bval, mmin, mmax, nsamples): """ Generate a random list of magnitudes distributed according to a truncated Gutenberg-Richter model :param float bval: b-value :param float mmin: Minimum Magnitude :param float mmax: Maximum Magnitude :param int nsamples: Number of samples :returns: Vector of generated magnitudes """ sampler = np.random.uniform(0.0, 1.0, nsamples) beta = bval * np.log(10.0) return (-1.0 / beta) * ( np.log(1.0 - sampler * (1 - np.exp(-beta * (mmax - mmin)))) ) + mmin
[docs]def generate_synthetic_magnitudes(aval, bval, mmin, mmax, nyears): """ Generates a synthetic catalogue for a specified number of years, with magnitudes distributed according to a truncated Gutenberg-Richter distribution :param float aval: a-value :param float bval: b-value :param float mmin: Minimum Magnitude :param float mmax: Maximum Magnitude :param int nyears: Number of years :returns: Synthetic catalogue (dict) with year and magnitude attributes """ nsamples = int(np.round(nyears * (10.0 ** (aval - bval * mmin)), 0)) year = np.random.randint(0, nyears, nsamples) # Get magnitudes mags = generate_trunc_gr_magnitudes(bval, mmin, mmax, nsamples) return {"magnitude": mags, "year": np.sort(year)}
# def downsample_completeness_table(comp_table, sample_width=0.1, mmax=None): # """ # Re-sample the completeness table to a specified sample_width # """ # new_comp_table = [] # for i in range(comp_table.shape[0] - 1): # mvals = np.arange(comp_table[i, 1], # comp_table[i + 1, 1], d_m) # FIXME: d_m is undefined! # new_comp_table.extend([[comp_table[i, 0], mval] for mval in mvals]) # # If mmax > last magnitude in completeness table # if mmax and (mmax > comp_table[-1, 1]): # new_comp_table.extend( # [[comp_table[-1, 0], mval] # for mval in np.arange(comp_table[-1, 1], mmax + d_m, d_m)]) # return np.array(new_comp_table)
[docs]def get_completeness_counts(catalogue, completeness, d_m): """ Returns the number of earthquakes in a set of magnitude bins of specified with, along with the corresponding completeness duration (in years) of the bin :param catalogue: Earthquake catalogue as instance of :class: openquake.hmtk.seisimicity.catalogue.Catalogue :param numpy.ndarray completeness: Completeness table [year, magnitude] :param float d_m: Bin size :returns: * cent_mag - array indicating center of magnitude bins * t_per - array indicating total duration (in years) of completeness * n_obs - number of events in completeness period """ mmax_obs = np.max(catalogue.data["magnitude"]) # thw line below was added by Nick Ackerley but it breaks the tests # catalogue.data["dtime"] = catalogue.get_decimal_time() if mmax_obs > np.max(completeness[:, 1]): cmag = np.hstack([completeness[:, 1], mmax_obs]) else: cmag = completeness[:, 1] cyear = np.hstack([catalogue.end_year + 1, completeness[:, 0]]) # When the magnitude value is on the bin edge numpy's histogram function # may assign randomly to one side or the other based on the floating # point value. As catalogues are rounded to the nearest 0.1 this occurs # frequently! So we offset the bin edge by a very tiny amount to ensure # that, for example, M = 4.099999999 is assigned to the bin M = 4.1 and # not 4.0 master_bins = np.arange(np.min(cmag) - 1.0e-7, np.max(cmag) + d_m, d_m) count_rates = np.zeros(len(master_bins) - 1) count_years = np.zeros_like(count_rates) for i in range(len(cyear) - 1): time_idx = np.logical_and( catalogue.data["dtime"] < cyear[i], catalogue.data["dtime"] >= cyear[i + 1], ) nyrs = cyear[i] - cyear[i + 1] sel_mags = catalogue.data["magnitude"][time_idx] m_idx = np.where(master_bins >= (cmag[i] - (d_m / 2.0)))[0] m_bins = master_bins[m_idx] count_rates[m_idx[:-1]] += np.histogram(sel_mags, bins=m_bins)[ 0 ].astype(float) count_years[m_idx[:-1]] += float(nyrs) # Removes any zero rates greater than last_loc = np.where(count_rates > 0)[0][-1] n_obs = count_rates[: (last_loc + 1)] t_per = count_years[: (last_loc + 1)] cent_mag = (master_bins[:-1] + master_bins[1:]) / 2.0 cent_mag = np.around(cent_mag[: (last_loc + 1)], 3) return cent_mag, t_per, n_obs