From e61ff9d3a39703714090661b21b3e4b82ef5e467 Mon Sep 17 00:00:00 2001
From: Graeme Weatherill <gweather@gfz-potsdam.de>
Date: Thu, 11 Mar 2021 15:42:56 +0000
Subject: [PATCH] Implements initial complete shakemap feature
---
shakyground2/earthquake.py | 96 +++--
shakyground2/rupture_mechanism.py | 50 ++-
shakyground2/shakemap.py | 286 +++++++++++++++
shakyground2/site_model.py | 14 +-
shakyground2/synthetic_rupture_generator.py | 388 ++++++++++++++++++++
shakyground2/valid.py | 19 +-
tests/data/__init__.py | 0
tests/earthquake_test.py | 39 +-
tests/rupture_mechanism_test.py | 59 ++-
tests/shakemap_test.py | 112 ++++++
tests/synthetic_rupture_generator_test.py | 335 +++++++++++++++++
11 files changed, 1327 insertions(+), 71 deletions(-)
create mode 100644 shakyground2/shakemap.py
create mode 100644 shakyground2/synthetic_rupture_generator.py
create mode 100644 tests/data/__init__.py
create mode 100644 tests/shakemap_test.py
create mode 100644 tests/synthetic_rupture_generator_test.py
diff --git a/shakyground2/earthquake.py b/shakyground2/earthquake.py
index 7dc90b8..6aff1c1 100644
--- a/shakyground2/earthquake.py
+++ b/shakyground2/earthquake.py
@@ -2,12 +2,13 @@
Contains the core classes to represent an earthquake
"""
import datetime
-from typing import Tuple
+from typing import Tuple, Optional, List
from math import radians, sin, tan
from openquake.hazardlib.geo import Point, PlanarSurface
from openquake.hazardlib.geo.surface.base import BaseSurface
from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
from shakyground2 import valid
+from shakyground2.rupture_mechanism import RuptureMechanism
class Earthquake(object):
@@ -28,6 +29,9 @@ class Earthquake(object):
Earthquake latitude in decimal degrees (-90, to 90)
depth:
Hypocentral depth (in km)
+ hypocenter:
+ Representation of the hypocenter as an instance of
+ :class:`openquake.hazardlib.geo.Point`
mag:
Earthquake magnitude
strike:
@@ -63,14 +67,15 @@ class Earthquake(object):
lat,
hypo_depth,
magnitude,
- strike=0.0,
- dip=90.0,
- rake=0.0,
+ strike=None,
+ dip=None,
+ rake=None,
aspect=1.0,
mag_scaling_relation=None,
upper_seismogenic_depth=0.0,
lower_seismogenic_depth=1000.0,
surface=None,
+ focal_mechanism=None,
hypocenter_position=(0.5, 0.5),
date=None,
time=None,
@@ -121,6 +126,7 @@ class Earthquake(object):
self.lon = valid.longitude(lon)
self.lat = valid.latitude(lat)
self.depth = valid.positive_float(hypo_depth, "hypocentre depth")
+ self.hypocenter = Point(self.lon, self.lat, self.depth)
self.mag = magnitude
self.strike = valid.strike(strike)
self.dip = valid.dip(dip)
@@ -138,8 +144,41 @@ class Earthquake(object):
self.time = time
assert isinstance(surface, BaseSurface) or surface is None
self.surface = surface
+ if self.surface:
+ # Can calculate rupture dimensions from the surface
+ self.area = self.surface.get_area()
+ self.width = self.surface.get_width()
+ self.length = self.area / self.width
+ else:
+ # Default rupture dimensions to none to none
+ self.area = self.width = self.length = None
self.hypo_pos = valid.hypocenter_position(hypocenter_position)
self._rupture = None
+ # Get a valid focal mechanism with two nodal planes
+ self.focal_mechanism = valid.focal_mechanism(focal_mechanism)
+ self.mechanism = self._get_mechanism()
+ self.tectonic_region = None
+
+ def _get_mechanism(self):
+ """
+ Defines the focal mechanism according to three different cases:
+ 1. A unique mechanism is defined explicitly from the strike, dip and rake
+ 2. A pair of equiprobable mechanisms is defined from the focal mechanism
+ 3. No mechanism is defined, in which case a global distribution is assumed
+
+ Returns:
+ Mechanism distribution as an instance of the
+ :class:`shakyground2.rupture_mechanism.RuptureMechanism`
+ """
+ if (self.strike is not None) and (self.dip is not None) and (self.rake is not None):
+ # Fixed and fully defined rupture mechanism
+ return RuptureMechanism.from_strike_dip_rake(self.strike, self.dip, self.rake)
+ elif self.focal_mechanism:
+ # Rupture mechanism defines from nodal planes
+ return RuptureMechanism.from_focal_mechanism(self.focal_mechanism)
+ else:
+ # Global distribution
+ return RuptureMechanism()
def __repr__(self):
# Returns a summary string of the event, with datetime if specified
@@ -156,6 +195,34 @@ class Earthquake(object):
self.id, self.lon, self.lat, self.depth, self.mag
)
+ def get_maximum_distance_bbox(self, max_distance: Optional[float] = None) -> List:
+ """
+ Defines a bounding box around the event up to a maximum distance.
+
+ Args:
+ max_distance: Maximum horizontal and vertical distance from the epicentre for the
+ bounding box to be defined. If None then a default bounding box
+ size is determined that is based on magnitude between Mw 3.0 (100 km)
+ and Mw 8.0 (1000 km)
+ Returns:
+ Bounding box as a list of [llon, llat, ulon, ulat]
+ """
+ if not max_distance:
+ # Scale from 100 km (for Mw 3 or less) to 1000 km for Mw >= 8.0
+ if self.mag <= 3.0:
+ max_distance = 100.0
+ elif self.mag >= 8.0:
+ max_distance = 1000.0
+ else:
+ # Interpolate
+ max_distance = 100.0 + (self.mag - 3.0) * (900.0 / 5.0)
+ # Define the bounding box from the specified maximum distance
+ north = self.hypocenter.point_at(max_distance, 0.0, 0.0)
+ south = self.hypocenter.point_at(max_distance, 0.0, 180.0)
+ east = self.hypocenter.point_at(max_distance, 0.0, 90.0)
+ west = self.hypocenter.point_at(max_distance, 0.0, 270.0)
+ return [west.longitude, south.latitude, east.longitude, north.latitude]
+
@property
def rupture(self):
"""
@@ -172,27 +239,6 @@ class Earthquake(object):
self.mag, self.rake, None, centroid, self.surface, 1.0, None
)
return self._rupture
- # Rupture is not defined, so needs to be constructed. Use simple case for now
- # TODO To be replaced by the synthetic rupture generator operation
- thickness = self.lsd - self.usd
- # Get rupture dimensions
- area, length, width = self.mag_scale_rel.get_rupture_dimensions(
- self.mag, aspect=self.aspect, dip=self.dip, thickness=thickness
- )
- # Build the rupture surface
- surface = self.build_planar_surface(
- centroid,
- self.strike,
- self.dip,
- length,
- width,
- self.lsd,
- self.usd,
- self.hypo_pos,
- )
- self._rupture = ParametricProbabilisticRupture(
- self.mag, self.rake, None, centroid, surface, 1.0, None
- )
return self._rupture
@staticmethod
diff --git a/shakyground2/rupture_mechanism.py b/shakyground2/rupture_mechanism.py
index 03c43e7..62fb339 100644
--- a/shakyground2/rupture_mechanism.py
+++ b/shakyground2/rupture_mechanism.py
@@ -4,7 +4,7 @@ Class to manage the distribution of rupture mechanisms
from __future__ import annotations
import numpy as np
from enum import Enum
-from typing import Optional, Dict, Tuple
+from typing import Optional, Dict, Tuple, List
from openquake.hazardlib.geo import NodalPlane
from openquake.hazardlib.pmf import PMF
from shakyground2 import valid
@@ -56,6 +56,10 @@ class RuptureMechanism(object):
for prob, nodal_plane in self.mechanism.data:
yield prob, nodal_plane
+ def __len__(self):
+ # Lenth corresponds to the number of mechanisms in the distribution
+ return len(self.mechanism.data)
+
def sample(self, nsamples: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""
Samples the focal mechanism distribution to return "nsamples" strike,
@@ -100,30 +104,48 @@ class RuptureMechanism(object):
)
@classmethod
- def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict) -> RuptureMechanism:
+ def from_focal_mechanism(cls, focal_mechanism: Dict) -> RuptureMechanism:
"""
- Constructs the rupture mechanism distribution from either one single
- valid rupture plane or two valud nodal planes
+ Constructs the focal mechanism from an evenly pair of nodal planes such as that
+ of a focal mechanism
Args:
- nodal_plane_1: Dictionary containing strike, dip and rake of the first nodel plane
- nodal_plane_1: Dictionary containing strike, dip and rake of the second nodel plane
+ focal_mechanism: Dictionary containing two :class:
+ `openquake.hazardlib.geo.NodalPlane` objects, each labeled as "nodal_plane_1"
+ and "nodal_plane_2"
"""
- strike_1 = valid.strike(nodal_plane_1["strike"])
- dip_1 = valid.dip(nodal_plane_1["dip"])
- rake_1 = valid.rake(nodal_plane_1["rake"])
- strike_2 = valid.strike(nodal_plane_2["strike"])
- dip_2 = valid.dip(nodal_plane_2["dip"])
- rake_2 = valid.rake(nodal_plane_2["rake"])
return cls(
PMF(
[
- (0.5, NodalPlane(strike_1, dip_1, rake_1)),
- (0.5, NodalPlane(strike_2, dip_2, rake_2)),
+ (0.5, focal_mechanism["nodal_plane_1"]),
+ (0.5, focal_mechanism["nodal_plane_2"]),
]
)
)
+ @classmethod
+ def from_nodal_planes(cls, nodal_planes: List, probabilities: List) -> RuptureMechanism:
+ """
+ Constructs the rupture mechanism distribution from a list of nodal planes and their
+ associated probabilities
+
+ Args:
+ nodal_planes: Set of nodal planes as a list of dictionaries, eac containing strike,
+ dip and rake
+ probabilities: List of probabilities of the nodal planes (must sum to 1)
+ """
+ assert len(nodal_planes) == len(
+ probabilities
+ ), "Number of nodal planes not equal to number of probabilities"
+ assert np.isclose(
+ sum(probabilities), 1.0
+ ), "Probabilities do not sum to 1.0 (sum = %.6f)" % sum(probabilities)
+ mechanism_distribution = []
+ for prob, npl in zip(probabilities, nodal_planes):
+ mechanism = valid.mechanism(npl["strike"], npl["dip"], npl["rake"])
+ mechanism_distribution.append((prob, mechanism))
+ return cls(PMF(mechanism_distribution))
+
@staticmethod
def build_mechanism_distribution(
strike: Optional[float] = None,
diff --git a/shakyground2/shakemap.py b/shakyground2/shakemap.py
new file mode 100644
index 0000000..9a460d2
--- /dev/null
+++ b/shakyground2/shakemap.py
@@ -0,0 +1,286 @@
+"""
+Implements the core shakemap class
+"""
+import h5py
+import numpy as np
+import warnings
+from typing import Dict, Optional, Tuple, List
+from openquake.hazardlib import const, imt
+from openquake.hazardlib.contexts import ContextMaker
+from shakyground2 import valid
+from shakyground2.earthquake import Earthquake
+from shakyground2.site_model import SiteModel
+from shakyground2.synthetic_rupture_generator import FiniteRuptureSampler
+
+
+class Shakemap(object):
+ """
+ The core class for calculating shakemaps. The general workflow contained within
+ takes the input earthquake and site model in order to build the source, path and site
+ context objects (i.e. containers that hold the data needed for the ground motion models).
+
+ The contexts objects are built upon instantiation of the shakemap class to avoid
+ recalculation if further shakemaps are needed.
+
+ The Shakemap class also caches the information into an hdf5 binary object, if desired by
+ the user. This means that the shakemap information for an earthquake ccould be retrieved at
+ a later time without the need for recalculation.
+ """
+
+ def __init__(
+ self,
+ earthquake: Earthquake,
+ site_model: SiteModel,
+ ground_motion_model: Dict,
+ tectonic_region: str,
+ cache_file: Optional[str] = None,
+ num_rupture_samples: int = 100,
+ rdim: float = 0.0,
+ synth_dist_weights: List = [0.25, 0.25, 0.25, 0.25],
+ synthetic_rupture_max_site_distance: float = 200.0,
+ synthetic_rupture_site_spacing: float = 0.05,
+ ):
+ """
+ Args:
+ earthquake:
+ Earthquake as an instance of the :class:shakyground2.earthquake.Earthquake
+ site_model:
+ Target sites used for calculation of the ground motion values, as instances of
+ :class:shakyground2.site_model.SiteModel
+ ground_motion_model:
+ Set of ground motion models and their respective weights as a dictionary
+ {"GMMs": [...], "weights": [...]}
+ tectonic_region:
+ Tectonic region to which the earthquake is assigned
+ cache_file:
+ Path to the hdf5 file for caching the results. It does not need to exist
+ but if the hdf5 file exists and contains an instance of the same event
+ (according to the earthquake.id) then the user will be warned here and the
+ existing shakemaps for that earthquake ID subsequently overwitten if the
+ user persists to implement the shakemap.
+ num_rupture_samples:
+ In the likely case that no rupture surface is available for the earthquake,
+ one will be generated by the :class:`shakyground2.synthetic_rupture_generator.
+ FiniteRuptureSampler`. This controls the number of samples to be used
+ rdim:
+ Exponent of the distance-dependent site weighing adopted by the :class:
+ `shakyground2.synthetic_rupture_generator.FiniteRuptureSampler`
+
+ synth_dist_weights: Can adjust the weight assigned to the different distance
+ metrics used by the FiniteRuptureSampler
+ synthetic_rupture max_site_distance: In the case that the target sites for the
+ FiniteRuptureSampler need to be built from the
+ bounding box, this defines the maximum
+ distances of the target sites to define the
+ bounding box
+ synthetic_rupture_site_spacing: In the case that the target sites for the
+ FiniteRuptureSampler need to be built from the
+ bounding box, this defines the site spacing of the
+ sites
+ """
+ self.earthquake = earthquake
+ self.site_model = site_model
+ self.tectonic_region = tectonic_region
+ self.ground_motion_model = ground_motion_model
+ self.num_rupture_samples = num_rupture_samples
+ self.rdim = rdim
+ self.synthetic_rupture_distance_weights = synth_dist_weights
+ if cache_file:
+ self.caching = True
+ self.cache_file = cache_file
+ # Check if the earthquake is cleady in the case and warn the user if so
+ fle = h5py.File(self.cache_file, "a")
+ if earthquake.id in list(fle):
+ warnings.warn(
+ "Earthquake %s already in cache file %s - "
+ "Running the shakemaps will overwrite this"
+ % (self.earthquake.id, self.cache_file)
+ )
+ fle.close()
+ else:
+ self.caching = False
+ self.cache_file = None
+ self.rctx = None
+ self.dctx = None
+ self.sctx = None
+ self.synth_rupture_max_site_dist = valid.positive_float(
+ synthetic_rupture_max_site_distance, "Max. Synthetic Rupture Site Distance"
+ )
+ self.synth_rupture_site_spacing = valid.positive_float(
+ synthetic_rupture_site_spacing, "Synthetic Rupture Site Spacing"
+ )
+ self._build_contexts()
+
+ def _build_contexts(self):
+ """
+ Construct the rupture, site and distances contexts from the earthquake, site and
+ ground motion models
+ """
+ cmaker = ContextMaker(self.tectonic_region, self.ground_motion_model["GMMs"])
+
+ if not self.earthquake.rupture:
+ # Use the finite rupture sampler to generate the rupture and corresponding
+ # distances from the available information
+ self.earthquake._rupture = FiniteRuptureSampler().get_finite_rupture(
+ self.num_rupture_samples,
+ self.earthquake,
+ rdim=self.rdim,
+ weights=self.synthetic_rupture_distance_weights,
+ maximum_site_distance=self.synth_rupture_max_site_dist,
+ site_spacing=self.synth_rupture_site_spacing,
+ )[0]
+ # For the sites and rupture context we can use the context maker to get all of the
+ # source, site and rupture distances
+ self.rctx, self.sctx, self.dctx = cmaker.make_contexts(
+ self.site_model.get_site_collection(), self.earthquake.rupture
+ )
+
+ def _cache_contexts(self, grp: h5py.Group):
+ """
+ If caching the shakemaps, then this stores the context information to the file
+ """
+ # Add the contexts to the group object
+ ctxt = grp.create_group("contexts")
+ rup_ctxt = ctxt.create_group("rupture")
+ dist_ctxt = ctxt.create_group("distances")
+ for gmm in self.ground_motion_model["GMMs"]:
+ for rup_attr in gmm.REQUIRES_RUPTURE_PARAMETERS:
+ if rup_attr not in rup_ctxt.attrs:
+ rup_ctxt.attrs[rup_attr] = getattr(self.rctx, rup_attr)
+ for attr in gmm.REQUIRES_DISTANCES:
+ if attr not in list(dist_ctxt):
+ distance = getattr(self.dctx, attr)
+ dist_dset = dist_ctxt.create_dataset(attr, distance.shape, dtype="f")
+ dist_dset[:] = distance
+ site_ctxt = ctxt.create_dataset(
+ "sites", self.site_model.site_array.shape, dtype=self.site_model.site_array.dtype
+ )
+ site_ctxt[:] = self.site_model.site_array
+ if self.site_model.bbox_properties:
+ # If the site model has bounding box properties then store these
+ # as attributes
+ site_ctxt.attrs["has_bbox"] = True
+ site_ctxt.attrs["llon"] = self.site_model.bbox_properties["bbox"][0]
+ site_ctxt.attrs["llat"] = self.site_model.bbox_properties["bbox"][1]
+ site_ctxt.attrs["ulon"] = self.site_model.bbox_properties["bbox"][2]
+ site_ctxt.attrs["ulat"] = self.site_model.bbox_properties["bbox"][3]
+ site_ctxt.attrs["spcx"] = self.site_model.bbox_properties["spcx"]
+ site_ctxt.attrs["spcy"] = self.site_model.bbox_properties["spcy"]
+ site_ctxt.attrs["ncol"] = self.site_model.bbox_properties["ncol"]
+ site_ctxt.attrs["nrow"] = self.site_model.bbox_properties["nrow"]
+ else:
+ site_ctxt.attrs["has_bbox"] = False
+ return
+
+ def get_shakemap(
+ self, intensity_measure_types: List
+ ) -> Tuple[np.ndarray, np.ndarray, Dict]:
+ """
+ Main function to constructs the shakemaps for the specified intensity measures,
+ caching the information to hdf5 if requested
+
+ Args:
+ intensity_measure_types: List of intensity measures for which the shakemaps will
+ be calculated
+
+ Returns:
+ aggregated_means: The mean of the ground motions from the different ground motion
+ models weighted by the assigned input weights
+ aggregated_stddevs: The total standard deviation of the ground motions from the
+ different ground motion models weighted by the assigned
+ input weights
+ shakemaps: Dictionary of individual shakemaps for each ground motion model
+ """
+ shakemaps = {}
+ shakemap_dtypes = np.dtype(
+ [
+ (intensity_measure_type, np.float64)
+ for intensity_measure_type in intensity_measure_types
+ ]
+ )
+ if self.caching:
+ # If caching, open (or create) the file and cache the contexts
+ fle = h5py.File(self.cache_file, "r+")
+ if self.earthquake.id in list(fle):
+ del fle[self.earthquake.id]
+ fle_eq = fle.create_group(self.earthquake.id)
+ self._cache_contexts(fle_eq)
+
+ # Pre-allocate the aggregated shakemaps with zeros
+ aggregated_means = np.zeros(self.site_model.shape, dtype=shakemap_dtypes)
+ aggregated_stddevs = np.zeros(self.site_model.shape, dtype=shakemap_dtypes)
+ for weight, gmm in zip(
+ self.ground_motion_model["weights"], self.ground_motion_model["GMMs"]
+ ):
+ gmm_str = gmm.__class__.__name__
+ shakemaps[gmm_str] = {
+ "weight": weight,
+ "mean": np.zeros(self.site_model.shape, dtype=shakemap_dtypes),
+ "stddev": np.zeros(self.site_model.shape, dtype=shakemap_dtypes),
+ }
+ for intensity_measure_type in intensity_measure_types:
+ input_imt = imt.from_string(intensity_measure_type)
+ try:
+ mean, [stddev] = gmm.get_mean_and_stddevs(
+ self.sctx, self.rctx, self.dctx, input_imt, [const.StdDev.TOTAL]
+ )
+ except KeyError:
+ warnings.warn(
+ "Ground motion model %s not defined for intensity "
+ "measure type %s" % (str(gmm), intensity_measure_type)
+ )
+ continue
+ aggregated_means[intensity_measure_type] += weight * mean
+ aggregated_stddevs[intensity_measure_type] += weight * stddev
+ shakemaps[gmm_str]["mean"][intensity_measure_type] = mean
+ shakemaps[gmm_str]["stddev"][intensity_measure_type] = stddev
+ if self.caching:
+ self._cache_shakemap(
+ fle_eq, shakemaps, aggregated_means, aggregated_stddevs, shakemap_dtypes
+ )
+ fle.close()
+ return aggregated_means, aggregated_stddevs, shakemaps
+
+ def _cache_shakemap(
+ self,
+ fle: h5py.Group,
+ shakemaps: Dict,
+ aggregated_means: np.ndarray,
+ aggregated_stddevs: np.ndarray,
+ shakemap_dtypes: np.dtype,
+ ):
+ """
+ If caching is required then the shakemaps are written to the hdf5 file
+
+ Args:
+ fle: HDF5 group object to store the shakemaps
+ shakemaps: Dictionary of individual shakemaps for each ground motion model
+ aggregated_means: The mean of the ground motions from the different ground motion
+ models weighted by the assigned input weights
+ aggregated_stddevs: The total standard deviation of the ground motions from the
+ different ground motion models weighted by the assigned
+ input weights
+ shakemap_dtypes: IMT-dependent data type of the shakemaps
+ """
+ shakemap_grp = fle.create_group("shakemaps")
+ for gmm_string in shakemaps:
+ gmm_grp = shakemap_grp.create_group(gmm_string)
+ gmm_grp.attrs["weight"] = shakemaps[gmm_string]["weight"]
+ mean_dset = gmm_grp.create_dataset(
+ "mean", shakemaps[gmm_string]["mean"].shape, dtype=shakemap_dtypes
+ )
+ mean_dset[:] = shakemaps[gmm_string]["mean"]
+ stddev_dset = gmm_grp.create_dataset(
+ "stddev", shakemaps[gmm_string]["stddev"].shape, dtype=shakemap_dtypes
+ )
+ stddev_dset[:] = shakemaps[gmm_string]["stddev"]
+ # Store the agregated results
+ aggregated_grp = fle.create_group("aggregated")
+ aggregated_mean_dset = aggregated_grp.create_dataset(
+ "mean", aggregated_means.shape, shakemap_dtypes
+ )
+ aggregated_mean_dset[:] = aggregated_means
+ aggregated_stddev_dset = aggregated_grp.create_dataset(
+ "stddev", aggregated_stddevs.shape, shakemap_dtypes
+ )
+ aggregated_stddev_dset[:] = aggregated_stddevs
diff --git a/shakyground2/site_model.py b/shakyground2/site_model.py
index e835e0a..7ce9daa 100644
--- a/shakyground2/site_model.py
+++ b/shakyground2/site_model.py
@@ -5,7 +5,7 @@ from __future__ import annotations
import warnings
import numpy as np
import pandas as pd
-from typing import Tuple, Optional, Dict
+from typing import Tuple, Optional, Dict, Union
from openquake.hazardlib.site import SiteCollection
@@ -56,6 +56,14 @@ class SiteModel(object):
# Returns the number of sites in the site array
return self.site_array.shape[0]
+ def __getitem__(self, key: Union[str, int]):
+ # Returns the named column from the site array (if provided with a string),
+ # or the corresponding row if provided with an integer
+ if isinstance(key, str) or isinstance(key, int):
+ return self.site_array[key]
+ else:
+ raise KeyError("Site model has no attribute %s" % key)
+
def __repr__(self):
# Returns an informative string regarding the composition of the sites,
# including the bounding box information when relevant
@@ -72,6 +80,10 @@ class SiteModel(object):
info_string += bbox_string
return info_string
+ @property
+ def shape(self):
+ return self.site_array.shape[0]
+
@property
def dataframe(self):
# Return the site array as a dataframe
diff --git a/shakyground2/synthetic_rupture_generator.py b/shakyground2/synthetic_rupture_generator.py
new file mode 100644
index 0000000..e628cb8
--- /dev/null
+++ b/shakyground2/synthetic_rupture_generator.py
@@ -0,0 +1,388 @@
+"""
+Generates a set of synthetic ruptures based on a point source configuration
+"""
+import numpy as np
+from typing import List, Tuple, Optional
+from scipy.stats import truncnorm
+from openquake.hazardlib.geo import Point, Mesh
+from openquake.hazardlib.geo.surface.base import BaseSurface
+from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
+from shakyground2.earthquake import Earthquake
+from shakyground2.site_model import SiteModel
+
+
+# The hypocentre distributions are defined from the Next Generation Attenuation
+# (NGA) West 2 rupture database
+NGAWEST_HYPO_DISTRIBUTIONS = {
+ "All": {
+ "Mean": np.array([0.46851357, 0.67460261]),
+ "Sigma": np.array([0.20691535, 0.22040932]),
+ "COV": np.array([[0.04293594, 0.00103987], [0.00103987, 0.04871868]]),
+ "DipRange": [30.0, 90.0],
+ },
+ "SS": {
+ "Mean": np.array([0.48545206, 0.64942746]),
+ "Sigma": np.array([0.22415657, 0.21677068]),
+ "COV": np.array([[0.05064814, -0.00603003], [-0.00603003, 0.04736544]]),
+ "DipRange": [80.0, 90.0],
+ },
+ "R": {
+ "Mean": np.array([0.4674859, 0.58483914]),
+ "Sigma": np.array([0.16275562, 0.22017015]),
+ "COV": np.array([[0.02673021, 0.0113362], [0.0113362, 0.04891558]]),
+ "DipRange": [25.0, 50.0],
+ },
+ "N": {
+ "Mean": np.array([0.50887254, 0.82404]),
+ "Sigma": np.array([0.22416128, 0.13647917]),
+ "COV": np.array([[0.05085368, -0.00332741], [-0.00332741, 0.01885098]]),
+ "DipRange": [50.0, 80.0],
+ },
+}
+
+
+# In case other definitions (thrust fault, undefined etc.) are used
+NGAWEST_HYPO_DISTRIBUTIONS["TF"] = NGAWEST_HYPO_DISTRIBUTIONS["R"]
+NGAWEST_HYPO_DISTRIBUTIONS["U"] = NGAWEST_HYPO_DISTRIBUTIONS["All"]
+
+
+class FiniteRuptureSampler(object):
+ """
+ Module to produce a set of finite rupture samples for an event
+ Process:
+ 1. From scaling relation and seismogenic thicknesses sample rupture
+ parameters (area, length, width)
+ 2. Sample hypocentre location from multivariate Gaussian distribution
+ 3. a. If mechanism is known, use strike and dip
+ b. If mechanism is not known sample strike and dip from distribution
+ 4. For each sample calculate Rjb, Rrup, Rx and Ry0
+ 5. For each distance calculate the median value
+ 6. For each rupture define a penalty function
+
+ P = sum_{i=rjb, rrup, ry0, rx} (Ri - median(Ri)) ** 2.
+ Take the rupture with the lowest penalty function!
+ """
+
+ def get_finite_rupture(
+ self,
+ nsamples: int,
+ earthquake: Earthquake,
+ site_model: Optional[SiteModel] = None,
+ rdim: float = 0.0,
+ weights: list = [0.25, 0.25, 0.25, 0.25],
+ maximum_site_distance: float = 200.0,
+ site_spacing: float = 0.05,
+ ):
+ """
+ Retrieves the finite rupture that produces the distances closest to the median
+ distance
+
+ Args:
+ nsamples: Number of rupture samples to use
+ earthquake: Earthquake object as instance of
+ :class:`shakyground2.earthquake.Earthquake`
+ site_model: If desired, can define the target sites to be used for the distance
+ calculation, as instance of :class:`shakyground2.site_model.SiteModel`
+ rdim: The penalty functions can be weighted according to the source-to-site
+ distance with the weighting determined via w = 1 / (distance ^ R). rdim
+ sets the value of R
+ weights: Can adjust the weight assigned to the different distance metrics. Defines
+ the list of weights for the four distance metrics: Rjb, Rrup, Rx and Ry0
+ maximum_site_distance: In the case that the target sites need to be built from the
+ bounding box, this defines the maximum distances of the
+ target sites to define the bounding box
+ site_spacing: In the case that the target sites need to be built from the
+ bounding box, this defines the site spacing
+
+ Returns:
+ central_rupture: The rupture corresponding to the `most central` as an instance of
+ :class:`openquake.hazardlib.source.rupture.ParametricProbabilisticRupture`
+ distances: Dictionary containing all of the possible distances calculated from
+ the rupture surface
+ """
+ # Sample the rupture surfaces
+ rupture_surfaces, strikes, dips, rakes, hypo_locs = self.sample_rupture_surfaces(
+ nsamples, earthquake
+ )
+
+ if site_model:
+ target_lons = site_model["lon"]
+ target_lats = site_model["lat"]
+ else:
+ # If no specific site model is defined then use a bounding box within 200 km of
+ # the hypocentre spaced every 0.05 degrees - longer distances are less likely
+ # to be significantly influenced by uncertainties in the rupture plane
+ llon, llat, ulon, ulat = earthquake.get_maximum_distance_bbox(maximum_site_distance)
+ target_lons, target_lats = np.meshgrid(
+ np.arange(llon, ulon + site_spacing, site_spacing),
+ np.arange(llat, ulat + site_spacing, site_spacing),
+ )
+ target_lons = target_lons.flatten()
+ target_lats = target_lats.flatten()
+
+ # Calculate the distances
+ distances, rhypo, repi = self.get_distances(
+ earthquake, rupture_surfaces, target_lons, target_lats
+ )
+ # Get the most central rupture closest to the median finite rupture distance
+ # from all finite ruptures
+ central_surface, central_distances, min_loc = self.get_central_rupture(
+ rupture_surfaces, distances, rhypo, rdim, weights
+ )
+
+ # Return the parametric probabilitic rupture and the corrsponding distance
+ central_rupture = ParametricProbabilisticRupture(
+ earthquake.mag,
+ rakes[min_loc],
+ None, # No tectonic region type
+ earthquake.hypocenter,
+ central_surface,
+ 1.0,
+ None,
+ )
+
+ return central_rupture, {
+ "lon": target_lons,
+ "lat": target_lats,
+ "rjb": central_distances[:, 0],
+ "rrup": central_distances[:, 1],
+ "rx": central_distances[:, 2],
+ "ry0": central_distances[:, 3],
+ "rhypo": rhypo,
+ "repi": repi,
+ }
+
+ def sample_rupture_surfaces(
+ self, nsamples: int, earthquake: Earthquake
+ ) -> Tuple[List, np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+ """
+ Mechanisms should be a list of Nodal plane objects
+
+ Args:
+ nsamples: Number of samples to take
+ earthquake: Earthquake object as instance of
+ :class:`shakyground2.earthquake.Earthquake`
+ Returns:
+ surfaces: List of rupture surfaces as instances of the :class:
+ `openquake.hazardlib.geo.PlanarSurface
+ strikes: Vector of strike values (in decimal degrees in the range 0 to 360)
+ dips: Vector of dip values (in decimal degrees in the range 0 to 90)
+ rakes: Vector of rake values (in decimal degrees in the range -180 to 180)
+ hypo_positions: 2D array of hypocentre positions (as along-strike and down-dip)
+ positions
+ """
+ strikes, dips, rakes = earthquake.mechanism.sample(nsamples)
+ dimensions = (
+ (earthquake.area is not None)
+ and (earthquake.length is not None)
+ and (earthquake.width is not None)
+ )
+ if dimensions:
+ areas = earthquake.area * np.ones(nsamples)
+ lengths = earthquake.length * np.ones(nsamples)
+ widths = earthquake.width * np.ones(nsamples)
+ else:
+ # Get the physical rupture dimensions by sampling from the magnitude scaling
+ # relation
+ areas, lengths, widths = self.sample_rupture_dimensions(
+ nsamples, earthquake, dips, rakes
+ )
+
+ # Get hypocentres
+ hypo_positions = self.sample_hypocenter_position(nsamples, self.get_sofs(rakes))
+
+ # Build the ruptures
+ surfaces = []
+ for strike, dip, length, width, hypo_pos in zip(
+ strikes, dips, lengths, widths, hypo_positions
+ ):
+ surfaces.append(
+ earthquake.build_planar_surface(
+ earthquake.hypocenter,
+ strike,
+ dip,
+ length,
+ width,
+ earthquake.lsd,
+ earthquake.usd,
+ hypo_pos,
+ )
+ )
+ return surfaces, strikes, dips, rakes, hypo_positions
+
+ @staticmethod
+ def get_sofs(rakes: np.ndarray) -> List:
+ """
+ Returns the style of faulting as a string given an input list of rakes
+ """
+ sofs = []
+ for rake in rakes:
+ if (rake > 45.0) and (rake <= 135.0):
+ sofs.append("R")
+ elif (rake > -135.0) and (rake <= -45.0):
+ sofs.append("N")
+ else:
+ sofs.append("SS")
+ return sofs
+
+ @staticmethod
+ def sample_rupture_dimensions(
+ nsamples: int, earthquake: Earthquake, dips: np.ndarray, rakes: np.ndarray
+ ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+ """
+ Samples the rupture dimension from the magnitude scaling relation accounting for
+ uncertainty where it is defined
+
+ Args:
+ nsamples: Number of samples to take
+ earthquake: Earthquake object as instance of
+ :class:`shakyground2.earthquake.Earthquake`
+ dips: Vector of dip values
+ rakes: Vector of rake values
+ Returns:
+ areas: Vector of sampled areas (in km ^ 2)
+ lengths: Vector of sampled lengths (in km)
+ widths: Vector of sampled widths (in km)
+ """
+ # Sample epsilon values from a truncated normal distribution between
+ # -3 and 3 standard deviations
+ msr_epsilons = truncnorm.rvs(-3.0, 3.0, loc=0.0, scale=1.0, size=nsamples)
+ # Generate rupture dimensions
+ lengths = np.empty(nsamples)
+ widths = np.empty(nsamples)
+ areas = np.empty(nsamples)
+ thickness = earthquake.lsd - earthquake.usd
+ for i, (dip, rake, epsilon) in enumerate(zip(dips, rakes, msr_epsilons)):
+ areas[i], lengths[i], widths[i] = earthquake.mag_scale_rel.get_rupture_dimensions(
+ earthquake.mag,
+ rake,
+ dip,
+ aspect=earthquake.aspect,
+ thickness=thickness,
+ epsilon=epsilon,
+ )
+ return areas, lengths, widths
+
+ @staticmethod
+ def sample_hypocenter_position(nsamples: int, sofs: List) -> np.ndarray:
+ """
+ Samples the hypocentral position from the pre-defined NGA West 2 hypocentre
+ distributions
+
+ Note here that the hypocentre position distributions are multivariate Gaussian, yet
+ they need to be truncated within the range 0 to 1. Sampling from a truncated
+ multivariate Gaussian distribution is not a trivial problem and requires a Markov Chain
+ Monte Carlo sampling approach that is too complex for application here. Instead,
+ invalid hypocentre positions are re-sampled repeatedly until all positions are valid
+
+ Args:
+ nsamples: Number of samples to take
+ sofs: List of styles-of-faulting
+ Returns:
+ 2D array of hypocentre position samples (along strike and down dip)
+ """
+ # Count the unique styles of faulting present in the sampling set
+ sof_indices = {}
+ sof_counter = {}
+ sof_array = np.array(sofs)
+ for key in NGAWEST_HYPO_DISTRIBUTIONS:
+ sof_indices[key] = np.where(sof_array == key)[0]
+ sof_counter[key] = len(sof_indices[key])
+
+ hypo_pos = np.zeros([nsamples, 2])
+ for key in sof_counter:
+ if not sof_counter[key]:
+ # Style-of-faulting not present in sample set
+ continue
+ mean = NGAWEST_HYPO_DISTRIBUTIONS[key]["Mean"]
+ covariance = NGAWEST_HYPO_DISTRIBUTIONS[key]["COV"]
+ # Draw N samples from a multivariate distribution, where N
+ # is the corresponding number of samples with the required style of
+ # faulting
+ pos = np.random.multivariate_normal(mean, covariance, sof_counter[key])
+ # Are all positions valid?
+ all_valid = np.logical_and(
+ np.logical_and(pos[:, 0] >= 0.05, pos[:, 0] <= 0.95),
+ np.logical_and(pos[:, 1] >= 0.05, pos[:, 1] <= 0.95),
+ )
+
+ while np.sum(all_valid) < sof_counter[key]:
+ # Some positions are invalid - re-sample these
+ idx = np.logical_not(all_valid)
+ pos[idx] = np.random.multivariate_normal(mean, covariance, sum(idx))
+ all_valid = np.logical_and(
+ np.logical_and(pos[:, 0] >= 0.05, pos[:, 0] <= 0.95),
+ np.logical_and(pos[:, 1] >= 0.05, pos[:, 1] <= 0.95),
+ )
+
+ hypo_pos[sof_indices[key], :] = pos
+ return hypo_pos
+
+ @staticmethod
+ def get_distances(
+ earthquake: Earthquake, surfaces: List, lons: np.ndarray, lats: np.ndarray
+ ) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
+ """
+ Calculates the source to site distances for each rupture surface in
+ the set of rupture surfaces
+
+ Args:
+ earthquake: Earthquake object as instance of
+ :class:`shakyground2.earthquake.Earthquake`
+ lons: Vector of target longitudes for distance calculation
+ lats: Vector of target latitudes for distance calculation
+ Returns:
+ distances: 2D array of finite fault rupture distances in km (Rjb, Rrup, Rx, Ry0)
+ rhypo: hypocentral distance in km
+ repi: epicentral distance in km
+ """
+ nsites = len(lons)
+ mesh = Mesh(lons, lats)
+ distances = np.zeros([len(surfaces), nsites, 4])
+ for i, surf in enumerate(surfaces):
+ distances[i, :, 0] = surf.get_joyner_boore_distance(mesh)
+ distances[i, :, 1] = surf.get_min_distance(mesh)
+ distances[i, :, 2] = surf.get_rx_distance(mesh)
+ distances[i, :, 3] = surf.get_ry0_distance(mesh)
+ # Include hypocentral and epicentral distance
+ rhypo = Point(earthquake.lon, earthquake.lat, earthquake.depth).distance_to_mesh(mesh)
+ repi = np.sqrt(rhypo ** 2.0 - earthquake.depth ** 2.0)
+ return distances, rhypo, repi
+
+ @staticmethod
+ def get_central_rupture(
+ rupture_surfaces: List,
+ distances: np.ndarray,
+ hypocentral_distances: np.ndarray,
+ rdim: float = 0.0,
+ weights: List = [0.25, 0.25, 0.25, 0.25],
+ ) -> Tuple[BaseSurface, np.ndarray, int]:
+ """
+ Returns the rupture closest to the median distances
+ """
+ nsites = distances.shape[1]
+ # If desired, it is possible to weight the distribution of distances
+ # such that the weight decays with distance from the centroid, with the
+ # exponent of the decay factor described by rdim (if 0 then no weighting is applied).
+ site_weights = 1.0 / (hypocentral_distances ** rdim)
+ site_weights = site_weights / np.sum(site_weights)
+
+ # For each of the four different distance types retrieve the median source-to-site
+ # distance for each site
+ median_distances = np.zeros([4, nsites])
+ for i in range(4):
+ median_distances[i, :] = np.percentile(distances[:, :, i], 50, axis=0)
+ # For each of the rupture surfaces determine the degree of divergence from the
+ # median distance and assign a penalty that is proportional to the divergence
+ penalty_function = np.zeros(distances.shape[0])
+ for i in range(distances.shape[2]):
+ site_penalty = np.zeros(distances.shape[0])
+ for k in range(distances.shape[1]):
+ site_penalty += site_weights[k] * (
+ np.sqrt((distances[:, k, i] - median_distances[i, k]) ** 2.0)
+ )
+ penalty_function += weights[i] * site_penalty
+ # Determine which rupture gives the minimal divergence from the median and return
+ # that rupture and corresponding distances
+ min_loc = np.argmin(penalty_function)
+ return rupture_surfaces[min_loc], distances[min_loc, :, :], min_loc
diff --git a/shakyground2/valid.py b/shakyground2/valid.py
index 07c769f..99424cd 100644
--- a/shakyground2/valid.py
+++ b/shakyground2/valid.py
@@ -2,7 +2,7 @@
Defines a set of input validation methods to check physically correct or
consistent quantities
"""
-from typing import Tuple, Optional, Union
+from typing import Tuple, Optional, Union, Dict
from openquake.hazardlib.geo.nodalplane import NodalPlane
from shakyground2.magnitude_scaling_relations import BaseScalingRelation, PEERScalingRelation
@@ -83,6 +83,23 @@ def mechanism(
return NodalPlane(*mechanism)
+def focal_mechanism(focal_mech: Optional[Dict]) -> Dict:
+ """
+ A focal mechanism is represented by two orthogonal planes, each plane described by a
+ strike, dip and rake
+ """
+ if not focal_mech:
+ return focal_mech
+ assert "nodal_plane_1" in focal_mech, "Focal mechanism missing nodal plane 1"
+ assert "nodal_plane_2" in focal_mech, "Focal mechanism missing nodal plane 2"
+ focal_mechanism = {}
+ for plane in ["nodal_plane_1", "nodal_plane_2"]:
+ focal_mechanism[plane] = mechanism(
+ focal_mech[plane]["strike"], focal_mech[plane]["dip"], focal_mech[plane]["rake"]
+ )
+ return focal_mechanism
+
+
def seismogenic_thickness(
upper_seismo_depth: float, lower_seismo_depth: float
) -> Tuple[float, float]:
diff --git a/tests/data/__init__.py b/tests/data/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/tests/earthquake_test.py b/tests/earthquake_test.py
index 038303d..108fd63 100644
--- a/tests/earthquake_test.py
+++ b/tests/earthquake_test.py
@@ -191,27 +191,20 @@ class EarthquakeTestCase(unittest.TestCase):
np.testing.assert_array_almost_equal(
plane0.corner_depths, eq0.rupture.surface.corner_depths
)
- # Instantiate earthquake object without surface
- eq1 = Earthquake(
- "XYZ",
- self.lon,
- self.lat,
- self.depth,
- self.mag,
- strike=strike,
- dip=dip,
- rake=0.0,
- upper_seismogenic_depth=usd,
- lower_seismogenic_depth=lsd,
- )
- self.assertAlmostEqual(eq1.rupture.mag, self.mag)
- self.assertAlmostEqual(eq1.rupture.rake, 0.0)
- np.testing.assert_array_almost_equal(
- plane0.corner_lons, eq1.rupture.surface.corner_lons
- )
- np.testing.assert_array_almost_equal(
- plane0.corner_lats, eq1.rupture.surface.corner_lats
- )
+
+ def test_maximum_distance_bbox(self):
+ # Tests the generation of a bounding box to a maximum distance
+ # Test case 1 - Defined maximum distance
+ eq0 = Earthquake("XYZ", 0.0, 0.0, 10.0, 6.0)
+ bbox = eq0.get_maximum_distance_bbox(111.195)
np.testing.assert_array_almost_equal(
- plane0.corner_depths, eq1.rupture.surface.corner_depths
- )
+ np.array(bbox), np.array([-1.0, -1.0, 1.0, 1.0]), 5
+ )
+ # Test case 2 - Default magnitude values
+ for mag, max_dist in zip([2.5, 6.0, 8.5], [200.0, 1280.0, 2000.0]):
+ eq0 = Earthquake("XYZ", 0.0, 0.0, 10.0, mag)
+ bbox = eq0.get_maximum_distance_bbox()
+ ew_distance = Point(bbox[0], 0.0, 0.0).distance(Point(bbox[2], 0.0, 0.0))
+ ns_distance = Point(0.0, bbox[1], 0.0).distance(Point(0.0, bbox[3], 0.0))
+ self.assertAlmostEqual(max_dist, ew_distance, 5)
+ self.assertAlmostEqual(max_dist, ns_distance, 5)
diff --git a/tests/rupture_mechanism_test.py b/tests/rupture_mechanism_test.py
index 79a7f3d..d3ad2b5 100644
--- a/tests/rupture_mechanism_test.py
+++ b/tests/rupture_mechanism_test.py
@@ -3,6 +3,7 @@ Unit tests for rupture mechanism
"""
import unittest
import numpy as np
+from openquake.hazardlib.geo import NodalPlane
from shakyground2.rupture_mechanism import RuptureMechanism
@@ -42,17 +43,21 @@ class RuptureMechanismTestCase(unittest.TestCase):
def test_rupture_mechanism_single_plane(self):
# Checks instantiation when a single combination of strike, dip and rake are input
mech = RuptureMechanism.from_strike_dip_rake(strike=20.0, dip=90.0, rake=10.0)
+ self.assertEqual(len(mech), 1)
probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
np.testing.assert_array_almost_equal(probs, np.array([1.0]))
np.testing.assert_array_almost_equal(strikes, np.array([20.0]))
np.testing.assert_array_almost_equal(dips, np.array([90.0]))
np.testing.assert_array_almost_equal(rakes, np.array([10.0]))
- def test_rupture_mechanism_from_nodal_planes(self):
+ def test_rupture_mechanism_from_focal_mechanism(self):
# Checks the creation of the class from a pair of nodal planes
- npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.0}
- npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
- mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
+ focal_mechanism = {
+ "nodal_plane_1": NodalPlane(strike=45.0, dip=80.0, rake=10.0),
+ "nodal_plane_2": NodalPlane(strike=135.0, dip=90.0, rake=-170.0),
+ }
+ mech = RuptureMechanism.from_focal_mechanism(focal_mechanism)
+ self.assertEqual(len(mech), 2)
probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
np.testing.assert_array_almost_equal(probs, np.array([0.5, 0.5]))
np.testing.assert_array_almost_equal(strikes, np.array([45.0, 135.0]))
@@ -61,9 +66,11 @@ class RuptureMechanismTestCase(unittest.TestCase):
def test_sample_mechanism(self):
# Test the random sampling
- npl1 = {"strike": 45.0, "dip": 80.0, "rake": 10.0}
- npl2 = {"strike": 135.0, "dip": 90.0, "rake": -170.0}
- mech = RuptureMechanism.from_nodal_planes(npl1, npl2)
+ focal_mechanism = {
+ "nodal_plane_1": NodalPlane(strike=45.0, dip=80.0, rake=10.0),
+ "nodal_plane_2": NodalPlane(strike=135.0, dip=90.0, rake=-170.0),
+ }
+ mech = RuptureMechanism.from_focal_mechanism(focal_mechanism)
np.random.seed(1000)
strikes, dips, rakes = mech.sample(5)
np.testing.assert_array_almost_equal(
@@ -73,3 +80,41 @@ class RuptureMechanismTestCase(unittest.TestCase):
np.testing.assert_array_almost_equal(
rakes, np.array([-170.0, 10.0, -170.0, 10.0, -170.0])
)
+
+ def test_rupture_mechanism_from_nodal_planes(self):
+ # Tests the creation of the class from a list of nodal planes
+ nodal_planes = [
+ {"strike": 0.0, "dip": 90.0, "rake": 0.0},
+ {"strike": 90.0, "dip": 90.0, "rake": 0.0},
+ {"strike": 0.0, "dip": 60.0, "rake": -90.0},
+ {"strike": 90.0, "dip": 60.0, "rake": -90.0},
+ ]
+ probabilities = [0.3, 0.3, 0.2, 0.2]
+ mech = RuptureMechanism.from_nodal_planes(nodal_planes, probabilities)
+ self.assertEqual(len(mech), 4)
+ probs, strikes, dips, rakes = self._mechanism_pmf_to_arrays(mech)
+ np.testing.assert_array_almost_equal(probs, np.array(probabilities))
+ np.testing.assert_array_almost_equal(strikes, np.array([0.0, 90.0, 0.0, 90.0]))
+ np.testing.assert_array_almost_equal(dips, np.array([90.0, 90.0, 60.0, 60.0]))
+ np.testing.assert_array_almost_equal(rakes, np.array([0.0, 0.0, -90.0, -90.0]))
+
+ def test_rupture_mechanism_from_nodal_planes_errors(self):
+ # Test that assertion errors are raised when incorrect inputs are used
+ # 1) Unequal number of planes and probabilities
+ nodal_planes = [
+ {"strike": 0.0, "dip": 90.0, "rake": 0.0},
+ {"strike": 90.0, "dip": 90.0, "rake": 0.0},
+ {"strike": 0.0, "dip": 60.0, "rake": -90.0},
+ {"strike": 90.0, "dip": 60.0, "rake": -90.0},
+ ]
+ probabilities = [0.3, 0.3, 4]
+ with self.assertRaises(AssertionError) as ae:
+ RuptureMechanism.from_nodal_planes(nodal_planes, probabilities)
+ self.assertEqual(
+ str(ae.exception), "Number of nodal planes not equal to number of probabilities"
+ )
+ # 2) Probabilities don't sum to 1
+ probabilities = [0.3, 0.3, 0.2, 0.1]
+ with self.assertRaises(AssertionError) as ae:
+ RuptureMechanism.from_nodal_planes(nodal_planes, probabilities)
+ self.assertEqual(str(ae.exception), "Probabilities do not sum to 1.0 (sum = 0.900000)")
diff --git a/tests/shakemap_test.py b/tests/shakemap_test.py
new file mode 100644
index 0000000..fdba4b7
--- /dev/null
+++ b/tests/shakemap_test.py
@@ -0,0 +1,112 @@
+"""
+Set of test cases for the shakemap class
+"""
+import os
+import unittest
+import h5py
+import numpy as np
+from shakyground2.earthquake import Earthquake
+from shakyground2.site_model import SiteModel
+from shakyground2.shakemap import Shakemap
+from openquake.hazardlib.gsim import get_available_gsims
+
+
+# List of available GMMs
+GSIM_LIST = get_available_gsims()
+
+
+# Path for storing the cache file that will be generated
+DATA_PATH = os.path.join(os.path.dirname(__file__), "data")
+
+
+class ShakemapTestCase(unittest.TestCase):
+ """
+ The Shakemap class effectively draws together different components of the shakyground
+ tools in order to build the inputs for the OpenQuake ground motion models (GMMs) and
+ execute them to retrieve the shakemap. Each component is already tested, including the
+ OpenQuake GMMs themselves, verifying the actual values in the shakemap is not worthwhile
+ here. Instead the tests verify that the shakemap workflow is executed without error and
+ that the outputs are of the correct form and shape.
+ """
+
+ def setUp(self):
+ self.cache_test_file = os.path.join(DATA_PATH, "tmp_cache_file1.hdf5")
+ self.focal_mechanism = {
+ "nodal_plane_1": {"strike": 0.0, "dip": 90.0, "rake": 0.0},
+ "nodal_plane_2": {"strike": 90.0, "dip": 90.0, "rake": 0.0},
+ }
+ self.site_model = SiteModel.from_bbox(
+ [28.0, 38.0, 32.0, 42.0], spcx=0.5, spcy=0.5, vs30=600.0, vs30measured=True
+ )
+ self.ground_motion_model = {
+ "GMMs": [
+ GSIM_LIST["BindiEtAl2014Rjb"](),
+ GSIM_LIST["CauzziEtAl2014"](),
+ GSIM_LIST["ChiouYoungs2014"](),
+ ],
+ "weights": [0.3, 0.4, 0.3],
+ }
+
+ def test_shakemap_creation_with_caching(self):
+ # Complete run through workflow with caching, verifying that the information stored in
+ # the cache file is the same (within the tolerated precision) as that output directly
+ eq_id = "XYZ001"
+ eq1 = Earthquake(eq_id, 30.0, 40.0, 10.0, 6.5, focal_mechanism=self.focal_mechanism)
+ sm1 = Shakemap(
+ eq1,
+ self.site_model,
+ self.ground_motion_model,
+ "Active Shallow Crust",
+ cache_file=self.cache_test_file,
+ )
+ intensity_measures = ["PGA", "SA(0.3)", "SA(1.0)"]
+ agg_means, agg_stddev, shakemaps = sm1.get_shakemap(intensity_measures)
+ # Verify that shakemaps have been stored in the cache correctly
+ fle = h5py.File(self.cache_test_file, "r")
+ for intensity_measure in intensity_measures:
+ np.testing.assert_array_almost_equal(
+ agg_means[intensity_measure],
+ fle["{:s}/aggregated/mean".format(eq_id)][:][intensity_measure],
+ 5,
+ )
+ np.testing.assert_array_almost_equal(
+ agg_stddev[intensity_measure],
+ fle["{:s}/aggregated/stddev".format(eq_id)][:][intensity_measure],
+ 5,
+ )
+ for gmm in shakemaps:
+ for intensity_measure in intensity_measures:
+ np.testing.assert_array_almost_equal(
+ shakemaps[gmm]["mean"][intensity_measure],
+ fle["{:s}/shakemaps/{:s}/mean".format(eq_id, gmm)][:][intensity_measure],
+ 5,
+ )
+ fle.close()
+
+ def test_shakemap_creation_no_caching(self):
+ # Complete run through workflow verifying that the outputs have the correct content
+ # and shape (actual values are not verified)
+ fle = h5py.File(self.cache_test_file, "a")
+ fle.close()
+ eq1 = Earthquake("XYZ001", 30.0, 40.0, 10.0, 6.5, focal_mechanism=self.focal_mechanism)
+ sm1 = Shakemap(eq1, self.site_model, self.ground_motion_model, "Active Shallow Crust")
+ intensity_measures = ["PGA", "SA(0.3)", "SA(1.0)"]
+ agg_means, agg_stddev, shakemaps = sm1.get_shakemap(intensity_measures)
+ # Assert that the aggregated means and stddevs have the correct shape and dtypes
+ self.assertTupleEqual(agg_means.shape, self.site_model.site_array.shape)
+ self.assertTupleEqual(agg_stddev.shape, self.site_model.site_array.shape)
+ self.assertTupleEqual(agg_means.dtype.names, ("PGA", "SA(0.3)", "SA(1.0)"))
+ self.assertListEqual(
+ list(shakemaps), ["BindiEtAl2014Rjb", "CauzziEtAl2014", "ChiouYoungs2014"]
+ )
+ # Verify correct aggregation
+ expected_aggregation_pga = (
+ shakemaps["BindiEtAl2014Rjb"]["mean"]["PGA"] * 0.3
+ + shakemaps["CauzziEtAl2014"]["mean"]["PGA"] * 0.4
+ + shakemaps["ChiouYoungs2014"]["mean"]["PGA"] * 0.3
+ )
+ np.testing.assert_array_almost_equal(agg_means["PGA"], expected_aggregation_pga)
+
+ def tearDown(self):
+ # Delete the cache file
+ os.remove(self.cache_test_file)
diff --git a/tests/synthetic_rupture_generator_test.py b/tests/synthetic_rupture_generator_test.py
new file mode 100644
index 0000000..cdaeb5e
--- /dev/null
+++ b/tests/synthetic_rupture_generator_test.py
@@ -0,0 +1,335 @@
+"""
+Test cases for the synthetic finite rupture generator. Note that the rupture generator
+uses random sampling so we choose to unittest using fixed seeding
+"""
+import unittest
+import numpy as np
+import pandas as pd
+from openquake.hazardlib.geo import Point, PlanarSurface
+from openquake.hazardlib.source.rupture import ParametricProbabilisticRupture
+from shakyground2.earthquake import Earthquake
+from shakyground2.site_model import SiteModel
+from shakyground2.magnitude_scaling_relations import Stafford2014
+from shakyground2.synthetic_rupture_generator import FiniteRuptureSampler
+
+
+class FiniteRuptureSamplerTestCase(unittest.TestCase):
+ """"""
+
+ def setUp(self):
+ self.sampler = FiniteRuptureSampler()
+ # Define a set of 5 parallel rupture surfaces
+ top_left = Point(29.5, 40.0, 1.0)
+ top_right = Point(30.5, 40.0, 1.0)
+ bottom_left = Point(29.5, 40.0, 15.0)
+ bottom_right = Point(30.5, 40.0, 15.0)
+ sfc1 = PlanarSurface.from_corner_points(
+ top_left.point_at(22.0, 0.0, 0.0),
+ top_right.point_at(22.0, 0.0, 0.0),
+ bottom_right.point_at(22.0, 0.0, 0.0),
+ bottom_left.point_at(22.0, 0.0, 0.0),
+ )
+ sfc2 = PlanarSurface.from_corner_points(
+ top_left.point_at(11.0, 0.0, 0.0),
+ top_right.point_at(11.0, 0.0, 0.0),
+ bottom_right.point_at(11.0, 0.0, 0.0),
+ bottom_left.point_at(11.0, 0.0, 0.0),
+ )
+ sfc3 = PlanarSurface.from_corner_points(top_left, top_right, bottom_right, bottom_left)
+ sfc4 = PlanarSurface.from_corner_points(
+ top_left.point_at(11.0, 0.0, 180.0),
+ top_right.point_at(11.0, 0.0, 180.0),
+ bottom_right.point_at(11.0, 0.0, 180.0),
+ bottom_left.point_at(11.0, 0.0, 180.0),
+ )
+ sfc5 = PlanarSurface.from_corner_points(
+ top_left.point_at(22.0, 0.0, 180.0),
+ top_right.point_at(22.0, 0.0, 180.0),
+ bottom_right.point_at(22.0, 0.0, 180.0),
+ bottom_left.point_at(22.0, 0.0, 180.0),
+ )
+ self.rupture_surfaces = [sfc1, sfc2, sfc3, sfc4, sfc5]
+
+ def test_get_sofs(self):
+ rakes = np.array([-90.0, 0.0, 90.0])
+ sofs = ["N", "SS", "R"]
+ self.assertListEqual(self.sampler.get_sofs(rakes), sofs)
+
+ def test_sample_rupture_dimensions_no_uncertainty(self):
+ # Tests the sampling of the rupture dimensons in the case when the
+ # default scaling relation is used, i.e. no uncertainty on the scaling
+ eq1 = Earthquake("XYZ", 30.0, 40.0, 10.0, 7.0, lower_seismogenic_depth=20.0)
+ rakes = np.array([-90.0, 0.0, 0.0])
+ dips = np.array([45.0, 90.0, 45.0])
+
+ target_areas = np.array([1000.0, 1000.0, 1000.0])
+ target_lengths = np.array([35.35533906, 50.0, 35.35533906])
+ target_widths = np.array([28.28427125, 20.0, 28.28427125])
+ areas, lengths, widths = self.sampler.sample_rupture_dimensions(3, eq1, dips, rakes)
+ np.testing.assert_array_almost_equal(areas, target_areas)
+ np.testing.assert_array_almost_equal(lengths, target_lengths)
+ np.testing.assert_array_almost_equal(widths, target_widths)
+
+ def test_sample_rupture_dimensions_with_uncertainty(self):
+ # Tests the sampling of rupture dimensions in the case when a scaling relation
+ # is used that has uncertainty, i.e. Stafford (2014)
+ eq1 = Earthquake(
+ "XYZ",
+ 30.0,
+ 40.0,
+ 10.0,
+ 7.0,
+ mag_scaling_relation=Stafford2014(),
+ lower_seismogenic_depth=20.0,
+ )
+
+ # Test with controlled seed
+ nsamples = 3
+ rakes = np.array([-90.0, 0.0, 90.0])
+ dips = np.array([45.0, 90.0, 45.0])
+ np.random.seed(1000)
+ areas, lengths, widths = self.sampler.sample_rupture_dimensions(
+ nsamples, eq1, dips, rakes
+ )
+ np.testing.assert_array_almost_equal(
+ areas, np.array([1083.31109826, 667.07807089, 1371.79933828])
+ )
+
+ np.testing.assert_array_almost_equal(
+ lengths, np.array([38.90911277, 38.68736142, 48.50043073])
+ )
+
+ np.testing.assert_array_almost_equal(
+ widths, np.array([27.84209202, 17.24279057, 28.28427125])
+ )
+
+ nsamples = 1000
+ # Consider strike-slip events with vertical ruptures
+ rakes = np.zeros(nsamples)
+ dips = 90.0 + np.zeros(nsamples)
+ areas, lengths, widths = self.sampler.sample_rupture_dimensions(
+ nsamples, eq1, dips, rakes
+ )
+ self.assertTrue(np.all(widths <= 20.0))
+
+ def test_sample_hypocentre_locations(self):
+ # Tests the sampling of hypocentres - no clear target distribution can be defined
+ # to random sampling is used here
+ sofs = ["R", "SS", "N"]
+ n = len(sofs)
+ np.random.seed(1000)
+ hypo_pos = self.sampler.sample_hypocenter_position(n, sofs)
+ target_hypo_pos = np.array(
+ [[0.55321552, 0.54237704], [0.6759978, 0.58733529], [0.58199734, 0.8697838]]
+ )
+ np.testing.assert_array_almost_equal(hypo_pos, target_hypo_pos)
+
+ def test_get_distances(self):
+ # Tests the distance calculation for a set of pre-defined rupture surfaces
+ eq1 = Earthquake("XYZ", 30.0, 40.0, 10.0, 7.0, lower_seismogenic_depth=20.0)
+
+ target_lons = np.array([29.5, 30.5])
+ target_lats = np.array([40.0, 40.0])
+ distances, rhypo, repi = self.sampler.get_distances(
+ eq1, self.rupture_surfaces, target_lons, target_lats
+ )
+ target_distances = np.array(
+ [
+ [[22.0, 22.0, 22.0, 0.1], [22.0, 22.0, 22.0, 0.0]],
+ [[11.0, 11.0, 11.0, 0.1], [11.0, 11.0, 11.0, 0.0]],
+ [[00.0, 1.0, 0.0, 0.0], [00.0, 1.0, 0.0, 0.0]],
+ [[11.0, 11.0, -11.0, 0.0], [11.0, 11.0, -11.0, 0.1]],
+ [[22.0, 22.0, -22.0, 0.0], [22.0, 22.0, -22.0, 0.1]],
+ ]
+ )
+ np.testing.assert_array_almost_equal(np.around(distances, 1), target_distances)
+ np.testing.assert_array_almost_equal(rhypo, np.array([43.74830548, 43.74830548]))
+ np.testing.assert_array_almost_equal(repi, np.array([42.59007199, 42.59007199]))
+
+ def test_get_central_rupture(self):
+ # Tests the function to get the "most central" rupture. In this case it is the
+ # second rupture of the set of 5 defined in the setUp.
+ eq1 = Earthquake("XYZ", 30.0, 40.0, 10.0, 7.0, lower_seismogenic_depth=20.0)
+
+ target_lons = np.array([29.5, 30.5])
+ target_lats = np.array([40.0, 40.0])
+ distances, rhypo, repi = self.sampler.get_distances(
+ eq1, self.rupture_surfaces, target_lons, target_lats
+ )
+ cent_rupt, distances_min_loc, min_loc = self.sampler.get_central_rupture(
+ self.rupture_surfaces, distances, rhypo
+ )
+ self.assertEqual(min_loc, 1)
+ np.testing.assert_array_almost_equal(distances_min_loc, distances[1, :, :])
+
+ def test_sample_rupture_surface(self):
+ # Tests the sampling of the rupture surfaces
+ focal_mechanism = {
+ "nodal_plane_1": {"strike": 0.0, "dip": 90.0, "rake": -10.0},
+ "nodal_plane_2": {"strike": 90.0, "dip": 80.0, "rake": 170.0},
+ }
+ eq1 = Earthquake(
+ "XYZ",
+ 30.0,
+ 40.0,
+ 10.0,
+ 7.0,
+ lower_seismogenic_depth=20.0,
+ focal_mechanism=focal_mechanism,
+ )
+ nsamples = 1000
+ # Generate the surfaces
+ sample_surfaces, strikes, dips, rakes, hypo_pos = self.sampler.sample_rupture_surfaces(
+ nsamples, eq1
+ )
+ # Given that the distributions are random, we adopt a set of sanity checks on the
+ # resulting values
+ self.assertEqual(len(sample_surfaces), nsamples)
+ surface_strikes = []
+ surface_dips = []
+ for surface in sample_surfaces:
+ self.assertIsInstance(surface, PlanarSurface)
+ surface_strikes.append(surface.get_strike())
+ surface_dips.append(surface.get_dip())
+ surface_strikes = np.array(surface_strikes)
+ surface_dips = np.array(surface_dips)
+ # The strikes and dips are not perfectly distributed at the initial angles
+ # so very that they are close and that the number of each strike and each dip
+ # is within 400 and 600
+
+ # Strikes should be distributed at 0.0 and 90.0
+ idx1 = np.logical_and(surface_strikes >= -1.0, surface_strikes <= 1.0)
+ idx2 = np.logical_and(surface_strikes >= 89.0, surface_strikes <= 91.0)
+ self.assertTrue(np.sum(idx1) >= 400 and np.sum(idx1) <= 600)
+ self.assertTrue(np.sum(idx2) >= 400 and np.sum(idx2) <= 600)
+
+ # Dips should be distributed at 80.0 and 90.0
+ idx1 = np.logical_and(surface_dips >= 79.0, surface_dips <= 81.0)
+ idx2 = np.logical_and(surface_dips >= 89.0, surface_dips <= 91.0)
+ self.assertTrue(np.sum(idx1) >= 400 and np.sum(idx1) <= 600)
+ self.assertTrue(np.sum(idx2) >= 400 and np.sum(idx2) <= 600)
+
+ # For the rakes, take this from the outputs of the function
+ unique_rakes, rake_counter = np.unique(rakes, return_counts=True)
+ np.testing.assert_array_almost_equal(unique_rakes, np.array([-10.0, 170.0]))
+ # The instances of each mechanism should be approximately equal
+ self.assertTrue(rake_counter[0] >= 400 and rake_counter[0] <= 600)
+ self.assertTrue(rake_counter[1] >= 400 and rake_counter[1] <= 600)
+
+ def _check_points_equal(self, pnt1, pnt2):
+ """
+ Determine if two points are equal
+ """
+ for attr in ["longitude", "latitude", "depth"]:
+ self.assertAlmostEqual(getattr(pnt1, attr), getattr(pnt2, attr), 6)
+
+ def _check_point_set_equal(self, pnt_set1, pnt_set2):
+ """
+ Determine if two sets of points are equal
+ """
+ for pnt1, pnt2 in zip(pnt_set1, pnt_set2):
+ self._check_points_equal(pnt1, pnt2)
+
+ def test_complete_finite_rupture_generator_unknown_sites(self):
+ np.random.seed(1000)
+ sampler = FiniteRuptureSampler()
+ focal_mechanism = {
+ "nodal_plane_1": {"strike": 0.0, "dip": 90.0, "rake": 0.0},
+ "nodal_plane_2": {"strike": 90.0, "dip": 90.0, "rake": 180.0},
+ }
+ eq1 = Earthquake(
+ "XYZ",
+ 30.0,
+ 40.0,
+ 10.0,
+ 7.0,
+ lower_seismogenic_depth=20.0,
+ focal_mechanism=focal_mechanism,
+ )
+
+ # Get the central rupture
+ cent_rupt, distance_set = sampler.get_finite_rupture(
+ 1000, eq1, maximum_site_distance=100.0, site_spacing=0.1
+ )
+ self.assertIsInstance(cent_rupt, ParametricProbabilisticRupture)
+ self.assertAlmostEqual(cent_rupt.mag, 7.0)
+ self.assertAlmostEqual(cent_rupt.rake, 0.0)
+ self._check_points_equal(cent_rupt.hypocenter, eq1.hypocenter)
+ # Target plane defined by set of corner points
+ top_left = Point(30.0, 39.780900, 0.0)
+ top_right = Point(30.0, 40.230561, 0.0)
+ bottom_right = Point(30.0, 40.230561, 20.0)
+ bottom_left = Point(30.0, 39.780900, 20.0)
+ self._check_point_set_equal(
+ [
+ cent_rupt.surface.top_left,
+ cent_rupt.surface.top_right,
+ cent_rupt.surface.bottom_right,
+ cent_rupt.surface.bottom_left,
+ ],
+ [top_left, top_right, bottom_right, bottom_left],
+ )
+
+ # Check that the distance set has the right keys and the right shape
+ self.assertListEqual(
+ list(distance_set), ["lon", "lat", "rjb", "rrup", "rx", "ry0", "rhypo", "repi"]
+ )
+ self.assertTrue(distance_set["rjb"].shape[0], 475)
+
+ def test_complete_finite_rupture_generator_known_sites(self):
+
+ # Define target sites - should be 5 altogether
+ target_lats = np.arange(39.0, 41.05, 0.5)
+ target_lons = 30.0 * np.ones(target_lats.shape)
+ vs30 = 800.0 * np.ones(target_lats.shape)
+ site_dframe = pd.DataFrame(
+ {
+ "lon": target_lons,
+ "lat": target_lats,
+ "vs30": vs30,
+ "vs30measured": np.ones(target_lats.shape, dtype=bool),
+ }
+ )
+ site_model = SiteModel.from_dataframe(site_dframe)
+ # Setup the calculation
+ np.random.seed(1000)
+ sampler = FiniteRuptureSampler()
+ focal_mechanism = {
+ "nodal_plane_1": {"strike": 0.0, "dip": 90.0, "rake": 0.0},
+ "nodal_plane_2": {"strike": 90.0, "dip": 90.0, "rake": 180.0},
+ }
+ eq1 = Earthquake(
+ "XYZ",
+ 30.0,
+ 40.0,
+ 10.0,
+ 7.0,
+ lower_seismogenic_depth=20.0,
+ focal_mechanism=focal_mechanism,
+ )
+ # Get the central rupture
+ cent_rupt, distance_set = sampler.get_finite_rupture(1000, eq1, site_model=site_model)
+ self.assertIsInstance(cent_rupt, ParametricProbabilisticRupture)
+ # Check outputs
+ self.assertAlmostEqual(cent_rupt.mag, 7.0)
+ self.assertAlmostEqual(cent_rupt.rake, 0.0)
+ self._check_points_equal(cent_rupt.hypocenter, eq1.hypocenter)
+ top_left = Point(30.0, 39.709612, 0.0)
+ top_right = Point(30.0, 40.159273, 0.0)
+ bottom_right = Point(30.0, 40.159273, 20.0)
+ bottom_left = Point(30.0, 39.709612, 20.0)
+ self._check_point_set_equal(
+ [
+ cent_rupt.surface.top_left,
+ cent_rupt.surface.top_right,
+ cent_rupt.surface.bottom_right,
+ cent_rupt.surface.bottom_left,
+ ],
+ [top_left, top_right, bottom_right, bottom_left],
+ )
+
+ # Check that the distance set has the right keys and the right shape
+ self.assertListEqual(
+ list(distance_set), ["lon", "lat", "rjb", "rrup", "rx", "ry0", "rhypo", "repi"]
+ )
+ self.assertTrue(distance_set["rjb"].shape[0], 5)
--
GitLab