# The Hazard Library
# Copyright (C) 2012-2026 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.simple_fault` defines
:class:`SimpleFaultSource`.
"""
import copy
import math
from collections import namedtuple
from openquake.baselib.general import round
from openquake.hazardlib import mfd
from openquake.hazardlib.source.base import ParametricSeismicSource
from openquake.hazardlib.geo.surface.simple_fault import SimpleFaultSurface
from openquake.hazardlib.geo.nodalplane import NodalPlane
from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
# Grouping all the "extra" hypo depth options to reduce num args in class
HypoData = namedtuple('HypoData', ['hypo_list', 'slip_list', 'depth_list'],
defaults=[(), (), ()])
[docs]class SimpleFaultSource(ParametricSeismicSource):
"""
Simple fault source typology represents seismicity occurring on a fault
surface with simple geometry.
:param upper_seismogenic_depth:
Minimum depth an earthquake rupture can reach, in km.
:param lower_seismogenic_depth:
Maximum depth an earthquake rupture can reach, in km.
:param fault_trace:
A :class:`~openquake.hazardlib.geo.line.Line` representing
the line of intersection between the fault plane and the Earth's
surface.
:param dip:
Angle between earth surface and fault plane in decimal degrees.
:param rake:
the direction of hanging wall relative to the foot wall.
:param rupture_slip_direction:
Angle describing rupture propagation direction in decimal degrees.
:param HypoData:
Instance used to group the input arguments for the two alternative ways
of specifying the hypocentres within a SimpleFaultSource.
Alternative approach 1 is described by hypo_list and slip_list:
:param hypo_list:
Array describing the relative position of the hypocentre on the
rupture surface. Each line represents a hypocentral position
defined in terms of the relative distance along strike and dip
(from the upper, left corner of the fault surface i.e. the
corner which results from the projection at depth of the first
vertex of the fault trace) and the corresponding weight.
Example 1: one single hypocentral position at the center of the
rupture will be described by the following array[(0.5, 0.5,
1.0)]. Example 2: two possible hypocenters are admitted for a
rupture. One hypocentre is located along the strike at 1/4 of
the fault length and at 1/4 of the fault width along the dip and
occurs with a weight of 0.3, the other one is at 3/4 of fault
length along strike and at 3/4 of fault width along strike with
a weight of 0.7. The numpy array would be entered as
numpy.array([[0.25, 0.25, 0.3], [0.75, 0.75, 0.7]]).
:param slip_list:
Array describing the rupture slip direction, which desribes the
rupture propagation direction on the rupture surface. Each line
represents a rupture slip direction and the corresponding
weight. Example 1: one single rupture slip direction with
angle 90 degree will be described by the following array[(90,
1.0)]. Example 2: two possible rupture slip directions are
admitted for a rupture. One slip direction is at 90 degree
with a weight of 0.7, the other one is at 135 degree with a
weight of 0.3. The numpy array would be entered as
numpy.array([[90, 0.7], [135, 0.3]]).
Alternative approach 2 is described by hypo_depth_list:
:param hypo_depth_list:
List of (probability, depth, fixed_dip_frac) describing a
hypocentre depth distribution. fixed_dip_frac is optional
(None if not set) and overrides the computed down-dip fraction
for that depth if provided. All depths must lie within the
seismogenic depth range or an error is raised. Cannot be
combined with alternative approach 1 (hypo_list/slip_list).
See also :class:`openquake.hazardlib.source.base.ParametricSeismicSource`
for description of other parameters.
:raises ValueError:
If :meth:`~openquake.hazardlib.geo.surface.simple_fault.SimpleFaultSurface.check_fault_data`
fails, if rake value is invalid and if rupture mesh spacing is too high
for the lowest magnitude value.
"""
code = b'S'
MODIFICATIONS = {
'adjust_aspect_ratio',
'adjust_dip',
'adjust_mfd_from_slip',
'recompute_mmax',
'set_aspect_ratio',
'set_bGR',
'set_dip',
'set_geometry',
'set_lower_seismogenic_depth',
'adjust_lower_seismogenic_depth',
'set_mmax_truncatedGR',
'set_msr',
'set_slip_rate',
}
def __init__(self, source_id, name, tectonic_region_type,
mfd, rupture_mesh_spacing,
magnitude_scaling_relationship,
rupture_aspect_ratio,
temporal_occurrence_model,
# simple fault specific parameters
upper_seismogenic_depth, lower_seismogenic_depth,
fault_trace, dip, rake,
hypo_data=None):
super().__init__(
source_id, name, tectonic_region_type, mfd, rupture_mesh_spacing,
magnitude_scaling_relationship, rupture_aspect_ratio,
temporal_occurrence_model)
NodalPlane.check_rake(rake)
SimpleFaultSurface.check_fault_data(
fault_trace, upper_seismogenic_depth, lower_seismogenic_depth,
dip, rupture_mesh_spacing)
self.fault_trace = fault_trace
self.upper_seismogenic_depth = upper_seismogenic_depth
self.lower_seismogenic_depth = lower_seismogenic_depth
self.dip = dip
self.rake = rake
min_mag, _max_mag = self.mfd.get_min_max_mag()
cols_rows = self._get_rupture_dimensions(float('inf'), float('inf'),
min_mag)
if hypo_data is not None:
self.hypo_list = hypo_data.hypo_list
self.slip_list = hypo_data.slip_list
self.hypo_depth_list = tuple(
e if len(e) == 3 else (e[0], e[1], None)
for e in hypo_data.depth_list)
else:
self.hypo_list = ()
self.slip_list = ()
self.hypo_depth_list = ()
if (len(self.hypo_list) and not len(self.slip_list) or
not len(self.hypo_list) and len(self.slip_list)):
raise ValueError('hypo_list and slip_list have to be both given '
'or neither given')
if self.hypo_depth_list and (len(self.hypo_list) or len(self.slip_list)):
raise ValueError(
'hypo_depth_list and hypo_list/slip_list cannot be used together'
)
if self.hypo_depth_list:
total = sum(wei for wei, _depth, _fdf in self.hypo_depth_list)
if abs(total - 1.0) > 1e-7:
raise ValueError(
f'hypo_depth_list probabilities sum to {total}, expected 1.0')
for _wei, depth, fdf in self.hypo_depth_list:
if not (upper_seismogenic_depth <= depth <= lower_seismogenic_depth):
raise ValueError(
f'hypo_depth_list {depth} km is outside the '
f'seismogenic zone [{upper_seismogenic_depth}, '
f'{lower_seismogenic_depth}] km')
if fdf is not None and not (0.0 < fdf <= 1.0):
raise ValueError(
f'hypo_depth_list fixedDipFrac {fdf} for depth {depth} km '
f'must be in the range [0, 1]')
if 1 in cols_rows:
raise ValueError('mesh spacing %s is too high to represent '
'ruptures of magnitude %s' %
(rupture_mesh_spacing, min_mag))
def _hypo_list_from_depths(self, first_row, rup_rows):
"""
Convert the hypo depth distribution to (dip_frac, weight) pairs
for a given floating rupture.
Depths outside [rup_top_depth, rup_bot_depth] are discarded and
the remaining weights are renormalised to sum to 1.
If fixedDipFrac is specified in the hypoDepthDist then it is
assigned to the corresponding depth instead of a fraction computed
using the floating rupture's depth range.
NOTE: Depths that resolve to the same dip_frac (whether from
fixedDipFrac or coincidental computed equality) are collapsed into
one pair with their weights summed.
NOTE: An error is raised if none of the hypo depths are within the
depth range of a given floating rup.
"""
# Compute top and bottom depth of this floating rupture
sin_dip = math.sin(math.radians(self.dip))
down_dip_delta = self.rupture_mesh_spacing * sin_dip
top_depth = (
self.upper_seismogenic_depth +
first_row * down_dip_delta)
bot_depth = (
self.upper_seismogenic_depth +
(first_row + rup_rows - 1) * down_dip_delta)
rup_depth_range = bot_depth - top_depth
# Keep hypocentres inside given floating rupture's depth range
weight_by_frac = {}
for prob, depth, fdf in self.hypo_depth_list:
if top_depth <= depth <= bot_depth:
# Use fixedDipFrac if provided and otherwise compute from depth
dip_frac = (fdf if fdf is not None
else (depth - top_depth) / rup_depth_range)
# Collapse depths with the same fraction summing their weights
weight_by_frac[dip_frac] = (
weight_by_frac.get(dip_frac, 0.0) + prob)
if not weight_by_frac:
raise ValueError(
f'No hypo_depth_list depths fall in depth range of floated '
f'rupture [{top_depth:.2f}, {bot_depth:.2f}] km')
total = sum(weight_by_frac.values())
return [(frac, w / total) for frac, w in weight_by_frac.items()]
def _iter_ruptures_hypo_depth(self, surface, occurrence_rate, mag,
first_row, rup_rows):
"""
Yield ruptures for given floating position on fault surface using
the absolute depth distribution in hypo_depth_list.
"""
# Iter over the depths in the hypo depth dist
for dip_frac, w in self._hypo_list_from_depths(first_row, rup_rows):
# Use the down-dip fraction to get the hypo on given floating rup
hypo = surface.get_hypo_location(
self.rupture_mesh_spacing, (0.5, dip_frac))
# Make the rup
yield ParametricProbabilisticRupture(
mag, self.rake, self.tectonic_region_type,
hypo, surface, occurrence_rate * w,
self.temporal_occurrence_model)
def _iter_ruptures_hypo_slip(self, surface, occurrence_rate, mag):
"""
Yield ruptures for given floating position on fault surface using
hypo_list and slip_list.
"""
for hypo in self.hypo_list:
for slip in self.slip_list:
hypocenter = surface.get_hypo_location(
self.rupture_mesh_spacing, hypo[:2])
yield ParametricProbabilisticRupture(
mag, self.rake, self.tectonic_region_type,
hypocenter, surface, occurrence_rate * hypo[2] * slip[1],
self.temporal_occurrence_model, slip[0])
[docs] def iter_ruptures(self, **kwargs):
"""
See :meth:
`openquake.hazardlib.source.base.BaseSeismicSource.iter_ruptures`.
Generates a ruptures using the "floating" algorithm: for all the
magnitude values of assigned MFD calculates the rupture size with
respect to MSR and aspect ratio and then places ruptures of that
size on the surface of the whole fault source. The occurrence
rate of each of those ruptures is the magnitude occurrence rate
divided by the number of ruptures that can be placed in a fault.
"""
step = kwargs.get('step', 1)
whole_fault_surface = SimpleFaultSurface.from_fault_data(
self.fault_trace, self.upper_seismogenic_depth,
self.lower_seismogenic_depth, self.dip, self.rupture_mesh_spacing)
whole_fault_mesh = whole_fault_surface.mesh
mesh_rows, mesh_cols = whole_fault_mesh.shape
fault_length = float((mesh_cols - 1) * self.rupture_mesh_spacing)
fault_width = float((mesh_rows - 1) * self.rupture_mesh_spacing)
for mag, mag_occ_rate in self.get_annual_occurrence_rates():
rup_cols, rup_rows = self._get_rupture_dimensions(
fault_length, fault_width, mag)
num_rup_along_length = mesh_cols - rup_cols + 1
num_rup_along_width = mesh_rows - rup_rows + 1
num_rup = num_rup_along_length * num_rup_along_width
occurrence_rate = mag_occ_rate / float(num_rup)
for first_row in range(num_rup_along_width)[::step]:
for first_col in range(num_rup_along_length)[::step]:
mesh = whole_fault_mesh[first_row: first_row + rup_rows,
first_col: first_col + rup_cols]
surface = SimpleFaultSurface(mesh)
if self.hypo_depth_list:
yield from self._iter_ruptures_hypo_depth(
surface, occurrence_rate, mag, first_row, rup_rows)
elif len(self.hypo_list) and len(self.slip_list):
yield from self._iter_ruptures_hypo_slip(
surface, occurrence_rate, mag)
else:
# Regular version using surf centroid
yield ParametricProbabilisticRupture(
mag, self.rake, self.tectonic_region_type,
mesh.get_middle_point(), surface, occurrence_rate,
self.temporal_occurrence_model)
[docs] def get_fault_surface_area(self):
"""
Computes the area covered by the surface of the fault.
:returns:
A float defining the area of the surface of the fault [km^2]
"""
sfc = SimpleFaultSurface.from_fault_data(
self.fault_trace,
self.upper_seismogenic_depth, self.lower_seismogenic_depth,
self.dip, 1.0)
return sfc.get_area()
[docs] def count_ruptures(self):
"""
See :meth:
`openquake.hazardlib.source.base.BaseSeismicSource.count_ruptures`.
"""
try:
whole_fault_surface = SimpleFaultSurface.from_fault_data(
self.fault_trace, self.upper_seismogenic_depth,
self.lower_seismogenic_depth, self.dip,
self.rupture_mesh_spacing)
except Exception as exc:
source_id = self.source_id
exc.args = (f'in {source_id=}: {exc.args[0]}',) + exc.args[1:]
raise
whole_fault_mesh = whole_fault_surface.mesh
mesh_rows, mesh_cols = whole_fault_mesh.shape
fault_length = float((mesh_cols - 1) * self.rupture_mesh_spacing)
fault_width = float((mesh_rows - 1) * self.rupture_mesh_spacing)
self._nr = []
n_hypo = len(self.hypo_depth_list) or len(self.hypo_list) or 1
n_slip = len(self.slip_list) or 1
for (mag, mag_occ_rate) in self.get_annual_occurrence_rates():
if mag_occ_rate == 0:
continue
rup_cols, rup_rows = self._get_rupture_dimensions(
fault_length, fault_width, mag)
num_rup_along_length = mesh_cols - rup_cols + 1
num_rup_along_width = mesh_rows - rup_rows + 1
self._nr.append(num_rup_along_length * num_rup_along_width *
n_hypo * n_slip)
counts = sum(self._nr)
return counts
def _get_rupture_dimensions(self, fault_length, fault_width, mag):
"""
Calculate rupture dimensions for a given magnitude.
:param fault_length:
The length of the fault as a sum of all segments, in km.
:param fault_width:
The width of the fault, in km.
:param mag:
Magnitude value to calculate rupture geometry for.
:returns:
A tuple of two integer items, representing rupture's dimensions:
number of mesh points along length and along width respectively.
The rupture is reshaped (conserving area, if possible) if one
of dimensions exceeds fault geometry. If both do, the rupture
is considered to cover the whole fault.
"""
area = self.magnitude_scaling_relationship.get_median_area(
mag, self.rake)
rup_length = math.sqrt(area * self.get_aspect_ratio(mag))
rup_width = area / rup_length
# clip rupture's length and width to fault's length and width
# if there is no way to fit the rupture otherwise
if area >= fault_length * fault_width:
rup_length = fault_length
rup_width = fault_width
# reshape rupture (conserving area) if its length or width
# exceeds fault's length or width
elif rup_width > fault_width:
rup_length = rup_length * (rup_width / fault_width)
rup_width = fault_width
elif rup_length > fault_length:
rup_width = rup_width * (rup_length / fault_length)
rup_length = fault_length
# round rupture dimensions with respect to mesh_spacing
# and compute number of points in the rupture along length
# and width (aka strike and dip)
rup_cols = int(round(rup_length / self.rupture_mesh_spacing) + 1)
rup_rows = int(round(rup_width / self.rupture_mesh_spacing) + 1)
return rup_cols, rup_rows
[docs] def modify_set_geometry(self, fault_trace, upper_seismogenic_depth,
lower_seismogenic_depth, dip, spacing):
"""
Modifies the current source geometry including trace, seismogenic
depths and dip
"""
# Check the new geometries are valid
SimpleFaultSurface.check_fault_data(
fault_trace, upper_seismogenic_depth, lower_seismogenic_depth,
dip, spacing
)
self.fault_trace = fault_trace
self.upper_seismogenic_depth = upper_seismogenic_depth
self.lower_seismogenic_depth = lower_seismogenic_depth
self.dip = dip
self.rupture_mesh_spacing = spacing
[docs] def modify_adjust_dip(self, increment):
"""
Modifies the dip by an incremental value
:param float increment:
Value by which to increase or decrease the dip (the resulting
dip must still be within 0.0 to 90.0 degrees)
"""
SimpleFaultSurface.check_fault_data(
self.fault_trace, self.upper_seismogenic_depth,
self.lower_seismogenic_depth, self.dip + increment,
self.rupture_mesh_spacing)
self.dip += increment
[docs] def modify_set_lower_seismogenic_depth(self, lsd):
"""
Modifies the lower seismogenic depth
:param float lsd:
New value of the lsd [km]
"""
SimpleFaultSurface.check_fault_data(
self.fault_trace, self.upper_seismogenic_depth, lsd,
self.dip, self.rupture_mesh_spacing)
self.lower_seismogenic_depth = lsd
[docs] def modify_adjust_lower_seismogenic_depth(self, increment):
"""
Modifies the lower seismogenic depth by adding an increment
:param float increment:
Value (in km) by which to increase or decrease the lower
seismogenic depth
"""
lsd = self.lower_seismogenic_depth + increment
SimpleFaultSurface.check_fault_data(
self.fault_trace, self.upper_seismogenic_depth, lsd,
self.dip, self.rupture_mesh_spacing)
self.lower_seismogenic_depth = lsd
[docs] def modify_set_dip(self, dip):
"""
Modifies the dip to the specified value
:param float dip:
New value of dip (must still be within 0.0 to 90.0 degrees)
"""
SimpleFaultSurface.check_fault_data(
self.fault_trace, self.upper_seismogenic_depth,
self.lower_seismogenic_depth, dip, self.rupture_mesh_spacing)
self.dip = dip
def __iter__(self):
mag_rates = self.get_annual_occurrence_rates()
if len(mag_rates) == 1: # not splittable
yield self
return
self.num_ruptures
for i, (mag, rate) in enumerate(mag_rates):
# This is needed in order to reproduce the logic in the
# `rupture_count` method
if rate == 0:
continue
src = copy.copy(self)
src.mfd = mfd.ArbitraryMFD([mag], [rate])
src._num_ruptures = self._nr[i]
yield src
@property
def polygon(self):
"""
The underlying polygon
`"""
return SimpleFaultSurface.surface_projection_from_fault_data(
self.fault_trace, self.upper_seismogenic_depth,
self.lower_seismogenic_depth, self.dip)
[docs] def wkt(self):
"""
:returns: the geometry as a WKT string
"""
return self.polygon.wkt
