# 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.area` defines :class:`AreaSource`.
"""
import math
from openquake.hazardlib import geo, mfd
from openquake.hazardlib.source.point import PointSource
from openquake.hazardlib.source.base import ParametricSeismicSource
from openquake.hazardlib.mfd.truncated_gr import TruncatedGRMFD
from openquake.hazardlib.mfd.alternative_characteristic_mfd import AlternativeCharacteristicMFD
from openquake.hazardlib.mfd.evenly_discretized import EvenlyDiscretizedMFD
[docs]class AreaSource(ParametricSeismicSource):
"""
Area source represents uniform seismicity occurring over a geographical
region.
:param polygon:
An instance of :class:`openquake.hazardlib.geo.polygon.Polygon`
that defines source's area.
:param area_discretization:
Float number, polygon area discretization spacing in kilometers.
See :meth:`openquake.hazardlib.source.area.AreaSource.iter_ruptures`.
Other parameters (except ``location``) are the same as for
:class:`~openquake.hazardlib.source.point.PointSource`.
"""
code = b'A'
MODIFICATIONS = {
'adjust_aspect_ratio',
'set_aspect_ratio',
'set_geometry',
'set_lower_seismogenic_depth',
'set_upper_seismogenic_depth',
'adjust_lower_seismogenic_depth',
'adjust_upper_seismogenic_depth',
'set_msr',
'set_hypo_depth_dist',
'set_recurrow',
}
def __init__(self, source_id, name, tectonic_region_type,
mfd, rupture_mesh_spacing,
magnitude_scaling_relationship, rupture_aspect_ratio,
temporal_occurrence_model,
# point-specific parameters (excluding location)
upper_seismogenic_depth, lower_seismogenic_depth,
nodal_plane_distribution, hypocenter_distribution,
# area-specific parameters
polygon, area_discretization):
super().__init__(
source_id, name, tectonic_region_type, mfd, rupture_mesh_spacing,
magnitude_scaling_relationship, rupture_aspect_ratio,
temporal_occurrence_model)
self.upper_seismogenic_depth = upper_seismogenic_depth
self.lower_seismogenic_depth = lower_seismogenic_depth
self.nodal_plane_distribution = nodal_plane_distribution
self.hypocenter_distribution = hypocenter_distribution
self.hypo_dip_fracs = getattr(
hypocenter_distribution, 'hypo_dip_fracs', None)
self.polygon = polygon
self.area_discretization = area_discretization
self.max_radius = 0
[docs] def modify_set_geometry(self, polygon):
"""
Modifies the current source geometry by replacing the original polygon
defining the boundary of the area
"""
self.polygon = polygon
[docs] def modify_set_lower_seismogenic_depth(self, lsd):
"""
Modifies the current source geometry by replacing the original
lower seismogenic depth with the passed depth.
Hypocenter depths deeper than the new LSD are dropped and the
remaining probabilities are renormalised (see
ParamatricSeismicSource._shrink_hypo_depths_to_lsd).
"""
self.lower_seismogenic_depth = lsd
self._shrink_hypo_depths_to_lsd(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
"""
self.lower_seismogenic_depth += increment
self._shrink_hypo_depths_to_lsd(self.lower_seismogenic_depth)
[docs] def modify_set_upper_seismogenic_depth(self, usd):
"""
Modifies the current source geometry by replacing the original
upper seismogenic depth with the passed depth
"""
self.upper_seismogenic_depth = usd
[docs] def modify_adjust_upper_seismogenic_depth(self, increment):
"""
Modifies the upper seismogenic depth by adding an increment
:param float increment:
Value (in km) by which to increase or decrease the upper
seismogenic depth
"""
self.upper_seismogenic_depth += increment
[docs] def modify_set_recurrow(self, recurrow):
"""
Modify the recurrence parameters by values given in a dict.
Params are read from the recurrow dict, preferring bg-prefixed keys
e.g., bg_b_value, bg_ref_mag, bg_rate) so a single alt3-style
recurRow branch can carry different but correlated bg and fault
parameters; otherwise the unprefixed alt1/alt2 keys are used.
How it works is determined by if rate_split_bg_frac is present:
- Alt1/Alt3 (no rate_split_bg_frac): Build the TE or AC MFD from the row
params and use it as it is.
- Alt2 (rate_split_bg_frac present): Build the same MFD, then
piecewise-scale bins at or above Mmax-1 by rate_split_bg_frac (bg
side of the alt2 partition between the area source and its faults
and make an EvenlyDiscretizedMFD.
NOTE: This is currently only intended for support of the BC Hydro
NVA SSC logic tree. We may expand it to become a more general
capability.
:param recurrow:
Dict of values to use in given type of MFD e.g.b_value, ref_mag,
rate). Alt3-style rows use prefixes of "bg" for the same keys.
"""
# Constants for BC Hydro NVA model
b_ac = 0.3
delta_mac = 1.0
gamma_eff = 0.9185 # Corrected gamma_eff from eq 1.2 of BCHydro AC memo
bin_width = 0.1
# Get params from recurRow, preferring bg-prefixed keys (alt3)
mmax = self.mmax
recur_model = self.recur_model
bval = float(recurrow.get("bg_b_value", recurrow.get("b_value")))
ref_mag = float(recurrow.get("bg_ref_mag", recurrow.get("ref_mag")))
rate = float(recurrow.get("bg_rate", recurrow.get("rate")))
# Build the MFD
if recur_model == "TE":
a_val = math.log10(rate) + bval * ref_mag
self.mfd = TruncatedGRMFD(
min_mag=ref_mag, max_mag=mmax,
bin_width=bin_width, a_val=a_val, b_val=bval,
)
else:
assert recur_model == "AC"
self.mfd = AlternativeCharacteristicMFD(
min_mag=ref_mag, max_mag=mmax, b_GR=bval, b_AC=b_ac,
bin_width=bin_width, gamma=gamma_eff, delta_m_AC=delta_mac,
total_rate=rate,
)
# Piecewise partition above Mmax-1 if a rate_split_bg_frac is active
bg_frac = getattr(self, 'rate_split_bg_frac', None)
if bg_frac is not None:
bins = self.mfd.get_annual_occurrence_rates()
threshold = mmax - 1.0
occ = [
(r * bg_frac if m >= threshold - bin_width / 4 else r)
for m, r in bins
]
self.mfd = EvenlyDiscretizedMFD(
min_mag=bins[0][0], bin_width=bin_width, occurrence_rates=occ)
[docs] def iter_ruptures(self, **kwargs):
"""
See :meth:
`openquake.hazardlib.source.base.BaseSeismicSource.iter_ruptures`
for description of parameters and return value.
Area sources are treated as a collection of point sources
(see :mod:`openquake.hazardlib.source.point`) with uniform parameters.
Ruptures of area source are just a union of ruptures
of those point sources. The actual positions of the implied
point sources form a uniformly spaced mesh on the polygon.
Polygon's method :meth:
`~openquake.hazardlib.geo.polygon.Polygon.discretize`
is used for creating a mesh of points on the source's area.
Constructor's parameter ``area_discretization`` is used as
polygon's discretization spacing (not to be confused with
rupture surface's mesh spacing which is as well provided
to the constructor).
The ruptures' occurrence rates are rescaled with respect to number
of points the polygon discretizes to.
"""
for src in self:
yield from src.iter_ruptures(**kwargs)
[docs] def count_ruptures(self):
"""
See
:meth:`openquake.hazardlib.source.base.BaseSeismicSource.count_ruptures`
for description of parameters and return value.
"""
polygon_mesh = self.polygon.discretize(self.area_discretization)
return (len(polygon_mesh) *
len(self.get_annual_occurrence_rates()) *
len(self.nodal_plane_distribution.data) *
len(self.hypocenter_distribution.data))
def __iter__(self):
"""
Split an area source into a generator of point sources.
MFDs will be rescaled appropriately for the number of points in the
area mesh.
"""
mesh = self.polygon.discretize(self.area_discretization)
num_points = len(mesh)
area_mfd = self.mfd
if isinstance(area_mfd, mfd.TruncatedGRMFD):
new_mfd = mfd.TruncatedGRMFD(
a_val=area_mfd.a_val - math.log10(num_points),
b_val=area_mfd.b_val,
bin_width=area_mfd.bin_width,
min_mag=area_mfd.min_mag,
max_mag=area_mfd.max_mag)
elif isinstance(area_mfd, mfd.EvenlyDiscretizedMFD):
new_occur_rates = [x/num_points for x in area_mfd.occurrence_rates]
new_mfd = mfd.EvenlyDiscretizedMFD(
min_mag=area_mfd.min_mag,
bin_width=area_mfd.bin_width,
occurrence_rates=new_occur_rates)
elif isinstance(area_mfd, mfd.ArbitraryMFD):
new_occur_rates = [x/num_points for x in area_mfd.occurrence_rates]
new_mfd = mfd.ArbitraryMFD(
magnitudes=area_mfd.magnitudes,
occurrence_rates=new_occur_rates)
elif isinstance(area_mfd, mfd.YoungsCoppersmith1985MFD):
new_mfd = mfd.YoungsCoppersmith1985MFD.from_characteristic_rate(
area_mfd.min_mag, area_mfd.b_val, area_mfd.char_mag,
area_mfd.char_rate / num_points, area_mfd.bin_width)
elif isinstance(area_mfd, mfd.AlternativeCharacteristicMFD):
new_mfd = mfd.AlternativeCharacteristicMFD(
min_mag=area_mfd.min_mag,
max_mag=area_mfd.max_mag,
bin_width=area_mfd.bin_width,
b_GR=area_mfd.b_GR,
b_AC=area_mfd.b_AC,
gamma=area_mfd.gamma,
delta_m_AC=area_mfd.delta_m_AC,
total_rate=area_mfd.total_rate / num_points)
else:
raise TypeError('Unknown MFD: %s' % area_mfd)
for i, (lon, lat) in enumerate(zip(mesh.lons, mesh.lats)):
pt = PointSource(
# Generate a new ID and name
source_id='%s:%s' % (self.source_id, i),
name=self.name,
tectonic_region_type=self.tectonic_region_type,
mfd=new_mfd,
rupture_mesh_spacing=self.rupture_mesh_spacing,
magnitude_scaling_relationship=
self.magnitude_scaling_relationship,
rupture_aspect_ratio=self.rupture_aspect_ratio,
upper_seismogenic_depth=self.upper_seismogenic_depth,
lower_seismogenic_depth=self.lower_seismogenic_depth,
location=geo.Point(lon, lat),
nodal_plane_distribution=self.nodal_plane_distribution,
hypocenter_distribution=self.hypocenter_distribution,
temporal_occurrence_model=self.temporal_occurrence_model)
pt._num_ruptures = pt.count_ruptures()
yield pt
[docs] def wkt(self):
"""
:returns: the geometry as a WKT string
"""
return self.polygon.wkt
