# 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 abc
import zlib
import numpy
from openquake.baselib import general
from openquake.hazardlib import mfd
from openquake.hazardlib.geo import Point
from openquake.hazardlib.geo.surface.planar import build_planar, PlanarSurface
from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
[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]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
min_mag = 0 # set in get_oqparams and CompositeSourceModel.filter
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
"""
return zlib.crc32(self.source_id.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
self.min_mag = 0 # set by the SourceConverter
[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)
numpy.random.seed(seed)
for trt_smr in self.trt_smrs:
for rup, num_occ in self._sample_ruptures(eff_num_ses):
rup.seed = seed
if hasattr(rup, 'occurrence_rate'):
# defined only for poissonian sources
rup.occurrence_rate *= self.smweight
seed += 1
yield rup, trt_smr, num_occ
# NB: overridden in MultiFaultSource
def _sample_ruptures(self, eff_num_ses):
tom = getattr(self, 'temporal_occurrence_model', None)
if tom: # time-independent source
yield from self.sample_ruptures_poissonian(eff_num_ses)
else: # time-dependent source (nonparametric)
mutex_weight = getattr(self, 'mutex_weight', 1)
for rup in self.iter_ruptures():
occurs = rup.sample_number_of_occurrences(eff_num_ses)
if mutex_weight < 1:
# consider only the occurrencies below the mutex_weight
occurs *= (numpy.random.random(eff_num_ses) < mutex_weight)
num_occ = occurs.sum()
if num_occ:
yield rup, num_occ
[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, 'source_file'):
# unbound UCERFSource
mags.add(numpy.nan)
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:
if rup.mag >= self.min_mag:
mags.add(rup.mag)
return sorted(mags)
[docs] def sample_ruptures_poissonian(self, eff_num_ses):
"""
:param eff_num_ses: number of stochastic event sets * number of samples
:yields: pairs (rupture, num_occurrences[num_samples])
"""
tom = self.temporal_occurrence_model
if not hasattr(self, 'nodal_plane_distribution'): # fault
ruptures = list(self.iter_ruptures())
rates = numpy.array([rup.occurrence_rate for rup in ruptures])
occurs = numpy.random.poisson(rates * tom.time_span * eff_num_ses)
for rup, num_occ in zip(ruptures, occurs):
if num_occ:
yield rup, num_occ
return
# else (multi)point sources and area sources
usd = self.upper_seismogenic_depth
lsd = self.lower_seismogenic_depth
rup_args = []
rates = []
for src in self:
lon, lat = src.location.x, src.location.y
for mag, mag_occ_rate in src.get_annual_occurrence_rates():
if mag < self.min_mag:
continue
for np_prob, np in src.nodal_plane_distribution.data:
for hc_prob, hc_depth in src.hypocenter_distribution.data:
args = (mag_occ_rate, np_prob, hc_prob,
mag, np, lon, lat, hc_depth, src)
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 = numpy.random.poisson(eff_rates)
for num_occ, args, rate in zip(occurs, rup_args, rates):
if num_occ:
_, np_prob, hc_prob, mag, np, lon, lat, hc_depth, src = args
hc = Point(lon, lat, hc_depth)
hdd = numpy.array([(1., hc.depth)])
[[[planar]]] = build_planar(
src.get_planin([(1., mag)], [(1., np)]),
hdd, lon, lat, usd, lsd)
rup = ParametricProbabilisticRupture(
mag, np.rake, src.tectonic_region_type, hc,
PlanarSurface.from_(planar), rate, tom)
yield rup, num_occ
[docs] @abc.abstractmethod
def get_one_rupture(self, ses_seed, rupture_mutex=False):
"""
Yields one random rupture from a source
"""
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) and
mag >= self.min_mag]
[docs] def get_min_max_mag(self):
"""
Get the minimum and maximum magnitudes of the ruptures generated
by the source from the underlying MFD.
"""
min_mag, max_mag = self.mfd.get_min_max_mag()
return max(self.min_mag, min_mag), max_mag
[docs] def get_one_rupture(self, ses_seed, rupture_mutex=False):
"""
Yields one random rupture from a source. IMPORTANT: this method
does not take into account the frequency of occurrence of the
ruptures
"""
# The Mutex case is admitted only for non-parametric ruptures
msg = 'Mutually exclusive ruptures are admitted only in case of'
msg += ' non-parametric sources'
assert (not rupture_mutex), msg
# Set random seed and get the number of ruptures
num_ruptures = self.count_ruptures()
seed = self.serial(ses_seed)
numpy.random.seed(seed)
idx = numpy.random.choice(num_ruptures)
# NOTE Would be nice to have a method generating a rupture given two
# indexes, one for magnitude and one setting the position
for i, rup in enumerate(self.iter_ruptures()):
if i == idx:
if hasattr(self, 'rup_id'):
rup.seed = self.seed
rup.idx = idx
return rup
[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)