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
        )