Source code for openquake.hazardlib.source.base

# The Hazard Library
# Copyright (C) 2012-2023 GEM Foundation
#
# This program 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.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Affero General Public License for more details.
#
# You should have received a copy of the GNU Affero General Public License
# along with this program.  If not, see <http://www.gnu.org/licenses/>.
"""
Module :mod:`openquake.hazardlib.source.base` defines a base class for
seismic sources.
"""
import re
import abc
import zlib
import numpy
from openquake.baselib import general
from openquake.hazardlib import mfd
from openquake.hazardlib.pmf import PMF
from openquake.hazardlib.calc.filters import magstr, split_source
from openquake.hazardlib.geo import Point
from openquake.hazardlib.geo.surface.planar import build_planar, PlanarSurface
from openquake.hazardlib.geo.surface.multi import MultiSurface
from openquake.hazardlib.source.rupture import (
    ParametricProbabilisticRupture, NonParametricProbabilisticRupture,
    EBRupture)


[docs]def get_code2cls(): """ :returns: a dictionary source code -> source class """ dic = {} for cls in general.gen_subclasses(BaseSeismicSource): if hasattr(cls, 'code'): dic[cls.code] = cls return dic
[docs]def is_poissonian(src): """ :returns: True if the underlying source is poissonian, false otherwise """ if src.code == b'F': # multiFault return hasattr(src, 'occur_rates') elif src.code == b'N': # nonParametric return False return True
[docs]def poisson_sample(src, eff_num_ses, seed): """ :param src: a poissonian source :param eff_num_ses: number of stochastic event sets * number of samples :param seed: stochastic seed :yields: triples (rupture, rup_id, num_occurrences) """ rng = numpy.random.default_rng(seed) tom = src.temporal_occurrence_model rupids = src.offset + numpy.arange(src.num_ruptures) if not hasattr(src, 'nodal_plane_distribution'): if src.code == b'F': # multifault s = src.get_sections() for i, rate in enumerate(src.occur_rates): # NB: rng.poisson called inside to save memory num_occ = rng.poisson(rate * tom.time_span * eff_num_ses) if num_occ == 0: # skip continue idxs = src.rupture_idxs[i] if len(idxs) == 1: sfc = s[idxs[0]] else: sfc = MultiSurface([s[idx] for idx in idxs]) hypo = s[idxs[0]].get_middle_point() rup = ParametricProbabilisticRupture( src.mags[i], src.rakes[i], src.tectonic_region_type, hypo, sfc, src.occur_rates[i], tom) yield rup, rupids[i], num_occ else: # simple or complex fault ruptures = list(src.iter_ruptures()) rates = numpy.array([rup.occurrence_rate for rup in ruptures]) occurs = rng.poisson(rates * tom.time_span * eff_num_ses) for rup, rupid, num_occ in zip(ruptures, rupids, occurs): if num_occ: yield rup, rupid, num_occ return # else (multi)point sources and area sources usd = src.upper_seismogenic_depth lsd = src.lower_seismogenic_depth rup_args = [] rates = [] for ps in split_source(src): if not hasattr(ps, 'location'): # unsplit containing a single source [ps] = src lon, lat = ps.location.x, ps.location.y for mag, mag_occ_rate in ps.get_annual_occurrence_rates(): for np_prob, np in ps.nodal_plane_distribution.data: for hc_prob, hc_depth in ps.hypocenter_distribution.data: args = (mag_occ_rate, np_prob, hc_prob, mag, np, lon, lat, hc_depth, ps) rup_args.append(args) rates.append(mag_occ_rate * np_prob * hc_prob) eff_rates = numpy.array(rates) * tom.time_span * eff_num_ses occurs = rng.poisson(eff_rates) for num_occ, args, rupid, rate in zip(occurs, rup_args, rupids, rates): if num_occ: _, np_prob, hc_prob, mag, np, lon, lat, hc_depth, ps = args hc = Point(lon, lat, hc_depth) hdd = numpy.array([(1., hc.depth)]) [[[planar]]] = build_planar( ps.get_planin([(1., mag)], [(1., np)]), hdd, lon, lat, usd, lsd) rup = ParametricProbabilisticRupture( mag, np.rake, ps.tectonic_region_type, hc, PlanarSurface.from_(planar), rate, tom) yield rup, rupid, num_occ
[docs]def timedep_sample(src, eff_num_ses, seed): """ :param src: a time-dependent source :param eff_num_ses: number of stochastic event sets * number of samples :param seed: stochastic seed :yields: triples (rupture, rup_id, num_occurrences) """ rng = numpy.random.default_rng(seed) rupids = src.offset + numpy.arange(src.num_ruptures) if src.code == b'F': # time-dependent multifault s = src.get_sections() for i, probs in enumerate(src.probs_occur): cdf = numpy.cumsum(probs) num_occ = numpy.digitize(rng.random(eff_num_ses), cdf).sum() if num_occ == 0: # ignore non-occurring ruptures continue idxs = src.rupture_idxs[i] if len(idxs) == 1: sfc = s[idxs[0]] else: sfc = MultiSurface([s[idx] for idx in idxs]) hypo = sfc.get_middle_point() pmf = PMF([(p, o) for o, p in enumerate(probs)]) yield (NonParametricProbabilisticRupture( src.mags[i], src.rakes[i], src.tectonic_region_type, hypo, sfc, pmf), rupids[i], num_occ) else: # time-dependent nonparametric mutex_weight = getattr(src, 'mutex_weight', 1) for rup, rupid in zip(src.iter_ruptures(), rupids): occurs = rup.sample_number_of_occurrences(eff_num_ses, rng) if mutex_weight < 1: # consider only the occurrencies below the mutex_weight occurs *= (rng.random(eff_num_ses) < mutex_weight) num_occ = occurs.sum() if num_occ: yield rup, rupid, num_occ
[docs]class BaseSeismicSource(metaclass=abc.ABCMeta): """ Base class representing a seismic source, that is a structure generating earthquake ruptures. :param source_id: Some (numeric or literal) source identifier. Supposed to be unique within the source model. :param name: String, a human-readable name of the source. :param tectonic_region_type: Source's tectonic regime. See :class:`openquake.hazardlib.const.TRT`. """ id = -1 # to be set trt_smr = 0 # set by the engine nsites = 1 # set when filtering the source splittable = True checksum = 0 # set in source_reader weight = 0.001 # set in contexts esites = 0 # updated in estimate_weight offset = 0 # set in fix_src_offset @abc.abstractproperty def MODIFICATIONS(self): pass @property def trt_smrs(self): """ :returns: a list of integers (usually of 1 element) """ trt_smr = self.trt_smr return (trt_smr,) if isinstance(trt_smr, int) else trt_smr
[docs] def serial(self, ses_seed): """ :returns: a random seed derived from source_id and ses_seed """ baseid = re.split('!;', self.source_id)[0] return zlib.crc32(baseid.encode('ascii'), ses_seed)
def __init__(self, source_id, name, tectonic_region_type): self.source_id = source_id self.name = name self.tectonic_region_type = tectonic_region_type self.trt_smr = -1 # set by the engine self.num_ruptures = 0 # set by the engine self.seed = None # set by the engine
[docs] def is_gridded(self): """ :returns: True if the source contains only gridded ruptures """ return False
[docs] @abc.abstractmethod def iter_ruptures(self, **kwargs): """ Get a generator object that yields probabilistic ruptures the source consists of. :returns: Generator of instances of sublclass of :class: `~openquake.hazardlib.source.rupture.BaseProbabilisticRupture`. """
[docs] def sample_ruptures(self, eff_num_ses, ses_seed): """ :param eff_num_ses: number of stochastic event sets * number of samples :yields: triples (rupture, trt_smr, num_occurrences) """ seed = self.serial(ses_seed) sample = poisson_sample if is_poissonian(self) else timedep_sample for rup, rupid, num_occ in sample(self, eff_num_ses, seed): if self.smweight < 1 and hasattr(rup, 'occurrence_rate'): # defined only for poissonian sources # needed to get convergency of the frequency to the rate # tested only in oq-risk-tests etna0 rup.occurrence_rate *= self.smweight ebr = EBRupture(rup, self.id, self.trt_smr, num_occ, rupid) ebr.seed = ebr.id + ses_seed yield ebr
[docs] def get_mags(self): """ :returns: the magnitudes of the ruptures contained in the source """ mags = set() if hasattr(self, 'get_annual_occurrence_rates'): for mag, rate in self.get_annual_occurrence_rates(): mags.add(mag) elif hasattr(self, 'mags'): # MultiFaultSource mags.update(mag for mag in self.mags if mag >= self.min_mag) else: # nonparametric for rup, pmf in self.data: mags.add(rup.mag) return sorted(mags)
[docs] def get_magstrs(self): """ :returns: the magnitudes of the ruptures contained as strings """ if hasattr(self, 'mags'): # MultiFaultSource mags = {magstr(mag) for mag in self.mags} elif hasattr(self, 'data'): # nonparametric mags = {magstr(item[0].mag) for item in self.data} else: mags = {magstr(item[0]) for item in self.get_annual_occurrence_rates()} return sorted(mags)
def __iter__(self): """ Override to implement source splitting """ yield self
[docs] @abc.abstractmethod def count_ruptures(self): """ Return the number of ruptures that will be generated by the source. """
[docs] @abc.abstractmethod def get_min_max_mag(self): """ Return minimum and maximum magnitudes of the ruptures generated by the source. """
[docs] def modify(self, modification, parameters): """ Apply a single modificaton to the source parameters Reflects the modification method and calls it passing ``parameters`` as keyword arguments. Modifications can be applied one on top of another. The logic of stacking modifications is up to a specific source implementation. :param modification: String name representing the type of modification. :param parameters: Dictionary of parameters needed for modification. :raises ValueError: If ``modification`` is missing from the attribute `MODIFICATIONS`. """ if modification not in self.MODIFICATIONS: raise ValueError('Modification %s is not supported by %s' % (modification, type(self).__name__)) meth = getattr(self, 'modify_%s' % modification) meth(**parameters)
[docs] def to_xml(self): """ Convert the source into an XML string, very useful for debugging """ from openquake.hazardlib import nrml, sourcewriter return nrml.to_string(sourcewriter.obj_to_node(self))
def __repr__(self): """ String representation of a source, displaying the source class name and the source id. """ return '<%s %s, weight=%.1f>' % ( self.__class__.__name__, self.source_id, self.weight)
[docs]class ParametricSeismicSource(BaseSeismicSource, metaclass=abc.ABCMeta): """ Parametric Seismic Source generates earthquake ruptures from source parameters, and associated probabilities of occurrence are defined through a magnitude frequency distribution and a temporal occurrence model. :param mfd: Magnitude-Frequency distribution for the source. See :mod:`openquake.hazardlib.mfd`. :param rupture_mesh_spacing: The desired distance between two adjacent points in source's ruptures' mesh, in km. Mainly this parameter allows to balance the trade-off between time needed to compute the :meth:`distance <openquake.hazardlib.geo.surface.base.BaseSurface.get_min_distance>` between the rupture surface and a site and the precision of that computation. :param magnitude_scaling_relationship: Instance of subclass of :class:`openquake.hazardlib.scalerel.base.BaseMSR` to describe how does the area of the rupture depend on magnitude and rake. :param rupture_aspect_ratio: Float number representing how much source's ruptures are more wide than tall. Aspect ratio of 1 means ruptures have square shape, value below 1 means ruptures stretch vertically more than horizontally and vice versa. :param temporal_occurrence_model: Instance of :class:`openquake.hazardlib.tom.PoissonTOM` defining temporal occurrence model for calculating rupture occurrence probabilities :raises ValueError: If either rupture aspect ratio or rupture mesh spacing is not positive (if not None). """ def __init__(self, source_id, name, tectonic_region_type, mfd, rupture_mesh_spacing, magnitude_scaling_relationship, rupture_aspect_ratio, temporal_occurrence_model): super().__init__(source_id, name, tectonic_region_type) if rupture_mesh_spacing is not None and not rupture_mesh_spacing > 0: raise ValueError('rupture mesh spacing must be positive') if rupture_aspect_ratio is not None and not rupture_aspect_ratio > 0: raise ValueError('rupture aspect ratio must be positive') self.mfd = mfd self.rupture_mesh_spacing = rupture_mesh_spacing self.magnitude_scaling_relationship = magnitude_scaling_relationship self.rupture_aspect_ratio = rupture_aspect_ratio self.temporal_occurrence_model = temporal_occurrence_model
[docs] def get_annual_occurrence_rates(self, min_rate=0): """ Get a list of pairs "magnitude -- annual occurrence rate". The list is taken from assigned MFD object (see :meth: `openquake.hazardlib.mfd.base.BaseMFD.get_annual_occurrence_rates`) with simple filtering by rate applied. :param min_rate: A non-negative value to filter magnitudes by minimum annual occurrence rate. Only magnitudes with rates greater than that are included in the result list. :returns: A list of two-item tuples -- magnitudes and occurrence rates. """ scaling_rate = getattr(self, 'scaling_rate', 1) return [(mag, occ_rate * scaling_rate) for mag, occ_rate in self.mfd.get_annual_occurrence_rates() if min_rate is None or occ_rate > min_rate]
[docs] def get_min_max_mag(self): """ Get the minimum and maximum magnitudes of the ruptures generated by the source from the underlying MFD. """ return self.mfd.get_min_max_mag()
[docs] def modify_set_msr(self, new_msr): """ Updates the MSR originally assigned to the source :param new_msr: An instance of the :class:`openquake.hazardlib.scalerel.BaseMSR` """ self.magnitude_scaling_relationship = new_msr
[docs] def modify_set_slip_rate(self, slip_rate: float): """ Updates the slip rate assigned to the source :param slip_rate: The value of slip rate [mm/yr] """ self.slip_rate = slip_rate
[docs] def modify_set_mmax_truncatedGR(self, mmax: float): """ Updates the mmax assigned. This works on for parametric MFDs.s :param mmax: The value of the new maximum magnitude """ # Check that the current src has a TruncatedGRMFD MFD msg = 'This modification works only when the source MFD is a ' msg += 'TruncatedGRMFD' assert self.mfd.__class__.__name__ == 'TruncatedGRMFD', msg self.mfd.max_mag
[docs] def modify_recompute_mmax(self, epsilon: float = 0): """ Updates the value of mmax using the msr and the area of the fault :param epsilon: Number of standard deviations to be added or substracted """ msr = self.magnitude_scaling_relationship area = self.get_fault_surface_area() # area in km^2 mag = msr.get_median_mag(area=area, rake=self.rake) std = msr.get_std_dev_mag(area=area, rake=self.rake) self.mfd.max_mag = mag + epsilon * std
[docs] def modify_adjust_mfd_from_slip(self, slip_rate: float, rigidity: float, constant_term: float = 9.1, recompute_mmax: float = None): """ :param slip_rate: A float defining slip rate [in mm] :param rigidity: A float defining material rigidity [in GPa] :param constant_term: Constant term of the equation used to compute log M0 from magnitude """ # Check that the current src has a TruncatedGRMFD MFD msg = 'This modification works only when the source MFD is a ' msg += 'TruncatedGRMFD' assert self.mfd.__class__.__name__ == 'TruncatedGRMFD', msg # Compute moment area = self.get_fault_surface_area() * 1e6 # area in m^2 rigidity *= 1e9 # rigidity in Pa slip_rate *= 1e-3 # slip rate in m mo = rigidity * area * slip_rate # Update the MFD min_mag = self.mfd.min_mag max_mag = self.mfd.max_mag bin_w = self.mfd.bin_width b_val = self.mfd.b_val self.mfd = mfd.TruncatedGRMFD.from_moment(min_mag, max_mag, bin_w, b_val, mo, constant_term)