Source code for openquake.hazardlib.near_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.nearfault` provides methods for near fault
PSHA calculation.
"""

import math
import numpy as np
from openquake.hazardlib.geo import geodetic as geod
import scipy.spatial.distance as dst


[docs]def get_xyz_from_ll(projected, reference):
    """
    This method computes the x, y and z coordinates of a set of points
    provided a reference point.

    :param projected:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the coordinates of target point to be projected
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the coordinates of the reference point.
    :returns: a 3D vector
    """

    azims = geod.azimuth(reference.longitude, reference.latitude,
                         projected.longitude, projected.latitude)
    depths = np.subtract(reference.depth, projected.depth)
    dists = geod.geodetic_distance(reference.longitude,
                                   reference.latitude,
                                   projected.longitude,
                                   projected.latitude)
    return np.array([dists * math.sin(math.radians(azims)),
                     dists * math.cos(math.radians(azims)),
                     depths])


[docs]def get_plane_equation(p0, p1, p2, reference):
    '''
    Define the equation of target fault plane passing through 3 given points
    which includes two points on the fault trace and one point on the
    fault plane but away from the fault trace. Note: in order to remain the
    consistency of the fault normal vector direction definition, the order
    of the three given points is strickly defined.

    :param p0:
        The fault trace and is the closer points from the starting point of
        fault trace.
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the one vertex of the fault patch.
    :param p1:
        The fault trace and is the further points from the starting point of
        fault trace.
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the one vertex of the fault patch.
    :param p2:
        The point on the fault plane but away from the fault trace.
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the one vertex of the fault patch.
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the origin of the cartesian system used the represent
        objects in a projected reference
    :returns:
        normal: normal vector of the plane (a,b,c)
        dist_to_plane: d in the plane equation, ax + by + cz = d
    '''

    p0_xyz = get_xyz_from_ll(p0, reference)
    p1_xyz = get_xyz_from_ll(p1, reference)
    p2_xyz = get_xyz_from_ll(p2, reference)

    p0 = np.array(p0_xyz)
    p1 = np.array(p1_xyz)
    p2 = np.array(p2_xyz)

    u = p1 - p0
    v = p2 - p0

    # vector normal to plane, ax+by+cy = d, normal=(a,b,c)
    normal = np.cross(u, v)

    # Define the d for the plane equation
    dist_to_plane = np.dot(p0, normal)

    return normal, dist_to_plane


[docs]def projection_pp(site, normal, dist_to_plane, reference):
    '''
    This method finds the projection of the site onto the plane containing
    the slipped area, defined as the Pp(i.e. 'perpendicular projection of
    site location onto the fault plane' Spudich et al. (2013) - page 88)
    given a site.

    :param site:
        Location of the site, [lon, lat, dep]
    :param normal:
        Normal to the plane including the fault patch,
        describe by a normal vector[a, b, c]
    :param dist_to_plane:
        D in the plane equation,  ax + by + cz = d
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of project reference point
    :returns:
        pp, the projection point, [ppx, ppy, ppz], in xyz domain
        , a numpy array.
    '''

    # Transform to xyz coordinate
    [site_x, site_y, site_z] = get_xyz_from_ll(site, reference)

    a = np.array([(1, 0, 0, -normal[0]),
                  (0, 1, 0, -normal[1]),
                  (0, 0, 1, -normal[2]),
                  (normal[0], normal[1], normal[2], 0)])
    b = np.array([site_x, site_y, site_z, dist_to_plane])

    x = np.linalg.solve(a, b)
    pp = np.array([x[0], x[1], x[2]])
    return pp


[docs]def vectors2angle(v1, v2):
    """ Returns the angle in radians between vectors 'v1' and 'v2'.

    :param v1:
        vector, a numpy array
    :param v2:
        vector, a numpy array
    :returns:
        the angle in radians between the two vetors
    """
    cosang = np.dot(v1, v2)
    sinang = np.linalg.norm(np.cross(v1, v2))
    return np.arctan2(sinang, cosang)


[docs]def average_s_rad(site, hypocenter, reference, pp,
                  normal, dist_to_plane, e, p0, p1, delta_slip):
    """
    Gets the average S-wave radiation pattern given an e-path as described in:
    Spudich et al. (2013) "Final report of the NGA-West2 directivity working
    group", PEER report, page 90- 92 and computes: the site to the direct point
    distance, rd, and the hypocentral distance, r_hyp.

    :param site:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the target site
    :param hypocenter:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of hypocenter
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location
        of the reference point for coordinate projection within the
        calculation. The suggested reference point is Epicentre.
    :param pp:
        the projection point pp on the patch plane,
        a numpy array
    :param normal:
        normal of the plane, describe by a normal vector[a, b, c]
    :param dist_to_plane:
        d is the constant term in the plane equation, e.g., ax + by + cz = d
    :param e:
        a float defining the E-path length, which is the distance from
        Pd(direction) point to hypocentre. In km.
    :param p0:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the starting point on fault segment
    :param p1:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of the ending point on fault segment.
    :param delta_slip:
        slip direction away from the strike direction, in decimal degrees.
        A positive angle is generated by a counter-clockwise rotation.
    :returns:
        fs, float value of the average S-wave radiation pattern.
        rd, float value of the distance from site to the direct point.
        r_hyp, float value of the hypocetre distance.
    """
    # Obtain the distance of Ps and Pp. If Ps is above the fault plane
    # zs is positive, and negative when Ps is below the fault plane
    site_xyz = get_xyz_from_ll(site, reference)
    zs = dst.pdist([pp, site_xyz])

    if site_xyz[0] * normal[0] + site_xyz[1] * normal[1] + site_xyz[2] * \
       normal[2] - dist_to_plane > 0:
        zs = -zs

    # Obtain the distance of Pp and hypocentre
    hyp_xyz = get_xyz_from_ll(hypocenter, reference)
    hyp_xyz = np.array(hyp_xyz).reshape(1, 3).flatten()
    l2 = dst.pdist([pp, hyp_xyz])

    rd = ((l2 - e) ** 2 + zs ** 2) ** 0.5
    r_hyp = (l2 ** 2 + zs ** 2) ** 0.5
    p0_xyz = get_xyz_from_ll(p0, reference)
    p1_xyz = get_xyz_from_ll(p1, reference)
    u = (np.array(p1_xyz) - np.array(p0_xyz))
    v = pp - hyp_xyz
    phi = vectors2angle(u, v) - np.deg2rad(delta_slip)
    ix = np.cos(phi) * (2 * zs * (l2 / r_hyp - (l2 - e) / rd) -
                        zs * np.log((l2 + r_hyp) / (l2 - e + rd)))
    inn = np.cos(phi) * (-2 * zs ** 2 * (1 / r_hyp - 1 / rd)
                         - (r_hyp - rd))
    iphi = np.sin(phi) * (zs * np.log((l2 + r_hyp) / (l2 - e + rd)))

    # Obtain the final average radiation pattern value
    fs = (ix ** 2 + inn ** 2 + iphi ** 2) ** 0.5 / e

    return fs, rd, r_hyp


[docs]def isochone_ratio(e, rd, r_hyp):
    """
    Get the isochone ratio as described in Spudich et al. (2013) PEER
    report, page 88.

    :param e:
        a float defining the E-path length, which is the distance from
        Pd(direction) point to hypocentre. In km.
    :param rd:
        float, distance from the site to the direct point.
    :param r_hyp:
        float, the hypocentre distance.
    :returns:
        c_prime, a float defining the isochone ratio
    """

    if e == 0.:
        c_prime = 0.8
    elif e > 0.:
        c_prime = 1. / ((1. / 0.8) - ((r_hyp - rd) / e))

    return c_prime


def _intersection(seg1_start, seg1_end, seg2_start, seg2_end):
    """
    Get the intersection point between two segments. The calculation is in
    Catestian coordinate system.

    :param seg1_start:
        A numpy array,
        representing one end point of a segment(e.g. segment1)
        segment.
    :param seg1_end:
        A numpy array,
        representing the other end point of the first segment(e.g. segment1)
    :param seg2_start:
        A numpy array,
        representing one end point of the other segment(e.g. segment2)
        segment.
    :param seg2_end:
        A numpy array,
        representing the other end point of the second segment(e.g. segment2)
    :returns:
        p_intersect, :a numpy ndarray.
        representing the location of intersection point of the two
        given segments
        vector1, a numpy array, vector defined by intersection point and
        seg2_end
        vector2, a numpy array, vector defined by seg2_start and seg2_end
        vector3, a numpy array, vector defined by seg1_start and seg1_end
        vector4, a numpy array, vector defined by intersection point
        and seg1_start
    """

    pa = np.array([seg1_start, seg2_start])
    pb = np.array([seg1_end, seg2_end])

    si = pb - pa

    ni = si / np.power(
        np.dot(np.sum(si ** 2, axis=1).reshape(2, 1),
               np.ones((1, 3))), 0.5)

    nx = ni[:, 0].reshape(2, 1)
    ny = ni[:, 1].reshape(2, 1)
    nz = ni[:, 2].reshape(2, 1)
    sxx = np.sum(nx ** 2 - 1)
    syy = np.sum(ny ** 2 - 1)
    szz = np.sum(nz ** 2 - 1)
    sxy = np.sum(nx * ny)
    sxz = np.sum(nx * nz)
    syz = np.sum(ny * nz)
    s = np.array([sxx, sxy, sxz, sxy, syy, syz, sxz, syz,
                 szz]).reshape(3, 3)

    cx = np.sum(pa[:, 0].reshape(2, 1) * (nx ** 2 - 1) +
                pa[:, 1].reshape(2, 1) * [nx * ny] +
                pa[:, 2].reshape(2, 1) * (nx * nz))
    cy = np.sum(pa[:, 0].reshape(2, 1) * [nx * ny] +
                pa[:, 1].reshape(2, 1) * [ny ** 2 - 1] +
                pa[:, 2].reshape(2, 1) * [ny * nz])
    cz = np.sum(pa[:, 0].reshape(2, 1) * [nx * nz] +
                pa[:, 1].reshape(2, 1) * [ny * nz] +
                pa[:, 2].reshape(2, 1) * [nz ** 2 - 1])
    c = np.array([cx, cy, cz]).reshape(3, 1)
    p_intersect = np.linalg.solve(s, c)

    vector1 = (p_intersect.flatten() - seg2_end) / \
        sum((p_intersect.flatten() - seg2_end) ** 2) ** 0.5
    vector2 = (seg2_start - seg2_end) / \
        sum((seg2_start - seg2_end) ** 2) ** 0.5
    vector3 = (seg1_end - seg1_start) / \
        sum((seg1_end - seg1_start) ** 2) ** 0.5
    vector4 = (p_intersect.flatten() - seg1_start) / \
        sum((p_intersect.flatten() - seg1_start) ** 2) ** 0.5

    return p_intersect, vector1, vector2, vector3, vector4


[docs]def directp(node0, node1, node2, node3, hypocenter, reference, pp):
    """
    Get the Direct Point and the corresponding E-path as described in
    Spudich et al. (2013). This method also provides a logical variable
    stating if the DPP calculation must consider the neighbouring patch.
    To define the intersection point(Pd) of PpPh line segment and fault plane,
    we obtain the intersection points(Pd) with each side of fault plan, and
    check which intersection point(Pd) is the one fitting the definition in
    the Chiou and Spudich(2014) directivity model.
    Two possible locations for Pd, the first case, Pd locates on the side of
    the fault patch when Pp is not inside the fault patch. The second case is
    when Pp is inside the fault patch, then Pd=Pp.

    For the first case, it follows three conditions:
    1. the PpPh and PdPh line vector are the same,
    2. PpPh >= PdPh,
    3. Pd is not inside the fault patch.

    If we can not find solution for all the four possible intersection points
    for the first case, we check if the intersection point fit the second case
    by checking if Pp is inside the fault patch.

    Because of the coordinate system mapping(from geographic system to
    Catestian system), we allow an error when we check the location. The allow
    error will keep increasing after each loop when no solution in the two
    cases are found, until the solution get obtained.

    :param node0:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment.
    :param node1:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node2:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param node3:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of one vertices on the target fault
        segment. Note, the order should be clockwise.
    :param hypocenter:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of floating hypocenter on each segment
        calculation. In the method, we take the direction point of the
        previous fault patch as hypocentre for the current fault patch.
    :param reference:
        :class:`~openquake.hazardlib.geo.point.Point` object
        representing the location of reference point for projection
    :param pp:
        the projection of the site onto the plane containing the fault
        slipped area. A numpy array.
    :returns:
        Pd, a numpy array, representing the location of direction point
        E, the distance from direction point to hypocentre.
        go_next_patch, flag indicates if the calculation goes on the next
        fault patch. 1: yes, 0: no.
    """
    # Find the intersection point Pd, by checking if the PdPh share the
    # same vector with PpPh,  and PpPh >= PdPh
    # Transform to xyz coordinate
    node0_xyz = get_xyz_from_ll(node0, reference)
    node1_xyz = get_xyz_from_ll(node1, reference)
    node2_xyz = get_xyz_from_ll(node2, reference)
    node3_xyz = get_xyz_from_ll(node3, reference)
    hypocenter_xyz = get_xyz_from_ll(hypocenter, reference)
    hypocenter_xyz = np.array(hypocenter_xyz).flatten()
    pp_xyz = pp
    e = []

    # Loop each segments on the patch to find Pd
    segment_s = [node0_xyz, node1_xyz, node2_xyz, node3_xyz]
    segment_e = [node1_xyz, node2_xyz, node3_xyz, node0_xyz]

    # set buffering bu
    buf = 0.0001
    atol = 0.0001

    loop = True
    exit_flag = False
    looptime = 0.
    while loop:
        x_min = np.min(np.array([node0_xyz[0], node1_xyz[0], node2_xyz[0],
                                node3_xyz[0]])) - buf
        x_max = np.max(np.array([node0_xyz[0], node1_xyz[0], node2_xyz[0],
                                node3_xyz[0]])) + buf
        y_min = np.min(np.array([node0_xyz[1], node1_xyz[1], node2_xyz[1],
                                node3_xyz[1]])) - buf
        y_max = np.max(np.array([node0_xyz[1], node1_xyz[1], node2_xyz[1],
                                node3_xyz[1]])) + buf
        n_seg = 0
        exit_flag = False
        for seg_s, seg_e in zip(segment_s, segment_e):
            seg_s = seg_s.flatten()
            seg_e = seg_e.flatten()
            p_intersect, vector1, vector2, vector3, vector4 = _intersection(
                seg_s, seg_e, pp_xyz, hypocenter_xyz)

            ppph = dst.pdist([pp, hypocenter_xyz])
            pdph = dst.pdist([p_intersect.flatten(), hypocenter_xyz])
            n_seg = n_seg + 1

            # Check that the direction of the hyp-pp and hyp-pd vectors
            # have are the same.
            if (np.allclose(vector1.flatten(), vector2,
                            atol=atol, rtol=0.)):
                if (np.allclose(vector3.flatten(), vector4, atol=atol,
                                rtol=0.)):

                    # Check if ppph >= pdph.
                    if (ppph >= pdph):
                        if (p_intersect[0] >= x_min) & (p_intersect[0] <=
                                                        x_max):
                            if (p_intersect[1] >= y_min) & (p_intersect[1]
                                                            <= y_max):
                                e = pdph
                                pd = p_intersect
                                exit_flag = True
                                break
            # when the pp located within the fault rupture plane, e = ppph
            if not e:
                if (pp_xyz[0] >= x_min) & (pp_xyz[0] <= x_max):
                    if (pp_xyz[1] >= y_min) & (pp_xyz[1] <= y_max):
                        pd = pp_xyz
                        e = ppph
                        exit_flag = True
        if exit_flag:
            break

        if not e:
            looptime += 1
            atol = 0.0001 * looptime
            buf = 0.0001 * looptime
    # if pd is located at 2nd fault segment, then the DPP calculation will
    # keep going on the next fault patch
    if n_seg == 2:
        go_next_patch = True
    else:
        go_next_patch = False

    return pd, e, go_next_patch