# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4
#
# Copyright (C) 2012-2023 GEM Foundation
#
# OpenQuake 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.
#
# OpenQuake 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 OpenQuake. If not, see <http://www.gnu.org/licenses/>.
"""
Module :mod:`openquake.hazardlib.geo.surface.simple_fault` defines
:class:`SimpleFaultSurface`.
"""
import math
import numpy
from openquake.baselib.node import Node
from openquake.hazardlib.geo.surface.base import BaseSurface
from openquake.hazardlib.geo.mesh import Mesh, RectangularMesh
from openquake.hazardlib.geo import utils as geo_utils
from openquake.hazardlib.geo.point import Point
from openquake.hazardlib.near_fault import get_plane_equation
[docs]def simple_fault_node(fault_trace, dip, upper_depth, lower_depth):
"""
:param fault_trace: an object with an attribute .points
:param dip: dip parameter
:param upper_depth: upper seismogenic depth
:param lower_depth: lower seismogenic depth
:returns: a Node of kind simpleFaultGeometry
"""
node = Node('simpleFaultGeometry')
line = []
for p in fault_trace.points:
line.append(p.longitude)
line.append(p.latitude)
node.append(Node('gml:LineString', nodes=[Node('gml:posList', {}, line)]))
node.append(Node('dip', {}, dip))
node.append(Node('upperSeismoDepth', {}, upper_depth))
node.append(Node('lowerSeismoDepth', {}, lower_depth))
return node
[docs]class SimpleFaultSurface(BaseSurface):
"""
Represent a fault surface as regular (uniformly spaced) 3D mesh of points.
:param mesh:
Instance of :class:`~openquake.hazardlib.geo.mesh.RectangularMesh`
representing surface geometry.
Another way to construct the surface object is to call
:meth:`from_fault_data`.
"""
def __init__(self, mesh):
self.mesh = mesh
assert 1 not in self.mesh.shape, (
"Mesh must have at least 2 nodes along both length and width.")
self.strike = self.dip = None
[docs] def get_dip(self):
"""
Return the fault dip as the average dip over the fault surface mesh.
The average dip is defined as the weighted mean inclination of top
row of mesh cells. See
:meth:`openquake.hazardlib.geo.mesh.RectangularMesh.get_mean_inclination_and_azimuth`
:returns:
The average dip, in decimal degrees.
"""
if self.dip is None:
# calculate weighted average dip and strike of only the top row
# of cells since those values are uniform along dip for simple
# faults
top_row = self.mesh[0:2]
self.dip, self.strike = top_row.get_mean_inclination_and_azimuth()
return self.dip
[docs] def get_strike(self):
"""
Return the fault strike as the average strike along the fault trace.
The average strike is defined as the weighted mean azimuth of top
row of mesh cells. See
:meth:`openquake.hazardlib.geo.mesh.RectangularMesh.get_mean_inclination_and_azimuth`
:returns:
The average strike, in decimal degrees.
"""
if self.strike is None:
self.get_dip() # this should cache strike value
return self.strike
[docs] @classmethod
def check_fault_data(cls, fault_trace, upper_seismogenic_depth,
lower_seismogenic_depth, dip, mesh_spacing):
"""
Verify the fault data and raise ``ValueError`` if anything is wrong.
This method doesn't have to be called by hands before creating the
surface object, because it is called from :meth:`from_fault_data`.
"""
if not len(fault_trace) >= 2:
raise ValueError("the fault trace must have at least two points")
if not fault_trace.horizontal():
raise ValueError("the fault trace must be horizontal")
tlats = [point.latitude for point in fault_trace.points]
tlons = [point.longitude for point in fault_trace.points]
if geo_utils.line_intersects_itself(tlons, tlats):
raise ValueError("fault trace intersects itself")
if not 0.0 < dip <= 90.0:
raise ValueError("dip must be between 0.0 and 90.0")
if not lower_seismogenic_depth > upper_seismogenic_depth:
raise ValueError("lower seismogenic depth must be greater than "
"upper seismogenic depth")
if not upper_seismogenic_depth >= fault_trace[0].depth:
raise ValueError("upper seismogenic depth must be greater than "
"or equal to depth of fault trace")
if not mesh_spacing > 0.0:
raise ValueError("mesh spacing must be positive")
[docs] @classmethod
def from_fault_data(cls, fault_trace, upper_seismogenic_depth,
lower_seismogenic_depth, dip, mesh_spacing):
"""
Create and return a fault surface using fault source data.
:param openquake.hazardlib.geo.line.Line fault_trace:
Geographical line representing the intersection between the fault
surface and the earth surface. The line must be horizontal (i.e.
all depth values must be equal). If the depths are not given, they
are assumed to be zero, meaning the trace intersects the surface at
sea level, e.g. fault_trace = Line([Point(1, 1), Point(1, 2)]).
:param upper_seismo_depth:
Minimum depth ruptures can reach, in km (i.e. depth
to fault's top edge).
:param lower_seismo_depth:
Maximum depth ruptures can reach, in km (i.e. depth
to fault's bottom edge).
:param dip:
Dip angle (i.e. angle between fault surface
and earth surface), in degrees.
:param mesh_spacing:
Distance between two subsequent points in a mesh, in km.
:returns:
An instance of :class:`SimpleFaultSurface` created using that data.
Uses :meth:`check_fault_data` for checking parameters.
"""
cls.check_fault_data(fault_trace, upper_seismogenic_depth,
lower_seismogenic_depth, dip, mesh_spacing)
# Loops over points in the top edge, for each point
# on the top edge compute corresponding point on the bottom edge, then
# computes equally spaced points between top and bottom points.
vdist_top = upper_seismogenic_depth - fault_trace[0].depth
vdist_bottom = lower_seismogenic_depth - fault_trace[0].depth
hdist_top = vdist_top / math.tan(math.radians(dip))
hdist_bottom = vdist_bottom / math.tan(math.radians(dip))
strike = fault_trace[0].azimuth(fault_trace[-1])
azimuth = (strike + 90.0) % 360
mesh = []
for point in fault_trace.resample(mesh_spacing):
top = point.point_at(hdist_top, vdist_top, azimuth)
bottom = point.point_at(hdist_bottom, vdist_bottom, azimuth)
mesh.append(top.equally_spaced_points(bottom, mesh_spacing))
# number of rows corresponds to number of points along dip
# number of columns corresponds to number of points along strike
surface_points = numpy.array(mesh).transpose().tolist()
mesh = RectangularMesh.from_points_list(surface_points)
assert 1 not in mesh.shape, (
"Mesh must have at least 2 nodes along both length and width."
" Possible cause: Mesh spacing could be too large with respect to"
" the fault length and width."
)
self = cls(mesh)
self.surface_nodes = [simple_fault_node(
fault_trace, dip,
upper_seismogenic_depth, lower_seismogenic_depth)]
return self
[docs] @classmethod
def get_fault_patch_vertices(cls, rupture_top_edge,
upper_seismogenic_depth,
lower_seismogenic_depth, dip, index_patch=1):
"""
Get surface main vertices.
Parameters are the same as for :meth:`from_fault_data`, excluding
fault_trace, and mesh spacing.
:param rupture_top_edge:
A instances of :class:`openquake.hazardlib.geo.line.Line`
representing the rupture surface's top edge.
:param index_patch:
Indicate the patch of the fault in order to output the vertices.
The fault patch numbering follows the same logic of the right-hand
rule i.e. patch with index 1 is the first patch along the trace.
:returns:
Four :class:~openquake.hazardlib.geo.point.Point objects
representing the four vertices of the target patch.
"""
# Similar to :meth:`from_fault_data`, we just don't resample edges
dip_tan = math.tan(math.radians(dip))
hdist_bottom = (
lower_seismogenic_depth - upper_seismogenic_depth) / dip_tan
strike = rupture_top_edge[0].azimuth(rupture_top_edge[-1])
azimuth = (strike + 90.0) % 360
# Collect coordinates of vertices on the top and bottom edge
lons = []
lats = []
deps = []
t_lon = []
t_lat = []
t_dep = []
for point in rupture_top_edge.points:
top_edge_point = point
bottom_edge_point = point.point_at(hdist_bottom, 0, azimuth)
lons.append(top_edge_point.longitude)
lats.append(top_edge_point.latitude)
deps.append(upper_seismogenic_depth)
t_lon.append(bottom_edge_point.longitude)
t_lat.append(bottom_edge_point.latitude)
t_dep.append(lower_seismogenic_depth)
all_lons = numpy.array(lons + list(reversed(t_lon)), float)
all_lats = numpy.array(lats + list(reversed(t_lat)), float)
all_deps = numpy.array(deps + list(reversed(t_dep)), float)
index1 = int(index_patch - 1)
index2 = int(index_patch)
index3 = int(2 * len(rupture_top_edge) - (index_patch + 1))
index4 = int(2 * len(rupture_top_edge) - index_patch)
p0 = Point(all_lons[index1], all_lats[index1], all_deps[index1])
p1 = Point(all_lons[index2], all_lats[index2], all_deps[index2])
p2 = Point(all_lons[index3], all_lats[index3], all_deps[index3])
p3 = Point(all_lons[index4], all_lats[index4], all_deps[index4])
return p0, p1, p2, p3
[docs] @classmethod
def hypocentre_patch_index(cls, hypocentre, rupture_top_edge,
upper_seismogenic_depth,
lower_seismogenic_depth, dip):
"""
This methods finds the index of the fault patch including
the hypocentre.
:param hypocentre:
:class:`~openquake.hazardlib.geo.point.Point` object
representing the location of hypocentre.
:param rupture_top_edge:
A instances of :class:`openquake.hazardlib.geo.line.Line`
representing the rupture surface's top edge.
:param upper_seismo_depth:
Minimum depth ruptures can reach, in km (i.e. depth
to fault's top edge).
:param lower_seismo_depth:
Maximum depth ruptures can reach, in km (i.e. depth
to fault's bottom edge).
:param dip:
Dip angle (i.e. angle between fault surface
and earth surface), in degrees.
:return:
An integer corresponding to the index of the fault patch which
contains the hypocentre.
"""
totaln_patch = len(rupture_top_edge)
indexlist = []
dist_list = []
for i, index in enumerate(range(1, totaln_patch)):
p0, p1, p2, p3 = cls.get_fault_patch_vertices(
rupture_top_edge, upper_seismogenic_depth,
lower_seismogenic_depth, dip, index_patch=index)
[normal, dist_to_plane] = get_plane_equation(p0, p1, p2,
hypocentre)
indexlist.append(index)
dist_list.append(dist_to_plane)
if numpy.allclose(dist_to_plane, 0., atol=25., rtol=0.):
return index
break
index = indexlist[numpy.argmin(dist_list)]
return index
[docs] @classmethod
def get_surface_vertexes(cls, fault_trace,
upper_seismogenic_depth,
lower_seismogenic_depth, dip):
"""
Get surface main vertexes.
Parameters are the same as for :meth:`from_fault_data`, excluding
mesh spacing.
:returns:
Instance of :class:`~openquake.hazardlib.geo.polygon.Polygon`
describing the surface projection of the simple fault with
specified parameters.
"""
# Similar to :meth:`from_fault_data`, we just don't resample edges
dip_tan = math.tan(math.radians(dip))
hdist_top = upper_seismogenic_depth / dip_tan
hdist_bottom = lower_seismogenic_depth / dip_tan
strike = fault_trace[0].azimuth(fault_trace[-1])
azimuth = (strike + 90.0) % 360
# Collect coordinates of vertices on the top and bottom edge
lons = []
lats = []
for point in fault_trace.points:
top_edge_point = point.point_at(hdist_top, 0, azimuth)
bottom_edge_point = point.point_at(hdist_bottom, 0, azimuth)
lons.append(top_edge_point.longitude)
lats.append(top_edge_point.latitude)
lons.append(bottom_edge_point.longitude)
lats.append(bottom_edge_point.latitude)
lons = numpy.array(lons, float)
lats = numpy.array(lats, float)
return lons, lats
[docs] @classmethod
def surface_projection_from_fault_data(cls, fault_trace,
upper_seismogenic_depth,
lower_seismogenic_depth, dip):
"""
Get a surface projection of the simple fault surface.
Parameters are the same as for :meth:`from_fault_data`, excluding
mesh spacing.
:returns:
Instance of :class:`~openquake.hazardlib.geo.polygon.Polygon`
describing the surface projection of the simple fault with
specified parameters.
"""
lons, lats = cls.get_surface_vertexes(fault_trace,
upper_seismogenic_depth,
lower_seismogenic_depth, dip)
return Mesh(lons, lats, depths=None).get_convex_hull()
[docs] def get_width(self):
"""
Return surface's width (that is surface extension along the
dip direction) in km.
The width is computed as the average width along the surface.
See
:meth:`openquake.hazardlib.geo.mesh.RectangularMesh.get_mean_width`
"""
# calculate width only along the first mesh column, because
# width is uniform for simple faults
left_column = self.mesh[:, 0:2]
return left_column.get_mean_width()