# -*- 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 faultPSHA calculation."""importmathimportnumpyasnpfromopenquake.hazardlib.geoimportgeodeticasgeodimportscipy.spatial.distanceasdst
[docs]defget_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)returnnp.array([dists*math.sin(math.radians(azims)),dists*math.cos(math.radians(azims)),depths])
[docs]defget_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-p0v=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 equationdist_to_plane=np.dot(p0,normal)returnnormal,dist_to_plane
[docs]defprojection_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]])returnpp
[docs]defvectors2angle(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))returnnp.arctan2(sinang,cosang)
[docs]defaverage_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 planesite_xyz=get_xyz_from_ll(site,reference)zs=dst.pdist([pp,site_xyz])ifsite_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 hypocentrehyp_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.5r_hyp=(l2**2+zs**2)**0.5p0_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_xyzphi=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 valuefs=(ix**2+inn**2+iphi**2)**0.5/ereturnfs,rd,r_hyp
[docs]defisochone_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 """ife==0.:c_prime=0.8elife>0.:c_prime=1./((1./0.8)-((r_hyp-rd)/e))returnc_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-pani=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.5vector2=(seg2_start-seg2_end)/ \
sum((seg2_start-seg2_end)**2)**0.5vector3=(seg1_end-seg1_start)/ \
sum((seg1_end-seg1_start)**2)**0.5vector4=(p_intersect.flatten()-seg1_start)/ \
sum((p_intersect.flatten()-seg1_start)**2)**0.5returnp_intersect,vector1,vector2,vector3,vector4
[docs]defdirectp(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 coordinatenode0_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=ppe=[]# Loop each segments on the patch to find Pdsegment_s=[node0_xyz,node1_xyz,node2_xyz,node3_xyz]segment_e=[node1_xyz,node2_xyz,node3_xyz,node0_xyz]# set buffering bubuf=0.0001atol=0.0001loop=Trueexit_flag=Falselooptime=0.whileloop:x_min=np.min(np.array([node0_xyz[0],node1_xyz[0],node2_xyz[0],node3_xyz[0]]))-bufx_max=np.max(np.array([node0_xyz[0],node1_xyz[0],node2_xyz[0],node3_xyz[0]]))+bufy_min=np.min(np.array([node0_xyz[1],node1_xyz[1],node2_xyz[1],node3_xyz[1]]))-bufy_max=np.max(np.array([node0_xyz[1],node1_xyz[1],node2_xyz[1],node3_xyz[1]]))+bufn_seg=0exit_flag=Falseforseg_s,seg_einzip(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=pdphpd=p_intersectexit_flag=Truebreak# when the pp located within the fault rupture plane, e = ppphifnote:if(pp_xyz[0]>=x_min)&(pp_xyz[0]<=x_max):if(pp_xyz[1]>=y_min)&(pp_xyz[1]<=y_max):pd=pp_xyze=ppphexit_flag=Trueifexit_flag:breakifnote:looptime+=1atol=0.0001*looptimebuf=0.0001*looptime# if pd is located at 2nd fault segment, then the DPP calculation will# keep going on the next fault patchifn_seg==2:go_next_patch=Trueelse:go_next_patch=Falsereturnpd,e,go_next_patch