Source code for openquake.hmtk.comparison.rate_grids
# -*- coding: utf-8 -*-# vim: tabstop=4 shiftwidth=4 softtabstop=4## LICENSE## Copyright (C) 2015-2023 GEM Foundation, G. Weatherill, M. Pagani## The Hazard Modeller's Toolkit 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.## You should have received a copy of the GNU Affero General Public License# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>## DISCLAIMER## The software Hazard Modeller's Toolkit (openquake.hmtk) provided herein# is released as a prototype implementation on behalf of# scientists and engineers working within the GEM Foundation (Global# Earthquake Model).## It is distributed for the purpose of open collaboration and in the# hope that it will be useful to the scientific, engineering, disaster# risk and software design communities.## The software is NOT distributed as part of GEM's OpenQuake suite# (https://www.globalquakemodel.org/tools-products) and must be considered as a# separate entity. The software provided herein is designed and implemented# by scientific staff. It is not developed to the design standards, nor# subject to same level of critical review by professional software# developers, as GEM's OpenQuake software suite.## Feedback and contribution to the software is welcome, and can be# directed to the hazard scientific staff of the GEM Model Facility# (hazard@globalquakemodel.org).## The Hazard Modeller's Toolkit (openquake.hmtk) is therefore distributed WITHOUT# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or# FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License# for more details.## The GEM Foundation, and the authors of the software, assume no# liability for use of the software."""Python tools for calculating activity rates on a grid from a source model"""importnumpyasnpfromopenquake.hazardlib.sourceconverterimportSourceConverterfromopenquake.hazardlibimportnrmlfromopenquake.hazardlib.source.complex_faultimportComplexFaultSourcefromopenquake.hazardlib.source.characteristicimportCharacteristicFaultSourcefromopenquake.hazardlib.source.simple_faultimportSimpleFaultSourcefromopenquake.hazardlib.source.areaimportAreaSourcefromopenquake.hazardlib.source.pointimportPointSourcefromopenquake.hazardlib.geo.meshimportMeshfromopenquake.hazardlib.geo.polygonimportPolygon
[docs]classRateGrid(object):""" Class for calculation of activity rate grids :param float xspc: Longitude spacing of grid :param float yspc: Latitude spacing of grid :param float zspc: Depth spacing (km) of grid :param np.ndarray xlim: Longitude cell bounds :param np.ndarray ylim: Latitude cell bounds :param np.ndarray zlim: Depth cell bounds :param int nx: Number of longitude cells :param int ny: Number of latitude cells :param int nz: Number of depth cells :param list source_model: Seismic source mode :param np.ndarray rates: Activity rates :param float area_discretisation: Discretisation step (km) of area sources """def__init__(self,limits,sources,area_discretisation=10.0):""" Instantiate class with grid configurations :param list limits: Grid configuration [west, east, xspc, south, north, yspc, upper, lower, zspc] """self.xspc=limits[2]self.yspc=limits[5]self.zspc=limits[8]self.xlim=np.arange(limits[0],limits[1]+self.xspc,self.xspc)self.ylim=np.arange(limits[3],limits[4]+self.yspc,self.yspc)self.zlim=np.arange(limits[6],limits[7]+self.zspc,self.zspc)self.nx=len(self.xlim)self.ny=len(self.ylim)self.nz=len(self.zlim)self.source_model=sourcesself.rates=np.zeros([self.nx-1,self.ny-1,self.nz-1],dtype=float)self.area_discretisation=area_discretisation
[docs]@classmethoddeffrom_model_files(cls,limits,input_model,investigation_time=1.0,simple_mesh_spacing=1.0,complex_mesh_spacing=5.0,mfd_width=0.1,area_discretisation=10.0,):""" Reads the hazard model from a file :param list limits: Grid configuration [west, east, xspc, south, north, yspc, upper, lower, zspc] :param str input_model: Path to input source model :param float investigation_time: Investigation time of Poisson model :param float simple_mesh_spacing: Rupture mesh spacing of simple fault (km) :param float complex_mesh_spacing: Rupture mesh spacing of complex fault (km) :param float mfd_width: Spacing (in magnitude units) of MFD :param float area_discretisation: Spacing of discretisation of area source (km) """converter=SourceConverter(investigation_time,simple_mesh_spacing,complex_mesh_spacing,mfd_width,area_discretisation,)sources=[]forgrpinnrml.to_python(input_model,converter):sources.extend(grp.sources)returncls(limits,sources,area_discretisation)
[docs]defnumber_sources(self):""" Returns the number of sources """returnlen(self.source_model)
[docs]defget_rates(self,mmin,mmax=np.inf):""" Returns the cumulative rates greater than Mmin :param float mmin: Minimum magnitude """nsrcs=self.number_sources()foriloc,sourceinenumerate(self.source_model):print("Source Number %s of %s, Name = %s, Typology = %s"%(iloc+1,nsrcs,source.name,source.__class__.__name__))ifisinstance(source,CharacteristicFaultSource):self._get_fault_rates(source,mmin,mmax)elifisinstance(source,ComplexFaultSource):self._get_fault_rates(source,mmin,mmax)elifisinstance(source,SimpleFaultSource):self._get_fault_rates(source,mmin,mmax)elifisinstance(source,AreaSource):self._get_area_rates(source,mmin,mmax)elifisinstance(source,PointSource):self._get_point_rates(source,mmin,mmax)else:print("Source type %s not recognised - skipping!"%source)continue
def_get_point_location(self,location):""" Returns the location in the output grid corresponding to the cell in which the epicentre lays :param location: Source hypocentre as instance of :class: openquake.hazardlib.geo.point.Point :returns: xloc - Location of longitude cell yloc - Location of latitude cell """if(location.longitude<self.xlim[0])or(location.longitude>self.xlim[-1]):returnNone,Nonexloc=int(((location.longitude-self.xlim[0])/self.xspc)+1e-7)if(location.latitude<self.ylim[0])or(location.latitude>self.ylim[-1]):returnNone,Noneyloc=int(((location.latitude-self.ylim[0])/self.yspc)+1e-7)returnxloc,ylocdef_get_point_rates(self,source,mmin,mmax=np.inf):""" Adds the rates for a point source :param source: Point source as instance of :class: openquake.hazardlib.source.point.PointSource :param float mmin: Minimum Magnitude :param float mmax: Maximum Magnitude """xloc,yloc=self._get_point_location(source.location)if(xlocisNone)or(ylocisNone):return# Get annual ratesannual_rate=source.get_annual_occurrence_rates()mags=np.array([val[0]forvalinannual_rate])annual_rate=np.array([val[1]forvalinannual_rate])idx=np.logical_and(mags>=mmin,mags<mmax)annual_rate=np.sum(annual_rate[idx])forhypo_depthinsource.hypocenter_distribution.data:zloc=int((hypo_depth[1]-self.zlim[0])/self.zspc)if(zloc<0)or(zloc>=(self.nz-1)):continueelse:self.rates[xloc,yloc,zloc]+=(float(hypo_depth[0])*annual_rate)def_get_area_rates(self,source,mmin,mmax=np.inf):""" Adds the rates from the area source by discretising the source to a set of point sources :param source: Area source as instance of :class: openquake.hazardlib.source.area.AreaSource """points=list(source)forpointinpoints:self._get_point_rates(point,mmin,mmax)def_get_fault_rates(self,source,mmin,mmax=np.inf):""" Adds the rates for a simple or complex fault source :param source: Fault source as instance of :class: openquake.hazardlib.source.simple_fault.SimpleFaultSource or openquake.hazardlib.source.complex_fault.ComplexFaultSource """forruptinlist(source.iter_ruptures()):valid_rupt=(rupt.mag>=mmin)and(rupt.mag<mmax)ifnotvalid_rupt:continuegrd=np.column_stack([rupt.surface.mesh.lons.flatten(),rupt.surface.mesh.lats.flatten(),rupt.surface.mesh.depths.flatten(),])npts=np.shape(grd)[0]counter=np.histogramdd(grd,bins=[self.xlim,self.ylim,self.zlim])[0]point_rate=rupt.occurrence_rate/float(npts)self.rates+=point_rate*counter
[docs]classRatePolygon(RateGrid):""" Calculates the rate of events within a polygon :param limits: Polygon as instance of :class: openquake.hazardlib.geo.polygon.Polygon :param float upper_depth: Upper seismic depth of the polygon (km) :param float lower_depth: Lower seismic depth of the polygon (km) :param list source_model: List of seismic sources :param float rates: Activity rate of polygon :param float area_discretisation: Discretisation spacing (km) of the area source """def__init__(self,limits,sources,area_discretisation=10.0):""" Instantiate class with grid configurations :param dict limits: Configuration as dictionary containing: * polygon - OpenQuake Polygon * uppper_depth - upper seismogenic depth (km) * lower_depth - lower seismogenic depth (km) """assertisinstance(limits["polygon"],Polygon)self.limits=limits["polygon"]self.upper_depth=limits["upper_depth"]self.lower_depth=limits["lower_depth"]self.source_model=sourcesself.rates=0.0self.area_discretisation=area_discretisationdef_get_point_rates(self,source,mmin,mmax=np.inf):""" Adds the rates for a point source :param source: Point source as instance of :class: openquake.hazardlib.source.point.PointSource :param float mmin: Minimum Magnitude :param float mmax: Maximum Magnitude """src_mesh=Mesh.from_points_list([source.location])in_poly=self.limits.intersects(src_mesh)[0]ifnotin_poly:returnelse:formag,rateinsource.get_annual_occurrence_rates():if(mag<mmin)or(mag>mmax):returnelse:forprob,depthinsource.hypocenter_distribution.data:if(depth<self.upper_depth)or(depth>self.lower_depth):continueelse:self.rates+=prob*ratedef_get_fault_rates(self,source,mmin,mmax=np.inf):""" Adds the rates for a simple or complex fault source :param source: Fault source as instance of :class: openquake.hazardlib.source.simple_fault.SimpleFaultSource or openquake.hazardlib.source.complex_fault.ComplexFaultSource """forrupinlist(source.iter_ruptures()):if(rup.mag<mmin)or(rup.mag>mmax):# Magnitude outside search rangecontinuedepths=rup.surface.mesh.depths.flatten()# Generate simple mesh from surfacerupt_mesh=Mesh(rup.surface.mesh.lons.flatten(),rup.surface.mesh.lats.flatten(),depths,)# Mesh points in polygonin_poly=self.limits.intersects(rupt_mesh)in_depth=np.logical_and(depths>=self.upper_depth,depths<=self.lower_depth)idx=np.logical_and(in_poly,in_depth)ifnp.any(idx):node_rate=rup.occurrence_rate/float(len(depths))self.rates+=node_rate*np.sum(idx)