Source code for openquake.hazardlib.geo.surface.base
# -*- 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.base` implements:class:`BaseSurface` and :class:`BaseSurface`."""importnumpyimportmathfromopenquake.hazardlib.geo.meshimportMeshfromopenquake.hazardlib.geoimportgeodetic,utils,Point,Line,\
RectangularMeshdef_get_finite_mesh(mesh):ok=numpy.isfinite(mesh.lons.flat)ifnumpy.all(ok):returnmeshok=numpy.reshape(ok,mesh.lons.shape)returnMesh(mesh.lons[ok],mesh.lats[ok],mesh.depths[ok])def_find_turning_points(mesh,tol=1.0):""" Identifies the turning points in a rectangular mesh based on the deviation in the azimuth between successive points on the upper edge. A turning point is flagged if the change in azimuth change is greater than the specified tolerance (in degrees) :param mesh: Mesh for downsampling as instance of :class: openquake.hazardlib.geo.mesh.RectangularMesh :param float tol: Maximum difference in azimuth (decimal degrees) between successive points to identify a turning point :returns: Column indices of turning points (as numpy array) """assertisinstance(mesh,RectangularMesh)idx1=numpy.isfinite(mesh.lons[0,:-1])idx2=numpy.isfinite(mesh.lons[0,1:])idx=numpy.where(numpy.logical_and(idx1,idx2))[0]azimuths=geodetic.azimuth(mesh.lons[0,idx],mesh.lats[0,idx],mesh.lons[0,idx+1],mesh.lats[0,idx+1])naz=len(azimuths)azim=azimuths[0]# Retain initial pointidx=[0]# Add more pointsforiinrange(1,naz):dff=_angle_difference(azimuths[i],azim)ifdff>tol:idx.append(i)azim=azimuths[i]# Add on last point - if not already in the setifidx[-1]!=mesh.lons.shape[1]-1:idx.append(mesh.lons.shape[1]-1)returnnumpy.array(idx)def_angle_difference(a_a,a_b):""" Computes the absolute difference between angle `a_a` and `a_b` """dff=a_a-a_breturn(dff+540)%360-180
[docs]defdownsample_mesh(mesh,tol=1.0):""" Returns a mesh sampled at a lower resolution - if the difference in azimuth is larger than the specified tolerance a turn is assumed :returns: Downsampled mesh as instance of :class: openquake.hazardlib.geo.mesh.RectangularMesh """idx=_find_turning_points(mesh,tol)ifmesh.depthsisnotNone:returnRectangularMesh(lons=mesh.lons[:,idx],lats=mesh.lats[:,idx],depths=mesh.depths[:,idx])else:returnRectangularMesh(lons=mesh.lons[:,idx],lats=mesh.lats[:,idx])
[docs]defdownsample_trace(mesh,tol=1.0):""" Downsamples the upper edge of a fault within a rectangular mesh, retaining node points only if changes in direction on the order of tol are found :returns: Downsampled edge as a numpy array of [long, lat, depth] """idx=_find_turning_points(mesh,tol)ifmesh.depthsisnotNone:returnnumpy.column_stack([mesh.lons[0,idx],mesh.lats[0,idx],mesh.depths[0,idx]])else:returnnumpy.column_stack([mesh.lons[0,idx],mesh.lats[0,idx]])
def_get_p1_p2(clsname,top_edge,i):# returns two points used in get_rx_distanceif(clsname=='KiteSurface'andnumpy.isnan(top_edge.lons[0,i])ornumpy.isnan(top_edge.lons[0,i+1])):raiseValueError('Rx calculation has less than two points')p1=Point(top_edge.lons[0,i],top_edge.lats[0,i],top_edge.depths[0,i])p2=Point(top_edge.lons[0,i+1],top_edge.lats[0,i+1],top_edge.depths[0,i+1])returnp1,p2
[docs]classBaseSurface:""" Base class for a surface in 3D-space. """def__init__(self,mesh=None):self.mesh=meshself.idx=None
[docs]defget_min_distance(self,mesh):""" Compute and return the minimum distance from the surface to each point of ``mesh``. This distance is sometimes called ``Rrup``. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to calculate minimum distance to. :returns: A numpy array of distances in km. """fmesh=_get_finite_mesh(self.mesh)returnfmesh.get_min_distance(mesh)
[docs]defget_closest_points(self,mesh):""" For each point from ``mesh`` find a closest point belonging to surface. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to find closest points to. :returns: :class:`~openquake.hazardlib.geo.mesh.Mesh` of the same shape as ``mesh`` with closest surface's points on respective indices. """fmesh=_get_finite_mesh(self.mesh)returnfmesh.get_closest_points(mesh)
[docs]defget_joyner_boore_distance(self,mesh):""" Compute and return Joyner-Boore (also known as ``Rjb``) distance to each point of ``mesh``. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to calculate Joyner-Boore distance to. :returns: Numpy array of closest distances between the projections of surface and each point of the ``mesh`` to the earth surface. """fmesh=_get_finite_mesh(self.mesh)returnfmesh.get_joyner_boore_distance(mesh)
[docs]defget_ry0_distance(self,mesh):""" Compute the minimum distance between each point of a mesh and the great circle arcs perpendicular to the average strike direction of the fault trace and passing through the end-points of the trace. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to calculate Ry0-distance to. :returns: Numpy array of distances in km. """# This computes ry0 by using an average strike directiontop_edge=self.mesh[0:1]mean_strike=self.get_strike()# Manage the case of kite fault surfacesia=0ib=-1if(self.__class__.__name__=='KiteSurface'):idx=numpy.nonzero(numpy.isfinite(top_edge.lons[0,:]))[0]ia=min(idx)ib=max(idx)# Computing the distances between the sites and the two lines# perpendicular to the strike passing trough the two extremes# of the top of the rupturedst1=geodetic.distance_to_arc(top_edge.lons[0,ia],top_edge.lats[0,ia],(mean_strike+90.)%360,mesh.lons,mesh.lats)dst2=geodetic.distance_to_arc(top_edge.lons[0,ib],top_edge.lats[0,ib],(mean_strike+90.)%360,mesh.lons,mesh.lats)# Get the shortest distance from the two linesidx=numpy.sign(dst1)==numpy.sign(dst2)dst=numpy.zeros_like(dst1)dst[idx]=numpy.fmin(numpy.abs(dst1[idx]),numpy.abs(dst2[idx]))ifnumpy.any(numpy.isnan(dst)):raiseValueError('NaN in Ry0')returndst
[docs]defget_rx_distance(self,mesh):""" Compute distance between each point of mesh and surface's great circle arc. Distance is measured perpendicular to the rupture strike, from the surface projection of the updip edge of the rupture, with the down dip direction being positive (this distance is usually called ``Rx``). In other words, is the horizontal distance to top edge of rupture measured perpendicular to the strike. Values on the hanging wall are positive, values on the footwall are negative. :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points to calculate Rx-distance to. :returns: Numpy array of distances in km. """top_edge=self.mesh[0:1]dists=[]ia=0ib=top_edge.lons.shape[1]-2ifself.__class__.__name__=='KiteSurface':idxs=numpy.nonzero(numpy.isfinite(top_edge.lons[0,:]))[0]ia=min(idxs)ib=sorted(idxs)[-2]iftop_edge.lons.shape[1]<3:p1,p2=_get_p1_p2(self.__class__.__name__,top_edge,i=0)azimuth=p1.azimuth(p2)dists.append(geodetic.distance_to_arc(p1.longitude,p1.latitude,azimuth,mesh.lons,mesh.lats))else:foriinrange(top_edge.lons.shape[1]-1):try:p1,p2=_get_p1_p2(self.__class__.__name__,top_edge,i)exceptValueError:continue# Swappingifi==0:p1,p2=p2,p1# Computing azimuth and distanceifi==iaori==ib:azimuth=p1.azimuth(p2)dst=geodetic.distance_to_semi_arc(p1.longitude,p1.latitude,azimuth,mesh.lons,mesh.lats)else:dst=geodetic.min_distance_to_segment(numpy.array([p1.longitude,p2.longitude]),numpy.array([p1.latitude,p2.latitude]),mesh.lons,mesh.lats)# Correcting the sign of the distanceifi==0:dst*=-1dists.append(dst)# Computing distancesdists=numpy.array(dists)iii=numpy.abs(dists).argmin(axis=0)dst=dists[iii,numpy.arange(dists.shape[1])]ifnumpy.any(numpy.isnan(dst)):raiseValueError('NaN in Rx')returndst
[docs]defget_top_edge_depth(self):""" Return minimum depth of surface's top edge. :returns: Float value, the vertical distance between the earth surface and the shallowest point in surface's top edge in km. """top_edge=self.mesh[0:1]iftop_edge.depthsisNone:return0else:dep=numpy.min(top_edge.depths)returndep
[docs]defget_top_edge_centroid(self):""" Return :class:`~openquake.hazardlib.geo.point.Point` representing the surface's top edge centroid. """top_edge=self.mesh[0:1]returntop_edge.get_middle_point()
[docs]defget_area(self):""" Compute area as the sum of the mesh cells area values. """fromopenquake.hazardlib.geo.surface.kite_faultimportKiteSurfaceifisinstance(self,KiteSurface):_,_,_,area=self.get_cell_dimensions()else:mesh=self.mesh_,_,_,area=mesh.get_cell_dimensions()returnnumpy.sum(area)
[docs]defget_bounding_box(self):""" Compute surface bounding box from surface mesh representation. That is extract longitudes and latitudes of mesh points and calls: :meth:`openquake.hazardlib.geo.utils.get_spherical_bounding_box` :return: A tuple of four items. These items represent western, eastern, northern and southern borders of the bounding box respectively. Values are floats in decimal degrees. """mesh=_get_finite_mesh(self.mesh)returnutils.get_spherical_bounding_box(mesh.lons.flatten(),mesh.lats.flatten())
[docs]defget_middle_point(self):""" Compute coordinates of surface middle point. The actual definition of ``middle point`` depends on the type of surface geometry. :return: instance of :class:`openquake.hazardlib.geo.point.Point` representing surface middle point. """returnself.mesh.get_middle_point()
def_boundaries(self,name):# name is one of lons, lats, mesharr=getattr(self.mesh,name)returnnumpy.concatenate((arr[0,:],arr[1:,-1],arr[-1,:-1][::-1],arr[:-1,0][::-1]))
[docs]defget_surface_boundaries(self):""" Returns the boundaries in the same format as a multiplanar surface, with two lists of lons and lats """returnself._boundaries('lons'),self._boundaries('lats')
[docs]defget_surface_boundaries_3d(self):""" Returns the boundaries as three lists of lons, lats, depths """return(self._boundaries('lons'),self._boundaries('lats'),self._boundaries('depths'))
[docs]defget_resampled_top_edge(self,angle_var=3.0):""" This methods computes a simplified representation of a fault top edge by removing the points that are not describing a change of direction, provided a certain tolerance angle. :param float angle_var: Number representing the maximum deviation (in degrees) admitted without the creation of a new segment :returns: A :class:`~openquake.hazardlib.geo.line.Line` representing the rupture surface's top edge. """mesh=self.meshtop_edge=[Point(mesh.lons[0][0],mesh.lats[0][0],mesh.depths[0][0])]foriinrange(len(mesh.triangulate()[1][0])-1):v1=numpy.asarray(mesh.triangulate()[1][0][i])v2=numpy.asarray(mesh.triangulate()[1][0][i+1])cosang=numpy.dot(v1,v2)sinang=numpy.linalg.norm(numpy.cross(v1,v2))angle=math.degrees(numpy.arctan2(sinang,cosang))ifabs(angle)>angle_var:top_edge.append(Point(mesh.lons[0][i+1],mesh.lats[0][i+1],mesh.depths[0][i+1]))top_edge.append(Point(mesh.lons[0][-1],mesh.lats[0][-1],mesh.depths[0][-1]))line_top_edge=Line(top_edge)returnline_top_edge
[docs]defget_hypo_location(self,mesh_spacing,hypo_loc=None):""" The method determines the location of the hypocentre within the rupture :param mesh: :class:`~openquake.hazardlib.geo.mesh.Mesh` of points :param mesh_spacing: The desired distance between two adjacent points in source's ruptures' mesh, in km. Mainly this parameter allows to balance the trade-off between time needed to compute the distance between the rupture surface and a site and the precision of that computation. :param hypo_loc: Hypocentre location as fraction of rupture plane, as a tuple of (Along Strike, Down Dip), e.g. a hypocentre located in the centroid of the rupture would be input as (0.5, 0.5), whereas a hypocentre located in a position 3/4 along the length, and 1/4 of the way down dip of the rupture plane would be entered as (0.75, 0.25). :returns: Hypocentre location as instance of :class:`~openquake.hazardlib.geo.point.Point` """mesh=self.meshcentroid=mesh.get_middle_point()ifhypo_locisNone:returncentroidtotal_len_y=(len(mesh.depths)-1)*mesh_spacingy_distance=hypo_loc[1]*total_len_yy_node=int(numpy.round(y_distance/mesh_spacing))total_len_x=(len(mesh.lons[y_node])-1)*mesh_spacingx_distance=hypo_loc[0]*total_len_xx_node=int(numpy.round(x_distance/mesh_spacing))hypocentre=Point(mesh.lons[y_node][x_node],mesh.lats[y_node][x_node],mesh.depths[y_node][x_node])returnhypocentre
[docs]defget_azimuth(self,mesh):""" This method computes the azimuth of a set of points in a :class:`openquake.hazardlib.geo.mesh` instance. The reference used for the calculation of azimuth is the middle point and the strike of the rupture. The value of azimuth computed corresponds to the angle measured in a clockwise direction from the strike of the rupture. :parameter mesh: An instance of :class:`openquake.hazardlib.geo.mesh` :return: An instance of `numpy.ndarray` """# Get info about the rupturestrike=self.get_strike()hypocenter=self.get_middle_point()# This is the azimuth from the north of each point Vs. the middle of# the ruptureazim=geodetic.azimuth(hypocenter.longitude,hypocenter.latitude,mesh.lons,mesh.lats)# Compute the azimuth from the fault strikerel_azi=(azim-strike)%360returnrel_azi
[docs]defget_azimuth_of_closest_point(self,mesh):""" Compute the azimuth between point in `mesh` and the corresponding closest point on the rupture surface. :param mesh: An instance of :class:`openquake.hazardlib.geo.mesh` :return: An :class:`numpy.ndarray` instance with the azimuth values. """mesh_closest=self.get_closest_points(mesh)returngeodetic.azimuth(mesh.lons,mesh.lats,mesh_closest.lons,mesh_closest.lats)