Merge remote-tracking branch 'upstream/master'

.gitignore

@@ -138,7 +138,7 @@ dmypy.json
 .pyre/
 
 # Vim
-.swp
+*.swp
 
 # OSX
 .DS_Store

shakyground2/earthquake.py

+"""
+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

shakyground2/valid.py

+"""
+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

tests/earthquake_test.py

+"""
+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)