Commit 37e7b0e2 authored by g-weatherill's avatar g-weatherill
Browse files

Changes abc.abstract property usage and stops forced calculation of MMI in shakemap_from_quakeml

parent 8d8b1e43
Pipeline #24419 passed with stage
in 14 minutes and 12 seconds
......@@ -10,7 +10,7 @@ shakemaps from the ground motion to macroseismic intensity. Includes:
import abc
import numpy as np
from typing import Dict, Tuple
from typing import Dict, Tuple, Set
from scipy.constants import g
from openquake.hazardlib.gsim.base import SitesContext, DistancesContext, RuptureContext
......@@ -23,14 +23,29 @@ class GMICE(metaclass=abc.ABCMeta):
macroseismic intensity
"""
# Distance attributes required by the GMICE (as a set)
REQUIRES_DISTANCES = abc.abstractproperty
@property
@abc.abstractmethod
def REQUIRES_DISTANCES(self) -> Set:
"""Shakemap source-to-site distances required by the GMICE"""
pass
@property
@abc.abstractmethod
def REQUIRES_SITES_PARAMETERS(self) -> Set:
"""Shakemap site properties required by the GMICE"""
pass
# Site attributes required by the GMICE (as a set)
REQUIRES_SITES_PARAMETERS = abc.abstractproperty
@property
@abc.abstractmethod
def REQUIRES_RUPTURE_PARAMETERS(self) -> Set:
"""Rupture attributes required by the GMICE"""
pass
# Rupture attributes required by the GMICE (as a set)
REQUIRES_RUPTURE_PARAMETERS = abc.abstractproperty
@property
@abc.abstractmethod
def COEFFS(self) -> Dict:
"""Coefficients of the GMICE as a dictionary organised by intensity measure type"""
pass
def __init__(self, sctx, rctx, dctx):
"""
......@@ -78,9 +93,6 @@ class GMICE(metaclass=abc.ABCMeta):
error propagation of the ground motion uncertainty
"""
# Coefficients of the GMICE
COEFFS = abc.abstractproperty
class WordenEtAl2012(GMICE):
"""
......@@ -95,14 +107,20 @@ class WordenEtAl2012(GMICE):
https://github.com/usgs/shakemap/blob/master/shakelib/gmice/wgrw12.py
"""
# No distance predictors required
REQUIRES_DISTANCES = {}
@property
def REQUIRES_DISTANCES(self) -> Set:
"""No distance predictors required"""
return {}
# No site predictors required
REQUIRES_SITES_PARAMETERS = {}
@property
def REQUIRES_SITES_PARAMETERS(self) -> Set:
"""No site predictors required"""
return {}
# No rupture predictors required
REQUIRES_RUPTURE_PARAMETERS = {}
@property
def REQUIRES_RUPTURE_PARAMETERS(self) -> Set:
"""No rupture predictors required"""
return {}
# Define the order of preference of the IMTs for conversion
ORDER_OF_PREFERENCE = ["PGV", "PGA", "SA(1.0)", "SA(0.3)", "SA(3.0)"]
......@@ -173,35 +191,51 @@ class WordenEtAl2012(GMICE):
] ** 2.0
return np.sqrt(stddevs)
# Coefficients for the bi-linear model according to Table 1
COEFFS = {
"PGA": {"c1": 1.78, "c2": 1.55, "c3": -1.6, "c4": 3.7, "t1": 1.57, "sigma_mmi": 0.73},
"PGV": {"c1": 3.78, "c2": 1.47, "c3": 2.89, "c4": 3.16, "t1": 0.53, "sigma_mmi": 0.65},
"SA(0.3)": {
"c1": 1.26,
"c2": 1.69,
"c3": -4.15,
"c4": 4.14,
"t1": 2.21,
"sigma_mmi": 0.84,
},
"SA(1.0)": {
"c1": 2.50,
"c2": 1.51,
"c3": 0.20,
"c4": 2.90,
"t1": 1.65,
"sigma_mmi": 0.80,
},
"SA(3.0)": {
"c1": 3.81,
"c2": 1.17,
"c3": 1.99,
"c4": 3.01,
"t1": 0.99,
"sigma_mmi": 0.95,
},
}
@property
def COEFFS(self) -> Dict:
"""Coefficients for the bi-linear model according to Table 1"""
return {
"PGA": {
"c1": 1.78,
"c2": 1.55,
"c3": -1.6,
"c4": 3.7,
"t1": 1.57,
"sigma_mmi": 0.73,
},
"PGV": {
"c1": 3.78,
"c2": 1.47,
"c3": 2.89,
"c4": 3.16,
"t1": 0.53,
"sigma_mmi": 0.65,
},
"SA(0.3)": {
"c1": 1.26,
"c2": 1.69,
"c3": -4.15,
"c4": 4.14,
"t1": 2.21,
"sigma_mmi": 0.84,
},
"SA(1.0)": {
"c1": 2.50,
"c2": 1.51,
"c3": 0.20,
"c4": 2.90,
"t1": 1.65,
"sigma_mmi": 0.80,
},
"SA(3.0)": {
"c1": 3.81,
"c2": 1.17,
"c3": 1.99,
"c4": 3.01,
"t1": 0.99,
"sigma_mmi": 0.95,
},
}
# Coefficients for the low intensity extension from
# https://github.com/usgs/shakemap/blob/master/shakelib/gmice/wgrw12.py
......@@ -228,14 +262,18 @@ class WordenEtAl2012MRDependent(WordenEtAl2012):
"""
# Requires rupture distance
REQUIRES_DISTANCES = {
"rrup",
}
@property
def REQUIRES_DISTANCES(self) -> Set:
return {
"rrup",
}
# Requires magnitude
REQUIRES_RUPTURE_PARAMETERS = {
"mag",
}
@property
def REQUIRES_RUPTURE_PARAMETERS(self) -> Set:
return {
"mag",
}
def get_mean_intensity_and_stddev(
self, imt: str, log_gmvs: np.ndarray, gmv_total_stddev: np.ndarray, clip: bool = True
......@@ -259,65 +297,68 @@ class WordenEtAl2012MRDependent(WordenEtAl2012):
mmi = np.clip(mmi, 1.0, 10.0)
return mmi, stddev
# Coefficients for the magnitude- and distance-dependent bi-linear model according to
# Tables 1 and 2
COEFFS = {
"PGA": {
"c1": 1.78,
"c2": 1.55,
"c3": -1.6,
"c4": 3.7,
"t1": 1.57,
"c5": -0.91,
"c6": 1.02,
"c7": -0.17,
"sigma_mmi": 0.66,
},
"PGV": {
"c1": 3.78,
"c2": 1.47,
"c3": 2.89,
"c4": 3.16,
"t1": 0.53,
"c5": 0.90,
"c6": 0.0,
"c7": -0.18,
"sigma_mmi": 0.63,
},
"SA(0.3)": {
"c1": 1.26,
"c2": 1.69,
"c3": -4.15,
"c4": 4.14,
"t1": 2.21,
"c5": -1.05,
"c6": 0.60,
"c7": 0.00,
"sigma_mmi": 0.82,
},
"SA(1.0)": {
"c1": 2.50,
"c2": 1.51,
"c3": 0.20,
"c4": 2.90,
"t1": 1.65,
"c5": 2.27,
"c6": -0.49,
"c7": -0.29,
"sigma_mmi": 0.75,
},
"SA(3.0)": {
"c1": 3.81,
"c2": 1.17,
"c3": 1.99,
"c4": 3.01,
"t1": 0.99,
"c5": 1.91,
"c6": -0.57,
"c7": -0.21,
"sigma_mmi": 0.89,
},
}
@property
def COEFFS(self) -> Dict:
"""Coefficients for the magnitude- and distance-dependent bi-linear model according to
Tables 1 and 2
"""
return {
"PGA": {
"c1": 1.78,
"c2": 1.55,
"c3": -1.6,
"c4": 3.7,
"t1": 1.57,
"c5": -0.91,
"c6": 1.02,
"c7": -0.17,
"sigma_mmi": 0.66,
},
"PGV": {
"c1": 3.78,
"c2": 1.47,
"c3": 2.89,
"c4": 3.16,
"t1": 0.53,
"c5": 0.90,
"c6": 0.0,
"c7": -0.18,
"sigma_mmi": 0.63,
},
"SA(0.3)": {
"c1": 1.26,
"c2": 1.69,
"c3": -4.15,
"c4": 4.14,
"t1": 2.21,
"c5": -1.05,
"c6": 0.60,
"c7": 0.00,
"sigma_mmi": 0.82,
},
"SA(1.0)": {
"c1": 2.50,
"c2": 1.51,
"c3": 0.20,
"c4": 2.90,
"t1": 1.65,
"c5": 2.27,
"c6": -0.49,
"c7": -0.29,
"sigma_mmi": 0.75,
},
"SA(3.0)": {
"c1": 3.81,
"c2": 1.17,
"c3": 1.99,
"c4": 3.01,
"t1": 0.99,
"c5": 1.91,
"c6": -0.57,
"c7": -0.21,
"sigma_mmi": 0.89,
},
}
# Coefficients for the low intensity extension from
# https://github.com/usgs/shakemap/blob/master/shakelib/gmice/wgrw12.py
......@@ -345,14 +386,20 @@ class AtkinsonKaka2007(GMICE):
the approach of the USGS and apply these coefficients to SA(3.0) rather than SA(2.0)
"""
# No distance predictors required
REQUIRES_DISTANCES = {}
@property
def REQUIRES_DISTANCES(self) -> Set:
"""No distance predictors required"""
return {}
# No site predictors required
REQUIRES_SITES_PARAMETERS = {}
@property
def REQUIRES_SITES_PARAMETERS(self) -> Set:
"""No site predictors required"""
return {}
# No rupture predictors required
REQUIRES_RUPTURE_PARAMETERS = {}
@property
def REQUIRES_RUPTURE_PARAMETERS(self) -> Set:
"""No rupture predictors required"""
return {}
# Define the order of preference of the IMTs for conversion
ORDER_OF_PREFERENCE = ["PGV", "PGA", "SA(1.0)", "SA(3.0)", "SA(0.3)"]
......@@ -416,34 +463,51 @@ class AtkinsonKaka2007(GMICE):
] ** 2.0
return np.sqrt(stddevs)
COEFFS = {
"PGA": {"c1": 2.65, "c2": 1.39, "c3": -1.91, "c4": 4.09, "t1": 1.69, "sigma_mmi": 1.01},
"PGV": {"c1": 4.37, "c2": 1.32, "c3": 3.54, "c4": 3.03, "t1": 0.48, "sigma_mmi": 0.8},
"SA(0.3)": {
"c1": 2.40,
"c2": 1.36,
"c3": -1.83,
"c4": 3.56,
"t1": 1.92,
"sigma_mmi": 0.88,
},
"SA(1.0)": {
"c1": 3.23,
"c2": 1.18,
"c3": 0.57,
"c4": 2.95,
"t1": 1.50,
"sigma_mmi": 0.84,
},
"SA(3.0)": {
"c1": 3.72,
"c2": 1.29,
"c3": 1.99,
"c4": 3.0,
"t1": 1.00,
"sigma_mmi": 0.86,
},
}
@property
def COEFFS(self) -> Dict:
"""Coefficients of the GMICE as defined in Table 4"""
return {
"PGA": {
"c1": 2.65,
"c2": 1.39,
"c3": -1.91,
"c4": 4.09,
"t1": 1.69,
"sigma_mmi": 1.01,
},
"PGV": {
"c1": 4.37,
"c2": 1.32,
"c3": 3.54,
"c4": 3.03,
"t1": 0.48,
"sigma_mmi": 0.8,
},
"SA(0.3)": {
"c1": 2.40,
"c2": 1.36,
"c3": -1.83,
"c4": 3.56,
"t1": 1.92,
"sigma_mmi": 0.88,
},
"SA(1.0)": {
"c1": 3.23,
"c2": 1.18,
"c3": 0.57,
"c4": 2.95,
"t1": 1.50,
"sigma_mmi": 0.84,
},
"SA(3.0)": {
"c1": 3.72,
"c2": 1.29,
"c3": 1.99,
"c4": 3.0,
"t1": 1.00,
"sigma_mmi": 0.86,
},
}
class AtkinsonKaka2007MRDependent(AtkinsonKaka2007):
......@@ -456,18 +520,19 @@ class AtkinsonKaka2007MRDependent(AtkinsonKaka2007):
Society of America, 97(2), pp 497 - 510, doi: 10.1785/0120060154
"""
# No distance predictors required
REQUIRES_DISTANCES = {
"rrup",
}
# No site predictors required
REQUIRES_SITES_PARAMETERS = {}
@property
def REQUIRES_DISTANCES(self) -> Set:
"""Requires rupture distance"""
return {
"rrup",
}
# No rupture predictors required
REQUIRES_RUPTURE_PARAMETERS = {
"mag",
}
@property
def REQUIRES_RUPTURE_PARAMETERS(self) -> Set:
"""Requires magnitude"""
return {
"mag",
}
# Define the order of preference of the IMTs for conversion
ORDER_OF_PREFERENCE = ["SA(3.0)", "SA(1.0)", "PGV", "SA(0.3)", "PGA"]
......@@ -495,60 +560,63 @@ class AtkinsonKaka2007MRDependent(AtkinsonKaka2007):
mmi = np.clip(mmi, 1.0, 10.0)
return mmi, stddev
COEFFS = {
"PGA": {
"c1": 2.65,
"c2": 1.39,
"c3": -1.91,
"c4": 4.09,
"t1": 1.69,
"c5": -1.96,
"c6": 0.02,
"c7": 0.98,
"sigma_mmi": 0.89,
},
"PGV": {
"c1": 4.37,
"c2": 1.32,
"c3": 3.54,
"c4": 3.03,
"t1": 0.48,
"c5": 0.47,
"c6": -0.19,
"c7": 0.26,
"sigma_mmi": 0.76,
},
"SA(0.3)": {
"c1": 2.40,
"c2": 1.36,
"c3": -1.83,
"c4": 3.56,
"t1": 1.92,
"c5": -0.11,
"c6": -0.20,
"c7": 0.64,
"sigma_mmi": 0.79,
},
"SA(1.0)": {
"c1": 3.23,
"c2": 1.18,
"c3": 0.57,
"c4": 2.95,
"t1": 1.50,
"c5": 1.92,
"c6": -0.39,
"c7": 0.04,
"sigma_mmi": 0.73,
},
"SA(3.0)": {
"c1": 3.72,
"c2": 1.29,
"c3": 1.99,
"c4": 3.0,
"t1": 1.00,
"c5": 2.24,
"c6": -0.33,
"c7": -0.31,
"sigma_mmi": 0.72,
},
}
@property
def COEFFS(self) -> Dict:
"""Coefficients of the GMICE as defined in Tables 4 and 5"""
return {
"PGA": {
"c1": 2.65,
"c2": 1.39,
"c3": -1.91,
"c4": 4.09,
"t1": 1.69,
"c5": -1.96,
"c6": 0.02,
"c7": 0.98,
"sigma_mmi": 0.89,
},
"PGV": {
"c1": 4.37,
"c2": 1.32,
"c3": 3.54,
"c4": 3.03,
"t1": 0.48,
"c5": 0.47,
"c6": -0.19,
"c7": 0.26,
"sigma_mmi": 0.76,
},
"SA(0.3)": {
"c1": 2.40,
"c2": 1.36,
"c3": -1.83,
"c4": 3.56,
"t1": 1.92,
"c5": -0.11,
"c6": -0.20,
"c7": 0.64,
"sigma_mmi": 0.79,
},
"SA(1.0)": {
"c1": 3.23,
"c2": 1.18,
"c3": 0.57,
"c4": 2.95,
"t1": 1.50,
"c5": 1.92,
"c6": -0.39,
"c7": 0.04,
"sigma_mmi": 0.73,
},
"SA(3.0)": {
"c1": 3.72,
"c2": 1.29,
"c3": 1.99,
"c4": 3.0,
"t1": 1.00,
"c5": 2.24,
"c6": -0.33,
"c7": -0.31,
"sigma_mmi": 0.72,
},
}
......@@ -125,7 +125,7 @@ def get_mmi_shakemaps(
def shakemaps_from_quakeml(
event_id: str,
imts: List = ["PGA", "SA(0.3)", "SA(1.0)"],
imts: List = ["PGA", "SA(0.3)", "SA(1.0)", "MMI"],
config: Optional[Dict] = None,
export_folder: Optional[str] = None,
cache_file: Optional[str] = None,
......@@ -213,9 +213,7 @@ def shakemaps_from_quakeml(
results.shakemaps = {"mean": {}, "stddevs": {}}
results.contours = {"mean": {}, "stddevs": {}}
results.statistics = {"mean": {}}
for imt in imts + [
"MMI",
]:
for imt in imts:
if imt == "MMI":
# Get MMI and its standard deviation from the mean and stddev shakemaps
mmi, stddev_mmi = get_mmi_shakemaps(
......
......@@ -110,20 +110,14 @@ class GeofonWorkflowTestCase(unittest.TestCase):
"""
def setUp(self):
self.imts = ["PGA", "SA(1.0)"]
self.imts = ["PGA", "SA(1.0)", "MMI"]
def test_complete_workflow_no_export_default_config(self):
event_id = os.path.join(DATA_PATH, "gfz2021eksc.xml")
results = shakemaps_from_quakeml(event_id, self.imts)
self.assertIsInstance(results, ShakemapWorkflowResult)
self.assertListEqual(list(results.shakemaps), ["mean", "stddevs"])
self.assertListEqual(
list(results.shakemaps["mean"]),
self.imts
+ [
"MMI",
],
)
self.assertListEqual(list(results.shakemaps["mean"]), self.imts)
for imt in self.imts:
for maptype in results.shakemaps:
......@@ -134,17 +128,9 @@ class GeofonWorkflowTestCase(unittest.TestCase):
event_id = os.path.join(DATA_PATH, "gfz2021eksc.xml")
results = shakemaps_from_quakeml(event_id, self.imts, config={"spcx": 0.1, "spcy": 0.1})
self.assertListEqual(list(results.shakemaps), ["mean", "stddevs"])
self.assertListEqual(
list(results.shakemaps["mean"]),
self.imts
+ [
"MMI",
],
)
self.assertListEqual(list(results.shakemaps[