Source code for openquake.hazardlib.geo.surface.simple_fault

# -*- 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()