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
"""
import numpy as np
from openquake.hazardlib.sourceconverter import SourceConverter
from openquake.hazardlib import nrml
from openquake.hazardlib.source.complex_fault import ComplexFaultSource
from openquake.hazardlib.source.characteristic import CharacteristicFaultSource
from openquake.hazardlib.source.simple_fault import SimpleFaultSource
from openquake.hazardlib.source.area import AreaSource
from openquake.hazardlib.source.point import PointSource
from openquake.hazardlib.geo.mesh import Mesh
from openquake.hazardlib.geo.polygon import Polygon


[docs]class RateGrid(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.): """ 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 = sources self.rates = np.zeros([self.nx - 1, self.ny - 1, self.nz - 1], dtype=float) self.area_discretisation = area_discretisation
[docs] @classmethod def from_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 = [] for grp in nrml.to_python(input_model, converter): sources.extend(grp.sources) return cls(limits, sources, area_discretisation)
[docs] def number_sources(self): """ Returns the number of sources """ return len(self.source_model)
[docs] def get_rates(self, mmin, mmax=np.inf): """ Returns the cumulative rates greater than Mmin :param float mmin: Minimum magnitude """ nsrcs = self.number_sources() for iloc, source in enumerate(self.source_model): print("Source Number %s of %s, Name = %s, Typology = %s" % ( iloc + 1, nsrcs, source.name, source.__class__.__name__)) if isinstance(source, CharacteristicFaultSource): self._get_fault_rates(source, mmin, mmax) elif isinstance(source, ComplexFaultSource): self._get_fault_rates(source, mmin, mmax) elif isinstance(source, SimpleFaultSource): self._get_fault_rates(source, mmin, mmax) elif isinstance(source, AreaSource): self._get_area_rates(source, mmin, mmax) elif isinstance(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]): return None, None xloc = int(((location.longitude - self.xlim[0]) / self.xspc) + 1E-7) if (location.latitude < self.ylim[0]) or\ (location.latitude > self.ylim[-1]): return None, None yloc = int(((location.latitude - self.ylim[0]) / self.yspc) + 1E-7) return xloc, yloc def _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 (xloc is None) or (yloc is None): return # Get annual rates annual_rate = source.get_annual_occurrence_rates() mags = np.array([val[0] for val in annual_rate]) annual_rate = np.array([val[1] for val in annual_rate]) idx = np.logical_and(mags >= mmin, mags < mmax) annual_rate = np.sum(annual_rate[idx]) for hypo_depth in source.hypocenter_distribution.data: zloc = int((hypo_depth[1] - self.zlim[0]) / self.zspc) if (zloc < 0) or (zloc >= (self.nz - 1)): continue else: 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) for point in points: 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 """ for rupt in list(source.iter_ruptures()): valid_rupt = (rupt.mag >= mmin) and (rupt.mag < mmax) if not valid_rupt: continue grd = 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]class RatePolygon(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.): """ 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) """ assert isinstance(limits["polygon"], Polygon) self.limits = limits["polygon"] self.upper_depth = limits["upper_depth"] self.lower_depth = limits["lower_depth"] self.source_model = sources self.rates = 0.0 self.area_discretisation = area_discretisation def _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] if not in_poly: return else: for (mag, rate) in source.get_annual_occurrence_rates(): if (mag < mmin) or (mag > mmax): return else: for (prob, depth) in source.hypocenter_distribution.data: if (depth < self.upper_depth) or\ (depth > self.lower_depth): continue else: self.rates += (prob * rate) 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 """ for rup in list(source.iter_ruptures()): if (rup.mag < mmin) or (rup.mag > mmax): # Magnitude outside search range continue depths = rup.surface.mesh.depths.flatten() # Generate simple mesh from surface rupt_mesh = Mesh(rup.surface.mesh.lons.flatten(), rup.surface.mesh.lats.flatten(), depths) # Mesh points in polygon in_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) if np.any(idx): node_rate = rup.occurrence_rate / float(len(depths)) self.rates += (node_rate * np.sum(idx))