Skip to content

GitLab

Commit dce96698 authored by g-weatherill's avatar g-weatherill
Browse files

Merge remote-tracking branch 'upstream/master'

parents da3e4ed3 90456796
Hide whitespace changes
Inline Side-by-side
.gitlab-ci.yml
View file @ dce96698
image: python:3-buster
image: debian:buster-slim
variables:
PIP_CACHE_DIR: "$CI_PROJECT_DIR/.cache/pip"
cache:
paths:
- .cache/pip
- venv/
before_script:
- python3 -V
- apt-get update
- apt-get install -y python3 git python3-pip
- pip3 install virtualenv
- virtualenv venv
- source venv/bin/activate
- pip3 install .
- pip3 install pytest-cov
- pip3 install .[tests]
- pip3 install .
formatting:
script:
......
Makefile
View file @ dce96698
SOURCES=shakyground2 *.py
SOURCES=shakyground2 tests *.py
LENGTH=96
check: $(SOURCES)
......
setup.py
View file @ dce96698
......@@ -10,7 +10,10 @@ setup(
description="",
license="AGPLv3+",
extras_require={"tests": test_requirements},
install_requires=[],
install_requires=[
"openquake.engine",
"geopandas",
],
packages=find_packages(),
python_requires=">=3.7",
)
shakyground2/earthquake.py
View file @ dce96698
......@@ -204,7 +204,7 @@ class Earthquake(object):
width: float,
lsd: float,
usd: float = 0.0,
hypo_loc: Tuple[float, float] = (0.5, 0.5)
hypo_loc: Tuple[float, float] = (0.5, 0.5),
):
"""
From a set of rupture properties returns a planar surface whose dimensions
......@@ -266,12 +266,8 @@ class Earthquake(object):
updip_dir = (dip - 90.0) % 360
mid_left = centroid.point_at(left_length, 0.0, (strike + 180.0) % 360.0)
mid_right = centroid.point_at(right_length, 0.0, strike)
top_left = mid_left.point_at(
updip_surface_length, -updip_depth_change, updip_dir
)
top_right = mid_right.point_at(
updip_surface_length, -updip_depth_change, updip_dir
)
top_left = mid_left.point_at(updip_surface_length, -updip_depth_change, updip_dir)
top_right = mid_right.point_at(updip_surface_length, -updip_depth_change, updip_dir)
bottom_left = mid_left.point_at(
downdip_surface_length, downdip_depth_change, downdip_dir
)
......@@ -288,7 +284,7 @@ class Earthquake(object):
extended_error_message = [
"Rupture surface failed to build with the following properties:",
"Strike: %.2f, Dip: %.2f, Length: %.2f, Width: %.2f, Hypo. Pos: %s"
% (strike, dip, length, width, str(hypo_loc))
% (strike, dip, length, width, str(hypo_loc)),
]
for pnt in [top_left, top_right, bottom_right, bottom_left]:
extended_error_message.append(str(pnt))
......
shakyground2/magnitude_scaling_relations.py
View file @ dce96698
......@@ -12,10 +12,17 @@ class BaseScalingRelation(metaclass=abc.ABCMeta):
Abstract base class representing an implementation of a magnitude scaling relation
"""
@abc.abstractmethod
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
def get_rupture_dimensions(
self,
magnitude: float,
rake: float = 0.0,
dip: float = 90.0,
aspect: float = 1.0,
thickness: float = 20.0,
epsilon: float = 0.0,
) -> Tuple[float, float, float]:
"""
Args:
magnitude:
......@@ -41,9 +48,15 @@ class PEERScalingRelation(BaseScalingRelation):
updated by the synthetic rupture generator
"""
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
def get_rupture_dimensions(
self,
magnitude: float,
rake: float = 0.0,
dip: float = 90.0,
aspect: float = 1.0,
thickness: float = 20.0,
epsilon: float = 0.0,
) -> Tuple[float, float, float]:
"""
Returns the area (km^2), length (km) and width (km) of the rupture using the PEER
Magnitude-Area scaling relation and constrained by the crustal thickness
......@@ -80,6 +93,7 @@ class StrasserEtAl2010Interface(PEERScalingRelation):
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@staticmethod
def get_area(magnitude: float) -> float:
"""
......@@ -97,6 +111,7 @@ class StrasserEtAl2010Inslab(PEERScalingRelation):
Interface and Intraslab Subduction-zone Earthquakes with Moment Magnitude", Seismological
Research Letters, 81: 941 - 950, doi:10.1785/gssrl.81.6.941
"""
@staticmethod
def get_area(magnitude: float) -> float:
"""
......@@ -117,27 +132,39 @@ class Stafford2014(BaseScalingRelation):
# Model coefficients for the specific style of faulting
COEFFS = {
# Coefficients from Table 1
1: {"U": {"a0": -27.4922, "a1": 4.6656, "a2": -0.2033},
1: {
"U": {"a0": -27.4922, "a1": 4.6656, "a2": -0.2033},
"SS": {"a0": -30.8395, "a1": 5.4184, "a2": -0.3044},
"N": {"a0": -36.9770, "a1": 6.3070, "a2": -0.1696},
"R": {"a0": -35.8239, "a1": 5.0680, "a2": -0.0457}},
"R": {"a0": -35.8239, "a1": 5.0680, "a2": -0.0457},
},
# Coefficients from Table 2
2: {"U": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
2: {
"U": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
"SS": {"b0": -2.300, "b1": 0.7167, "sigma": 0.2337},
"N": {"b0": -4.1055, "b1": 1.0370, "sigma": 0.2509},
"R": {"b0": -3.8300, "b1": 0.9982, "sigma": 0.2285}},
"R": {"b0": -3.8300, "b1": 0.9982, "sigma": 0.2285},
},
# Coefficients from Table 3
3: {"U": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
3: {
"U": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
"SS": {"gamma0": -9.3137, "sigma": 0.3138, "rho": 0.3104},
"N": {"gamma0": -9.2483, "sigma": 0.3454, "rho": 0.4336},
"R": {"gamma0": -9.2749, "sigma": 0.2534, "rho": 0.1376}},
"R": {"gamma0": -9.2749, "sigma": 0.2534, "rho": 0.1376},
},
# Coefficients from Table 4
4: {"U": 0.7574, "SS": 0.7574, "N": 0.8490, "R": 0.7496}
4: {"U": 0.7574, "SS": 0.7574, "N": 0.8490, "R": 0.7496},
}
def get_rupture_dimensions(self, magnitude: float, rake: float = 0.0, dip: float = 90.0,
aspect: float = 1.0, thickness: float = 20.0,
epsilon: float = 0.0) -> Tuple[float, float, float]:
def get_rupture_dimensions(
self,
magnitude: float,
rake: float = 0.0,
dip: float = 90.0,
aspect: float = 1.0,
thickness: float = 20.0,
epsilon: float = 0.0,
) -> Tuple[float, float, float]:
"""
Gets the rupture dimensions from for the given magnitude subject to the physical
constaints.
......@@ -193,8 +220,9 @@ class Stafford2014(BaseScalingRelation):
# normal
return "N"
def get_rupture_width(self, magnitude: float, dip: float, sof: str, thickness: float) ->\
Tuple[float, float, float, float]:
def get_rupture_width(
self, magnitude: float, dip: float, sof: str, thickness: float
) -> Tuple[float, float, float, float]:
"""
Returns the parameters needed to define the censored rupture width
distribution defined in equations 8 to 14
......@@ -206,44 +234,45 @@ class Stafford2014(BaseScalingRelation):
"""
# Gets the probability of a full width rupture
rw_max = thickness / sin(radians(dip))
z_i = self.COEFFS[1][sof]["a0"] + \
self.COEFFS[1][sof]["a1"] * magnitude +\
self.COEFFS[1][sof]["a2"] * rw_max
z_i = (
self.COEFFS[1][sof]["a0"]
+ self.COEFFS[1][sof]["a1"] * magnitude
+ self.COEFFS[1][sof]["a2"] * rw_max
)
# Probability of rupturing full seismogenic thickness from logistic regression
p_i = 1.0 / (1.0 + exp(-z_i))
# Median width from equation 16
ln_rw = self.COEFFS[2][sof]["b0"] + \
self.COEFFS[2][sof]["b1"] * magnitude
ln_rw = self.COEFFS[2][sof]["b0"] + self.COEFFS[2][sof]["b1"] * magnitude
# Equations 18 - 20
phi_rw = (log(rw_max) - ln_rw) / self.COEFFS[2][sof]["sigma"]
phi_rw_ncdf = norm.cdf(phi_rw)
ln_rw_trunc = ln_rw - self.COEFFS[2][sof]["sigma"] *\
(norm.pdf(phi_rw) / phi_rw_ncdf)
ln_rw_trunc = ln_rw - self.COEFFS[2][sof]["sigma"] * (norm.pdf(phi_rw) / phi_rw_ncdf)
mean_rw = p_i * log(rw_max) + (1.0 - p_i) * ln_rw_trunc
# Equations 21 - 22
stddev_rw = self._get_rupture_width_sigma(self.COEFFS[2][sof]["sigma"],
phi_rw,
phi_rw_ncdf,
p_i)
stddev_rw = self._get_rupture_width_sigma(
self.COEFFS[2][sof]["sigma"], phi_rw, phi_rw_ncdf, p_i
)
return mean_rw, stddev_rw, rw_max, p_i
@staticmethod
def _get_rupture_width_sigma(sigma: float, phi_rw: float, phi_rw_ncdf: float,
p_i: float) -> float:
def _get_rupture_width_sigma(
sigma: float, phi_rw: float, phi_rw_ncdf: float, p_i: float
) -> float:
"""
Returns the variabiliy in the rupture width described by equations 21 and 22
"""
denom = sqrt(2.0 * pi) * phi_rw_ncdf
if phi_rw_ncdf >= 0.0:
elem1 = sqrt(pi / 2.) * (1.0 + chi2.cdf(phi_rw, 3))
elem1 = sqrt(pi / 2.0) * (1.0 + chi2.cdf(phi_rw, 3))
else:
elem1 = sqrt(pi / 2.) * (1.0 - chi2.cdf(phi_rw, 3))
elem1 = sqrt(pi / 2.0) * (1.0 - chi2.cdf(phi_rw, 3))
elem2 = exp(-(phi_rw ** 2)) / denom
sigma_truncated = sqrt(((sigma ** 2.) / denom) * (elem1 - elem2))
sigma_truncated = sqrt(((sigma ** 2.0) / denom) * (elem1 - elem2))
return (1.0 - p_i) * sigma_truncated
def get_rupture_area(self, magnitude: float, sof: str, rw_max: float,
sigma_lnrw: float) -> Tuple[float, float]:
def get_rupture_area(
self, magnitude: float, sof: str, rw_max: float, sigma_lnrw: float
) -> Tuple[float, float]:
"""
Returns the rupture area conditioned upon the maximum rupture width and style of
faulting
......@@ -258,15 +287,16 @@ class Stafford2014(BaseScalingRelation):
ln_ra: Natural logarithm of rupture area
sigma_ra: Standard deviation of the natural logarithm of rupture area
"""
mw_crit = (log(rw_max) - self.COEFFS[2][sof]["b0"]) /\
self.COEFFS[2][sof]["b1"]
ln_ra = self.COEFFS[3][sof]["gamma0"] + log(10.) * magnitude
mw_crit = (log(rw_max) - self.COEFFS[2][sof]["b0"]) / self.COEFFS[2][sof]["b1"]
ln_ra = self.COEFFS[3][sof]["gamma0"] + log(10.0) * magnitude
if magnitude > mw_crit:
# Equation 23
ln_ra -= ((log(10.) / 4.) * (magnitude - mw_crit))
ln_ra -= (log(10.0) / 4.0) * (magnitude - mw_crit)
# Get the sigma log rupture area (equation 28)
sigma = self.COEFFS[3][sof]["sigma"]
sigma_ra = sqrt(
(sigma ** 2.) + (sigma_lnrw ** 2.) +
(2.0 * self.COEFFS[3][sof]["rho"] * sigma_lnrw * sigma))
(sigma ** 2.0)
+ (sigma_lnrw ** 2.0)
+ (2.0 * self.COEFFS[3][sof]["rho"] * sigma_lnrw * sigma)
)
return ln_ra, sigma_ra
shakyground2/rupture_mechanism.py
View file @ dce96698
......@@ -17,6 +17,7 @@ class DIP_RANGES(Enum):
for each style-of-faulting is provided and they are assumed to be uniformly
distributed.
"""
R = (20.0, 45.0)
SS = (75.0, 91.0)
N = (45.0, 75.0)
......@@ -40,6 +41,7 @@ class RuptureMechanism(object):
of :class:`openquake.hazardlib.pmf.PMF`
"""
def __init__(self, mechanism: Optional[PMF] = None):
if mechanism:
assert isinstance(mechanism, PMF)
......@@ -76,8 +78,12 @@ class RuptureMechanism(object):
return strikes, dips, rakes
@classmethod
def from_strike_dip_rake(cls, strike: Optional[float] = None, dip: Optional[float] = None,
rake: Optional[float] = None) -> RuptureMechanism:
def from_strike_dip_rake(
cls,
strike: Optional[float] = None,
dip: Optional[float] = None,
rake: Optional[float] = None,
) -> RuptureMechanism:
"""
Constructs the rupture mechanism distribution from a simple strike, dip
and rake combination, permitting None values if undefined
......@@ -87,13 +93,14 @@ class RuptureMechanism(object):
dip: Dip of fault (in decimal degrees between 0 and 90)
rake: Rake of fault (in decimal degrees between -180 and 180)
"""
return cls(cls.build_mechanism_distribution(valid.strike(strike),
valid.dip(dip),
valid.rake(rake)))
return cls(
cls.build_mechanism_distribution(
valid.strike(strike), valid.dip(dip), valid.rake(rake)
)
)
@classmethod
def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict)\
-> RuptureMechanism:
def from_nodal_planes(cls, nodal_plane_1: Dict, nodal_plane_2: Dict) -> RuptureMechanism:
"""
Constructs the rupture mechanism distribution from either one single
valid rupture plane or two valud nodal planes
......@@ -108,13 +115,21 @@ class RuptureMechanism(object):
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))]))
return cls(
PMF(
[
(0.5, NodalPlane(strike_1, dip_1, rake_1)),
(0.5, NodalPlane(strike_2, dip_2, rake_2)),
]
)
)
@staticmethod
def build_mechanism_distribution(strike: Optional[float] = None,
dip: Optional[float] = None,
rake: Optional[float] = None) -> PMF:
def build_mechanism_distribution(
strike: Optional[float] = None,
dip: Optional[float] = None,
rake: Optional[float] = None,
) -> PMF:
"""
Builds a mechanism distribution from a partial (or complete) characterisation of the
rupture mechanism
......
shakyground2/site_model.py
View file @ dce96698
......@@ -23,8 +23,8 @@ SITE_PROPERTIES = {
"xvf": np.float64, # Distance (km) to the volcanic front
"backarc": np.bool, # Site is in the subduction backarc (True) or else
"region": np.int32, # Region to which the site belongs
"geology": (np.string_, 20) # Geological classification for the site
}
"geology": (np.string_, 20), # Geological classification for the site
}
# In some cases reasonable default values can be used for relevant ground motion models
......@@ -218,15 +218,11 @@ class SiteModel(object):
if key not in dframe.columns:
# Use the default
if SITE_PROPERTIES[key] in ((np.bytes_, 20),):
site_array[key] = cls._get_string_array(
nsites, value, SITE_PROPERTIES[key]
)
site_array[key] = cls._get_string_array(nsites, value, SITE_PROPERTIES[key])
elif SITE_PROPERTIES[key] == np.bool:
site_array[key] = cls._get_boolean_array(nsites, value)
else:
site_array[key] = value * np.ones(
nsites, dtype=SITE_PROPERTIES[key]
)
site_array[key] = value * np.ones(nsites, dtype=SITE_PROPERTIES[key])
return cls(site_array)
@staticmethod
......@@ -277,9 +273,7 @@ class SiteModel(object):
if japan:
c1 = 412.0 ** 2.0
c2 = 1360.0 ** 2.0
return np.exp(
(-5.23 / 2.0) * np.log((np.power(vs30, 2.0) + c1) / (c2 + c1))
)
return np.exp((-5.23 / 2.0) * np.log((np.power(vs30, 2.0) + c1) / (c2 + c1)))
else:
c1 = 571.0 ** 4.0
c2 = 1360.0 ** 4.0
......
shakyground2/valid.py
View file @ dce96698
......@@ -4,8 +4,7 @@ consistent quantities
"""
from typing import Tuple, Optional, Union
from openquake.hazardlib.geo.nodalplane import NodalPlane
from shakyground2.magnitude_scaling_relations import (BaseScalingRelation,
PEERScalingRelation)
from shakyground2.magnitude_scaling_relations import BaseScalingRelation, PEERScalingRelation
def longitude(lon: float) -> float:
......@@ -93,8 +92,10 @@ def seismogenic_thickness(
usd = positive_float(upper_seismo_depth, "upper seismogenic depth")
lsd = positive_float(lower_seismo_depth, "lower seismogenic depth")
if lsd < usd:
raise ValueError("Lower seismogenic depth %.2f km shallower than upper seismogenic "
"depth %.2f km" % (lsd, usd))
raise ValueError(
"Lower seismogenic depth %.2f km shallower than upper seismogenic "
"depth %.2f km" % (lsd, usd)
)
return usd, lsd
......@@ -105,8 +106,9 @@ def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
"""
along_strike, down_dip = hypo_pos
if along_strike < 0.0 or along_strike > 1.0:
raise ValueError("Along strike position %.3f should be in the range 0 to 1"
% along_strike)
raise ValueError(
"Along strike position %.3f should be in the range 0 to 1" % along_strike
)
if down_dip < 0.0 or down_dip > 1.0:
raise ValueError("Down dip position %.3f should be in the range 0 to 1" % down_dip)
return along_strike, down_dip
......@@ -121,6 +123,7 @@ def scaling_relation(msr: Optional[BaseScalingRelation]):
# If no relation is defined then use the default
return PEERScalingRelation()
if not isinstance(msr, BaseScalingRelation):
raise TypeError("Magnitude Scaling Relation %s not instance of BaseScalingRelation" %
str(msr))
raise TypeError(
"Magnitude Scaling Relation %s not instance of BaseScalingRelation" % str(msr)
)
return msr
tests/earthquake_test.py
View file @ dce96698
......@@ -11,7 +11,6 @@ from shakyground2.magnitude_scaling_relations import PEERScalingRelation
class EarthquakeTestCase(unittest.TestCase):
def setUp(self):
self.lon = 35.0
self.lat = 40.0
......@@ -27,10 +26,18 @@ class EarthquakeTestCase(unittest.TestCase):
self.assertEqual(str(eq0), "XYZ (35.00000E, 40.00000N, 10.00 km) M 6.00")
# Minimal instantiation - no datetime
eq0 = Earthquake("XYZ", self.lon, self.lat, self.depth, self.mag,
date=date(2010, 10, 6), time=time(11, 35, 45))
self.assertEqual(str(eq0),
"XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00")
eq0 = Earthquake(
"XYZ",
self.lon,
self.lat,
self.depth,
self.mag,
date=date(2010, 10, 6),
time=time(11, 35, 45),
)
self.assertEqual(
str(eq0), "XYZ 2010-10-06 11:35:45 (35.00000E, 40.00000N, 10.00 km) M 6.00"
)
def test_create_planar_surface_simple(self):
# Test the creation of a simple planar surface that requires no
......@@ -41,10 +48,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip,
length, width, lsd, usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
# Length should be approximately 10 km, but allow a slight precision difference
self.assertAlmostEqual(plane0.length, 10.0, 1)
self.assertAlmostEqual(plane0.width, 10.0, 5)
......@@ -66,10 +73,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
# Verify correct plane from top left and bottom right points
self.assertAlmostEqual(plane0.top_left.longitude, -self.target_width / 2.0, 7)
......@@ -88,10 +95,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
# Verify correct plane from top left and bottom right points
self.assertAlmostEqual(plane0.top_left.longitude, -self.target_width / 2.0, 7)
......@@ -110,17 +117,16 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
# Verify correct plane from top left and bottom right points
self.assertAlmostEqual(plane0.top_left.longitude, -10.0 * self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.top_left.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.top_left.depth, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.longitude,
10.0 * self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.bottom_right.longitude, 10.0 * self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.bottom_right.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
......@@ -133,10 +139,10 @@ class EarthquakeTestCase(unittest.TestCase):
dip = 60.0
usd = 0.0
lsd = 10.0
area, length, width = PEERScalingRelation().get_rupture_dimensions(mag, dip=dip,
thickness=lsd - usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
self.assertAlmostEqual(width, 10.0 / sin(radians(dip)))
# Verify correct plane from top left and bottom right points - calculated by hand
......@@ -154,34 +160,58 @@ class EarthquakeTestCase(unittest.TestCase):
strike = 0.0
dip = 90.0
# Create a simple surface
area, length, width = PEERScalingRelation().get_rupture_dimensions(self.mag, dip=dip,
thickness=lsd - usd)
area, length, width = PEERScalingRelation().get_rupture_dimensions(
self.mag, dip=dip, thickness=lsd - usd
)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width,
lsd, usd)
plane0 = Earthquake.build_planar_surface(centroid, strike, dip, length, width, lsd, usd)
# Instantiate earthquake object with known surface
eq0 = 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,
surface=plane0)
eq0 = 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,
surface=plane0,
)
# Generate rupture
self.assertAlmostEqual(eq0.rupture.mag, self.mag)
self.assertAlmostEqual(eq0.rupture.rake, 0.0)
np.testing.assert_array_almost_equal(plane0.corner_lons,
eq0.rupture.surface.corner_lons)
np.testing.assert_array_almost_equal(plane0.corner_lats,
eq0.rupture.surface.corner_lats)