Source code for openquake.hmtk.seismicity.smoothing.smoothed_seismicity

# -*- coding: utf-8 -*-
# vim: tabstop=4 shiftwidth=4 softtabstop=4

#
# LICENSE
#
# Copyright (C) 2010-2023 GEM Foundation, G. Weatherill, M. Pagani,
# D. Monelli.
#
# 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.

"""
Module :mod: openquake.hmtk.seismicity.smoothing.smoothed_seismicity implements the
:class: openquake.hmtk.seismicity.smoothing.smoothed_seismicity.SmoothedSeismicity,
a general class for implementing seismicity smoothing algorithms
"""
import csv

from math import fabs, log
import numpy as np
from openquake.hazardlib.geo.point import Point
from openquake.hazardlib.geo.polygon import Polygon
from openquake.hmtk.seismicity.smoothing import utils
from openquake.hmtk.seismicity.smoothing.kernels.isotropic_gaussian import (
    IsotropicGaussian,
)
from openquake.hmtk.registry import CatalogueFunctionRegistry


[docs]class Grid(dict):
[docs] @classmethod def make_from_list(cls, grid_limits): new = cls() new.update( { "xmin": grid_limits[0], "xmax": grid_limits[1], "xspc": grid_limits[2], "ymin": grid_limits[3], "ymax": grid_limits[4], "yspc": grid_limits[5], "zmin": grid_limits[6], "zmax": grid_limits[7], "zspc": grid_limits[8], } ) return new
[docs] @classmethod def make_from_catalogue(cls, catalogue, spacing, dilate): """ Defines the grid on the basis of the catalogue """ new = cls() cat_bbox = get_catalogue_bounding_polygon(catalogue) if dilate > 0: cat_bbox = cat_bbox.dilate(dilate) # Define Grid spacing new.update( { "xmin": np.min(cat_bbox.lons), "xmax": np.max(cat_bbox.lons), "xspc": spacing, "ymin": np.min(cat_bbox.lats), "ymax": np.max(cat_bbox.lats), "yspc": spacing, "zmin": 0.0, "zmax": np.max(catalogue.data["depth"]), "zspc": np.max(catalogue.data["depth"]), } ) if new["zmin"] == new["zmax"] == new["zspc"] == 0: new["zmax"] = new["zspc"] = 1 return new
[docs] def as_list(self): return [ self["xmin"], self["xmax"], self["xspc"], self["ymin"], self["ymax"], self["yspc"], self["zmin"], self["zmax"], self["zspc"], ]
[docs] def as_polygon(self): return Polygon( [ Point(self["xmin"], self["ymax"]), Point(self["xmax"], self["ymax"]), Point(self["xmax"], self["ymin"]), Point(self["xmin"], self["ymin"]), ] )
[docs] def dilate(self, width): polygon = self.as_polygon().dilate(width) self.update( { "xmin": np.min(polygon.lons), "xmax": np.max(polygon.lons), "ymin": np.min(polygon.lats), "ymax": np.max(polygon.lats), } ) return self
def _get_adjustment(mag, year, mmin, completeness_year, t_f, mag_inc=0.1): """ If the magnitude is greater than the minimum in the completeness table and the year is greater than the corresponding completeness year then return the Weichert factor :param float mag: Magnitude of an earthquake :param float year: Year of earthquake :param np.ndarray completeness_table: Completeness table :param float mag_inc: Magnitude increment :param float t_f: Weichert adjustment factor :returns: Weichert adjustment factor is event is in complete part of catalogue (0.0 otherwise) """ if len(completeness_year) == 1: if (mag >= mmin) and (year >= completeness_year[0]): # No adjustment needed - event weight == 1 return 1.0 else: # Event should not be counted return False kval = int(((mag - mmin) / mag_inc)) + 1 if (kval >= 1) and (year >= completeness_year[kval - 1]): return t_f else: return False
[docs]def get_catalogue_bounding_polygon(catalogue): """ Returns a polygon containing the bounding box of the catalogue """ upper_lon = np.max(catalogue.data["longitude"]) upper_lat = np.max(catalogue.data["latitude"]) lower_lon = np.min(catalogue.data["longitude"]) lower_lat = np.min(catalogue.data["latitude"]) return Polygon( [ Point(lower_lon, upper_lat), Point(upper_lon, upper_lat), Point(upper_lon, lower_lat), Point(lower_lon, lower_lat), ] )
[docs]class SmoothedSeismicity(object): """ Class to implement an analysis of Smoothed Seismicity, including the grid counting of data and the smoothing. :param np.ndarray grid: Observed count in each cell [Long., Lat., Depth., Count] :param catalogue: Valid instance of the :class: openquake.hmtk.seismicity.catalogue.Catalogue :param bool use_3d: Decide if analysis is 2-D (False) or 3-D (True). If 3-D then distances will use hypocentral distance, otherwise epicentral distance :param float bval: b-value :param float beta: Beta value for exponential form (beta = bval * log(10.)) :param np.ndarray data: Smoothed seismicity output :param dict grid_limits: Limits ot the grid used for defining the cells """ def __init__(self, grid_limits, use_3d=False, bvalue=None): """ Instatiate class with a set of grid limits :param grid_limits: It could be a float (in that case the grid is computed from the catalogue with the given spacing). Or an array of the form: [xmin, xmax, spcx, ymin, ymax, spcy, zmin, spcz] :param bool use_3d: Choose whether to use hypocentral distances for smoothing or only epicentral :param float bval: b-value for analysis """ self.grid = None self.catalogue = None self.use_3d = use_3d self.bval = bvalue if self.bval: self.beta = self.bval * log(10.0) else: self.beta = None self.data = None self.grid_limits = grid_limits self.kernel = None
[docs] def run_analysis( self, catalogue, config, completeness_table=None, smoothing_kernel=None ): """ Runs an analysis of smoothed seismicity in the manner originally implemented by Frankel (1995) :param catalogue: Instance of the openquake.hmtk.seismicity.catalogue.Catalogue class catalogue.data dictionary containing the following - 'year' - numpy.ndarray vector of years 'longitude' - numpy.ndarray vector of longitudes 'latitude' - numpy.ndarray vector of latitudes 'depth' - numpy.ndarray vector of depths :param dict config: Configuration settings of the algorithm: * 'Length_Limit' - Maximum number of bandwidths for use in smoothing (Float) * 'BandWidth' - Bandwidth (km) of the Smoothing Kernel (Float) * 'increment' - Output incremental (True) or cumulative a-value (False) :param np.ndarray completeness_table: Completeness of the catalogue assuming evenly spaced magnitudes from most recent bin to oldest bin [year, magnitude] :param smoothing_kernel: Smoothing kernel as instance of :class: `openquake.hmtk.seismicity.smoothing.kernels.base.BaseSmoothingKernel` :returns: Full smoothed seismicity data as np.ndarray, of the form [Longitude, Latitude, Depth, Observed, Smoothed] """ self.catalogue = catalogue if smoothing_kernel: self.kernel = smoothing_kernel else: self.kernel = IsotropicGaussian() # If no grid limits are specified then take from catalogue if isinstance(self.grid_limits, list): self.grid_limits = Grid.make_from_list(self.grid_limits) assert self.grid_limits["xmax"] >= self.grid_limits["xmin"] assert self.grid_limits["xspc"] > 0.0 assert self.grid_limits["ymax"] >= self.grid_limits["ymin"] assert self.grid_limits["yspc"] > 0.0 elif isinstance(self.grid_limits, float): self.grid_limits = Grid.make_from_catalogue( self.catalogue, self.grid_limits, config["Length_Limit"] * config["BandWidth"], ) completeness_table, mag_inc = utils.get_even_magnitude_completeness( completeness_table, self.catalogue ) end_year = self.catalogue.end_year # Get Weichert factor t_f, _ = utils.get_weichert_factor( self.beta, completeness_table[:, 1], completeness_table[:, 0], end_year, ) # Get the grid self.create_3D_grid(self.catalogue, completeness_table, t_f, mag_inc) if config["increment"]: # Get Hermann adjustment factors fval, fival = utils.hermann_adjustment_factors( self.bval, completeness_table[0, 1], config["increment"] ) self.data[:, -1] = fval * fival * self.data[:, -1] # Apply smoothing smoothed_data, sum_data, sum_smooth = self.kernel.smooth_data( self.data, config, self.use_3d ) print( "Smoothing Total Rate Comparison - " "Observed: %.6g, Smoothed: %.6g" % (sum_data, sum_smooth) ) self.data = np.column_stack([self.data, smoothed_data]) return self.data
[docs] def create_2D_grid_simple( self, longitude, latitude, year, magnitude, completeness_table, t_f=1.0, mag_inc=0.1, ): """ Generates the grid from the limits using an approach closer to that of Frankel (1995) :param numpy.ndarray longitude: Vector of earthquake longitudes :param numpy.ndarray latitude: Vector of earthquake latitudes :param numpy.ndarray year: Vector of earthquake years :param numpy.ndarray magnitude: Vector of earthquake magnitudes :param numpy.ndarray completeness_table: Completeness table :param float t_f: Weichert adjustment factor :returns: Two-dimensional spatial grid of observed rates """ assert mag_inc > 0.0 xlim = np.ceil( (self.grid_limits["xmax"] - self.grid_limits["xmin"]) / self.grid_limits["xspc"] ) ylim = np.ceil( (self.grid_limits["ymax"] - self.grid_limits["ymin"]) / self.grid_limits["yspc"] ) ncolx = int(xlim) ncoly = int(ylim) grid_count = np.zeros(ncolx * ncoly, dtype=float) for iloc in range(0, len(longitude)): dlon = ( longitude[iloc] - self.grid_limits["xmin"] ) / self.grid_limits["xspc"] if (dlon < 0.0) or (dlon > xlim): # Earthquake outside longitude limits continue xcol = int(dlon) if xcol == ncolx: # If longitude is directly on upper grid line then retain xcol = ncolx - 1 dlat = ( fabs(self.grid_limits["ymax"] - latitude[iloc]) / self.grid_limits["yspc"] ) if (dlat < 0.0) or (dlat > ylim): # Earthquake outside latitude limits continue ycol = int(dlat) # Correct for floating precision if ycol == ncoly: # If latitude is directly on upper grid line then retain ycol = ncoly - 1 kmarker = (ycol * int(xlim)) + xcol adjust = _get_adjustment( magnitude[iloc], year[iloc], completeness_table[0, 1], completeness_table[:, 0], t_f, mag_inc, ) if adjust: grid_count[kmarker] = grid_count[kmarker] + adjust return grid_count
def _xyz(self): x_bins = np.arange( self.grid_limits["xmin"], self.grid_limits["xmax"], self.grid_limits["xspc"], ) if x_bins[-1] < self.grid_limits["xmax"]: x_bins = np.hstack([x_bins, x_bins[-1] + self.grid_limits["xspc"]]) y_bins = np.arange( self.grid_limits["ymin"], self.grid_limits["ymax"], self.grid_limits["yspc"], ) if y_bins[-1] < self.grid_limits["ymax"]: y_bins = np.hstack([y_bins, y_bins[-1] + self.grid_limits["yspc"]]) z_bins = np.arange( self.grid_limits["zmin"], self.grid_limits["zmax"] + self.grid_limits["zspc"], self.grid_limits["zspc"], ) if z_bins[-1] < self.grid_limits["zmax"]: z_bins = np.hstack([z_bins, z_bins[-1] + self.grid_limits["zspc"]]) # Define centre points of grid cells grid_x, grid_y = np.meshgrid( (x_bins[1:] + x_bins[:-1]) / 2.0, (y_bins[1:] + y_bins[:-1]) / 2.0 ) return grid_x, grid_y, z_bins
[docs] def create_3D_grid(self, catalogue, completeness_table, t_f=1.0, mag_inc=0.1): """ Counts the earthquakes observed in a three dimensional grid :param catalogue: Instance of the openquake.hmtk.seismicity.catalogue.Catalogue class catalogue.data dictionary containing the following - 'year' - numpy.ndarray vector of years 'longitude' - numpy.ndarray vector of longitudes 'latitude' - numpy.ndarray vector of latitudes 'depth' - numpy.ndarray vector of depths :param np.ndarray completeness_table: Completeness of the catalogue assuming evenly spaced magnitudes from most recent bin to oldest bin [year, magnitude] :param float t_f: Weichert adjustment factor :param float mag_inc: Increment of the completeness magnitude (rendered 0.1) :returns: Three-dimensional spatial grid of observed rates (or two dimensional if only one depth layer is considered) """ gridx, gridy, z_bins = self._xyz() n_x, n_y = np.shape(gridx) gridx = np.reshape(gridx, [n_x * n_y, 1]) gridy = np.reshape(np.flipud(gridy), [n_x * n_y, 1]) # Only one depth range idx = np.logical_and( catalogue.data["depth"] >= z_bins[0], catalogue.data["depth"] < z_bins[1], ) mid_depth = (z_bins[0] + z_bins[1]) / 2.0 data_grid = np.column_stack( [ gridx, gridy, mid_depth * np.ones(n_x * n_y, dtype=float), self.create_2D_grid_simple( catalogue.data["longitude"][idx], catalogue.data["latitude"][idx], catalogue.data["year"][idx], catalogue.data["magnitude"][idx], completeness_table, t_f, mag_inc, ), ] ) if len(z_bins) < 3: # Only one depth range self.data = data_grid return # Multiple depth layers - append to grid for iloc in range(1, len(z_bins) - 1): idx = np.logical_and( catalogue.data["depth"] >= z_bins[iloc], catalogue.data["depth"] < z_bins[iloc + 1], ) mid_depth = (z_bins[iloc] + z_bins[iloc + 1]) / 2.0 temp_grid = np.column_stack( [ gridx, gridy, mid_depth * np.ones(n_x * n_y, dtype=float), self.create_2D_grid_simple( catalogue.data["longitude"][idx], catalogue.data["latitude"][idx], catalogue.data["year"][idx], catalogue.data["magnitude"][idx], completeness_table, t_f, mag_inc, ), ] ) data_grid = np.vstack([data_grid, temp_grid]) self.data = data_grid
[docs] def write_to_csv(self, filename): """ Exports to simple csv :param str filename: Path to file for export """ fid = open(filename, "wt") # Create header list header_info = [ "Longitude", "Latitude", "Depth", "Observed Count", "Smoothed Rate", "b-value", ] writer = csv.DictWriter(fid, fieldnames=header_info) headers = dict((name0, name0) for name0 in header_info) # Write to file writer.writerow(headers) for row in self.data: # institute crude compression by omitting points with no seismicity # and taking advantage of the %g format if row[4] == 0: continue row_dict = { "Longitude": "%g" % row[0], "Latitude": "%g" % row[1], "Depth": "%g" % row[2], "Observed Count": "%d" % row[3], "Smoothed Rate": "%.6g" % row[4], "b-value": "%g" % self.bval, } writer.writerow(row_dict) fid.close()
SMOOTHED_SEISMICITY_METHODS = CatalogueFunctionRegistry()
[docs]@SMOOTHED_SEISMICITY_METHODS.add( "run", completeness=True, b_value=float, use_3d=bool, grid_limits=Grid, Length_Limit=float, BandWidth=float, increment=bool, ) class IsotropicGaussianMethod(object):
[docs] def run(self, catalogue, config, completeness=None): ss = SmoothedSeismicity( config["grid_limits"], config["use_3d"], config["b_value"] ) return ss.run_analysis( catalogue, config, completeness_table=completeness )