Commit 9d49fab8 authored by Graeme Weatherill's avatar Graeme Weatherill
Browse files

Implements Earthquake Object and Input validators

parent 76917315
......@@ -138,7 +138,7 @@ dmypy.json
.pyre/
# Vim
.swp
*.swp
# OSX
.DS_Store
"""
Contains the core classes to represent an earthquake
"""
import datetime
from typing import Tuple
from math import radians, sin, tan, sqrt
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
class DummyScalingRelation(object):
"""
A simple magnitude-area scaling relation to use as a placeholder until
updated by the synthetic rupture generator
"""
@staticmethod
def get_rupture_dimensions(
magnitude: float, rake: float = 0.0, dip: float = 90.0, aspect: float = 1.0,
thickness: float = 20.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:
A = 10.0 ^ (magnitude - 4.0)
Args:
magnitude:
Earthquake magnitude
rake:
Rake of the earthquake rupture in degrees (in the range -180 to 180)
dip:
Dip of the earthquake rupture in degrees (in the range 0.0 to 90.0)
aspect:
Aspect ratio of the rupture
thickness:
Seismogenic thickness of the crust (in km)
Returns:
Dimensions of the rupture as a tuple of rupture area (in km^2), rupture length
(in km) and downdip rupture width (in km)
"""
# PEER defined magnitude-area scaling relation returns area in km^2
area = 10.0 ** (magnitude - 4.0)
width = sqrt(area / aspect)
# dz is the vertical dimension of the rupture (in km)
dz = width * sin(radians(dip))
if dz > thickness:
# Rupture is wider than the crustal thickness, so rescale it to the
# maximum width and break the aspect ratio
width = thickness / sin(radians(dip))
length = area / width
return area, length, width
class Earthquake(object):
"""
Shakemap event class requires a minimum of an ID, longitude, latitude,
depth and magnitude. Other attributes can be input to control the rupture
orientation and mechanism.
Can input a rupture geometry directly or else the rupture geometry will be
built from the other properties input
Attributes:
id:
Unique ID for the earthquake
lon:
Earthquake longitude in decimal degrees (-180 to 180)
lat:
Earthquake latitude in decimal degrees (-90, to 90)
depth:
Hypocentral depth (in km)
mag:
Earthquake magnitude
strike:
Strike of the rupture (in degrees from north)
dip:
Dip of the rupture (in degrees from horizontal)
rake:
Rake of the earthquake rupture in decimal degrees according to the
Aki & Richards (1980) convention (-180 to 180)
aspect:
Aspect ratio (length / width) of the rupture
mag_scale_rel:
Magnitude scaling relation
lsd:
Lower seismogenic depth (in km)
usd:
Upper seismogenic depth (in km)
datetime.date date:
Date of the earthquake (year, month, day) as datetime.date object
datetime.time time:
Time of the earthquake (hour, minute, second) as datetime.time object
surface:
Rupture surface as instance of :class:`openquake.hazardlib.geo.BaseSurface`
hypo_pos:
Hypocentre location within the rupture plane as a tuple of
(fraction of along strike length, fraction of down-dip width)
"""
def __init__(
self,
id,
lon,
lat,
hypo_depth,
magnitude,
strike=0.0,
dip=90.0,
rake=0.0,
aspect=1.0,
mag_scaling_relation=None,
upper_seismogenic_depth=0.0,
lower_seismogenic_depth=1000.0,
surface=None,
hypocenter_position=(0.5, 0.5),
date=None,
time=None,
):
self.id = id
self.lon = valid.longitude(lon)
self.lat = valid.latitude(lat)
self.depth = valid.positive_float(hypo_depth, "hypocentre depth")
self.mag = magnitude
self.strike = valid.strike(strike)
self.dip = valid.dip(dip)
self.rake = valid.rake(rake)
self.aspect = valid.positive_float(aspect, "aspect ratio")
self.usd, self.lsd = valid.seismogenic_thickness(
upper_seismogenic_depth, lower_seismogenic_depth
)
self.mag_scale_rel = mag_scaling_relation or DummyScalingRelation()
# Date and time should be parsed as datetime.date and datetime.time
# objects if defined, otherwise none
assert isinstance(date, datetime.date) or date is None
self.date = date
assert isinstance(time, datetime.time) or time is None
self.time = time
assert isinstance(surface, BaseSurface) or surface is None
self.surface = surface
self.hypo_pos = valid.hypocenter_position(hypocenter_position)
self._rupture = None
def __repr__(self):
# Returns a summary string of the event, with datetime if specified
if self.date:
if self.time:
datetime_string = str(datetime.datetime.combine(self.date, self.time))
else:
datetime_string = str(self.date)
return "{:s} {:s} ({:.5f}E, {:.5f}N, {:.2f} km) M {:.2f}".format(
self.id, datetime_string, self.lon, self.lat, self.depth, self.mag
)
else:
return "{:s} ({:.5f}E, {:.5f}N, {:.2f} km) M {:.2f}".format(
self.id, self.lon, self.lat, self.depth, self.mag
)
@property
def rupture(self):
"""
If a rupture is provided then it is returned, otherwise it will build the rupture
from the available information
"""
if self._rupture:
return self._rupture
centroid = Point(self.lon, self.lat, self.depth)
if self.surface:
# Rupture surface has been input, so build the OpenQuake rupture object from
# existing properties
self._rupture = ParametricProbabilisticRupture(
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
def build_planar_surface(
centroid: Point,
strike: float,
dip: float,
length: float,
width: float,
lsd: float,
usd: float = 0.0,
hypo_loc: Tuple[float, float] = (0.5, 0.5)
):
"""
From a set of rupture properties returns a planar surface whose dimensions
are constrained by the seismogenic thickness of the crust
Args:
centroid:
Centroid of the rupture as instance of `class`:openquake.hazardlib.geo.Point`
length:
Rupture length (in km)
width:
Down-dip rupture width (in km)
Returns:
Rupture plane as instance of :class:`openquake.hazardlib.geo.PlanarSurface`
"""
rdip = radians(dip)
thickness = lsd - usd
# Determine whether the upper edge of the plane would be above the Earth's surface
updip_width = hypo_loc[1] * width
downdip_width = (1.0 - hypo_loc[1]) * width
updip_depth_change = updip_width * sin(rdip)
downdip_depth_change = downdip_width * sin(rdip)
if centroid.depth < updip_depth_change:
# This would move the rupture above the top surface so translate
# the rupture down until the upper depth is at the top surface
offset = updip_depth_change - centroid.depth
updip_depth_change = centroid.depth
downdip_depth_change += offset
# Now to address the case that the bottom edge exceeds the seismogenic
# thickness
if downdip_depth_change > (lsd - centroid.depth):
if (updip_depth_change + downdip_depth_change) > thickness:
# Determine excess width and translate rupture updip
offset = (centroid.depth + downdip_depth_change) - lsd
offset_area = length * (offset / sin(rdip))
rw_max = thickness / sin(rdip)
length += offset_area / rw_max
updip_depth_change = centroid.depth
downdip_depth_change = lsd - centroid.depth
else:
# This would move the rupture below the lower surface, so relocate it
offset = (centroid.depth + downdip_depth_change) - lsd
downdip_depth_change = lsd - centroid.depth
updip_depth_change += offset
if dip % 90.0:
updip_surface_length = updip_depth_change / tan(rdip)
downdip_surface_length = downdip_depth_change / tan(rdip)
else:
# Vertical rupture, so no change in surface distance
updip_surface_length = 0.0
downdip_surface_length = 0.0
# Now deal with strike parts
left_length = hypo_loc[0] * length
right_length = (1.0 - hypo_loc[0]) * length
# Build corner points
downdip_dir = (dip + 90.0) % 360
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
)
bottom_left = mid_left.point_at(
downdip_surface_length, downdip_depth_change, downdip_dir
)
bottom_right = mid_right.point_at(
downdip_surface_length, downdip_depth_change, downdip_dir
)
try:
surface = PlanarSurface.from_corner_points(
top_left, top_right, bottom_right, bottom_left
)
except Exception as e:
# If an exception is raised then something was wrong in
# the geometry. Return the user information to help debug
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))
]
for pnt in [top_left, top_right, bottom_right, bottom_left]:
extended_error_message.append(str(pnt))
extended_error_message.append(str(e))
raise ValueError("\n".join(extended_error_message))
return surface
"""
Defines a set of input validation methods to check physically correct or
consistent quantities
"""
from typing import Tuple, Optional, Union
from openquake.hazardlib.geo.nodalplane import NodalPlane
def longitude(lon: float) -> float:
"""
Verify that the longitude is between -180 and 180 degrees
"""
if lon < -180.0 or lon > 180.0:
raise ValueError("Longitude %.4f is not between -180 and 180" % lon)
return lon
def latitude(lat: float) -> float:
"""
Verify that the latitude is between -90 and 90 degrees
"""
if lat < -90.0 or lat > 90.0:
raise ValueError("Latitude %.4f is not between -90 and 90" % lat)
return lat
def positive_float(value: float, key: str) -> float:
"""
Verify the value is a positive float (or zero) and raise error otherwise
"""
if value < 0.0:
raise ValueError("%s must be positive float or zero, %.6f input" % (key, value))
return value
def strike(value: Optional[float]) -> float:
"""
Verify that the strike is within the range 0 to 360 degrees, or else return None if
unspecified
"""
if value is not None and (value < 0 or value >= 360.0):
raise ValueError("Strike %.2f not in the range 0 to 360 degrees" % value)
return value
def dip(value: Optional[float]) -> float:
"""
Verify that the dip is within the range 0 to 90 degrees, or else return None if unspecified
"""
if value is not None and (value < 0 or value > 90.0):
raise ValueError("Dip %.2f not in the range 0 to 90 degrees" % value)
return value
def rake(value: Optional[float]) -> float:
"""
Verify that the rake is within the range -180 to 180 degrees, according to the Aki &
Richards (1980) convention, or else return None if unspecified
"""
if value is not None and (value < -180.0 or value > 180.0):
raise ValueError("Rake %.2f not in the range -180 to 180 degrees" % value)
return value
def mechanism(
istrike: float, idip: float, irake: float
) -> Union[Tuple[Optional[float], Optional[float], Optional[float]], NodalPlane]:
"""
Verifies that a valid focal mechanism is defined. A valid focal mechanism requires a
strike, dip and rake value, which are limited to the range (0, 360), (0, 90) and
(-180, 180) respectively. The Aki & Richards (1980) definition of rake is applied.
Note that this only checks validity of the mechanism in terms of input value,
not in terms of the complete physical properties of the mechanism
"""
mechanism = (strike(istrike), dip(idip), rake(irake))
if None in mechanism:
# Not a valid mechanism, return tuple
return (istrike, idip, irake)
else:
# Is a valid mechanism, so return openquake.hazardlib.nodalplane.NodalPlane object
return NodalPlane(*mechanism)
def seismogenic_thickness(
upper_seismo_depth: float, lower_seismo_depth: float
) -> Tuple[float, float]:
"""
Verifies that a valid seismogenic thickness of the crust is defined
"""
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))
return usd, lsd
def hypocenter_position(hypo_pos: Tuple[float, float]) -> Tuple[float, float]:
"""
Verifies that a hypocenter position is valid within the range [0, 1] for both the
along-strike and down-dip cases
"""
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)
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
"""
Test suite for earthquake classes
"""
import unittest
import numpy as np
from datetime import date, time
from math import radians, sin, pi
from openquake.hazardlib.geo import Point
from shakyground2.earthquake import DummyScalingRelation, Earthquake
class DummyScalingRelationTestCase(unittest.TestCase):
"""
Simple tests for the dummy magnitude scaling relation
"""
def setUp(self):
self.msr = DummyScalingRelation()
def test_simple_case_within_thickness(self):
area, length, width = self.msr.get_rupture_dimensions(6.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 10.0)
self.assertAlmostEqual(width, 10.0)
self.assertAlmostEqual(length * width, area)
def test_aspect_greater_than_one(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, aspect=2.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length / width, 2.0)
self.assertAlmostEqual(length * width, area)
def test_rupture_rescale(self):
area, length, width = self.msr.get_rupture_dimensions(6.0, thickness=5.0)
self.assertAlmostEqual(area, 100.0)
self.assertAlmostEqual(length, 20.0)
self.assertAlmostEqual(width, 5.0)
self.assertAlmostEqual(length / width, 4.0)
self.assertAlmostEqual(length * width, area)
def test_dipping_rupture_rescaled(self):
area, length, width = self.msr.get_rupture_dimensions(7.0, dip=60.0, thickness=15.0)
self.assertAlmostEqual(area, 1000.0)
self.assertAlmostEqual(length * width, area)
self.assertAlmostEqual(width, 15.0 / sin(radians(60.0)))
class EarthquakeTestCase(unittest.TestCase):
def setUp(self):
self.lon = 35.0
self.lat = 40.0
self.depth = 10.0
self.mag = 6.0
# Centroid is along equator, so degrees offset should
# be proportional to ratio of width / Earth's circumference
self.target_width = (10.0 / (2.0 * pi * 6371)) * 360.0 # Assumes Earth radius 6371 km
def test_minimal_instantiation(self):
# Minimal instantiation - no datetime
eq0 = Earthquake("XYZ", self.lon, self.lat, self.depth, self.mag)
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")
def test_create_planar_surface_simple(self):
# Test the creation of a simple planar surface that requires no
# translation or re-scaling
centroid = Point(0.0, 0.0, 5.0)
mag = 6.0
strike = 90.0
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = DummyScalingRelation.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)
# Verify correct plane from top left and bottom right points
self.assertAlmostEqual(plane0.top_left.longitude, -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, self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.bottom_right.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
def test_create_planar_surface_translation_top(self):
# Test the creation of a simple planar surface requiring translation
# downward to avoid breaking the surface - vertical fault case
centroid = Point(0.0, 0.0, 3.0)
mag = 6.0
strike = 90.0
dip = 90.0
usd = 0.0
lsd = 20.0
area, length, width = DummyScalingRelation.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)
self.assertAlmostEqual(plane0.top_left.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.top_left.depth, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.longitude, self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.bottom_right.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
def test_create_planar_surface_translation_bottom(self):
# Test the creation of a simple planar surface requiring translation
# downward to avoid breaking the lower seismogenic surface - vertical fault case
centroid = Point(0.0, 0.0, 7.0)
mag = 6.0
strike = 90.0
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.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)
self.assertAlmostEqual(plane0.top_left.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.top_left.depth, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.longitude, self.target_width / 2.0, 7)
self.assertAlmostEqual(plane0.bottom_right.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
def test_create_planar_surface_rescale(self):
# Test the creation of a simple planar surface that exceeds the seismogenic
# thickness and required re-scaling - vertical fault case
centroid = Point(0.0, 0.0, 5.0)
mag = 7.0 # Area should now be 1000 km^2
strike = 90.0
dip = 90.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.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.latitude, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.depth, 10.0, 7)
def test_create_planar_surface_rescale_dipping(self):
# Test the creation of a simple planar surface that exceeds the seismogenic
# thickness and required re-scaling - dipping fault case
centroid = Point(0.0, 0.0, 5.0)
mag = 7.0 # Area should now be 1000 km^2
strike = 90.0
dip = 60.0
usd = 0.0
lsd = 10.0
area, length, width = DummyScalingRelation.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
self.assertAlmostEqual(plane0.top_left.longitude, -0.402398, 6)
self.assertAlmostEqual(plane0.top_left.latitude, 0.022483, 6)
self.assertAlmostEqual(plane0.top_left.depth, 0.0, 7)
self.assertAlmostEqual(plane0.bottom_right.longitude, 0.402398, 6)
self.assertAlmostEqual(plane0.bottom_right.latitude, -0.022483, 6)