# The Hazard Library
# Copyright (C) 2012-2014, 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.geo.surface.simple_fault` defines
:class:`SimpleFaultSurface`.
"""
import math
import numpy
from openquake.hazardlib.geo.surface.base import BaseQuadrilateralSurface
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]class SimpleFaultSurface(BaseQuadrilateralSurface):
"""
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):
super(SimpleFaultSurface, self).__init__()
self.mesh = mesh
assert not 1 in self.mesh.shape, (
"Mesh must have at least 2 nodes along both length and width."
)
self.strike = self.dip = None
[docs] def _create_mesh(self):
"""
Return a mesh provided to object's constructor.
"""
return self.mesh
[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.get_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
@classmethod
[docs] 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.on_surface():
raise ValueError("the fault trace must be defined on the surface")
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 >= 0.0:
raise ValueError("upper seismo depth must be non-negative")
if not mesh_spacing > 0.0:
raise ValueError("mesh spacing must be positive")
@classmethod
[docs] 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.
: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
vdist_bottom = lower_seismogenic_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."
)
return cls(mesh)
@classmethod
[docs] 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
@classmethod
[docs] 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)
for index in 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)
if (numpy.allclose(dist_to_plane, 0., atol=25., rtol=0.)):
return index
@classmethod
[docs] 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
@classmethod
[docs] 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.get_mesh()[:, 0:2]
return left_column.get_mean_width()