Commit e61ff9d3 authored by Graeme Weatherill's avatar Graeme Weatherill
Browse files

Implements initial complete shakemap feature

parent 90456796
Pipeline #21206 canceled with stage
......@@ -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)
Hypocentral depth (in km)
Representation of the hypocenter as an instance of
Earthquake magnitude
......@@ -63,14 +67,15 @@ class Earthquake(object):
hypocenter_position=(0.5, 0.5),
......@@ -121,6 +126,7 @@ class Earthquake(object):
self.lon = valid.longitude(lon) = valid.latitude(lat)
self.depth = valid.positive_float(hypo_depth, "hypocentre depth")
self.hypocenter = Point(self.lon,, 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
# 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
Mechanism distribution as an instance of the
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)
# 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.lon,, 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.
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)
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
# 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]
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(
self._rupture = ParametricProbabilisticRupture(
self.mag, self.rake, None, centroid, surface, 1.0, None
return self._rupture
......@@ -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
yield prob, nodal_plane
def __len__(self):
# Lenth corresponds to the number of mechanisms in the distribution
return len(
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):
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
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(
(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"]),
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
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(
), "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))
def build_mechanism_distribution(
strike: Optional[float] = None,
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__(
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,
Earthquake as an instance of the :class:shakyground2.earthquake.Earthquake
Target sites used for calculation of the ground motion values, as instances of
Set of ground motion models and their respective weights as a dictionary
{"GMMs": [...], "weights": [...]}
Tectonic region to which the earthquake is assigned
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 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.
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
Exponent of the distance-dependent site weighing adopted by the :class:
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
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 in list(fle):
"Earthquake %s already in cache file %s - "
"Running the shakemaps will overwrite this"
% (, self.cache_file)
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"
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(
# 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"]:
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"]
site_ctxt.attrs["has_bbox"] = False
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
intensity_measure_types: List of intensity measures for which the shakemaps will
be calculated
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 in list(fle):
del fle[]
fle_eq = fle.create_group(
# 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)
mean, [stddev] = gmm.get_mean_and_stddevs(
self.sctx, self.rctx, self.dctx, input_imt, [const.StdDev.TOTAL]
except KeyError:
"Ground motion model %s not defined for intensity "
"measure type %s" % (str(gmm), intensity_measure_type)
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:
fle_eq, shakemaps, aggregated_means, aggregated_stddevs, shakemap_dtypes
return aggregated_means, aggregated_stddevs, shakemaps
def _cache_shakemap(
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
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
......@@ -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 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]
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
def shape(self):
return self.site_array.shape[0]
def dataframe(self):
# Return the site array as a dataframe
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
"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