Add local earthquake model support using pygmm

navarroc · navarroc · commit 4f3bc00e9017 · 2024-07-24T17:03:52.000-05:00
diff --git a/pyincore/models/hazard/__init__.py b/pyincore/models/hazard/__init__.py
@@ -6,6 +6,7 @@
 
 
 from pyincore.models.hazard.earthquake import Earthquake
+from pyincore.models.hazard.localearthquake import LocalEarthquake
 from pyincore.models.hazard.hazard import Hazard
 from pyincore.models.hazard.hazarddataset import HazardDataset
 from pyincore.models.hazard.flood import Flood
diff --git a/pyincore/models/hazard/localearthquake.py b/pyincore/models/hazard/localearthquake.py
@@ -0,0 +1,166 @@
+
+import pygmm
+import numpy as np
+import math
+
+from pyincore.models.hazard.earthquake import Earthquake
+from pyincore.utils.hazardutil import HazardUtil
+from pyincore.models.units import Units
+
+
+class LocalEarthquake(Earthquake):
+
+    def __init__(self, metadata):
+        super().__init__(metadata)
+
+        attenuations = metadata["attenuations"]
+
+        # Currently, we only support one attenuation
+        self.attenuation = list(attenuations.keys())[0]
+
+        self.eq_parameters = self.get_pygmm_parameters(metadata)
+        self.default_units = {"pga": "g", "pgv": "cm/s", "pgd": "cm", "sa": "g"}
+
+    def read_hazard_values(self, payload: list, hazard_service=None, **kwargs):
+        """Retrieve bulk earthquake hazard values either from the Hazard service or read it from local Dataset
+
+        Args:
+            payload (list):
+            hazard_service (obj): Hazard service.
+            kwargs (dict): Keyword arguments.
+        Returns:
+            obj: Hazard values.
+
+        """
+
+        print("read value from model: "+self.attenuation)
+        model_class = getattr(pygmm, self.attenuation)
+        # magnitude = self.magnitude
+        eq_parameters = self.eq_parameters
+
+        response = []
+
+        # TODO shouldn't need this
+        magnitude = eq_parameters["mag"]
+        for req in payload:
+            hazard_values = []
+            for index, req_demand_type in enumerate(req["demands"]):
+                req_demand = req_demand_type.lower()
+                req_demand_period = 0.0
+
+                if "sa" in req_demand:
+                    req_demand_parts = req_demand.split(" ")
+                    req_demand_period = req_demand_parts[0]
+                    req_demand = req_demand_parts[1]
+
+                latitude = float(req["loc"].split(",")[0])
+                longitude = float(req["loc"].split(",")[1])
+                # longitude, latitude
+
+                local_site = [latitude, longitude]
+
+                rupture_distance_ergo = HazardUtil.compute_rupture_distance(magnitude, self.depth, self.mechanism,
+                                                                            self.azimuth_angle, self.dip_angle,
+                                                                            local_site, self.source_site,)
+                joyner_boore_distance = HazardUtil.compute_joyner_boore_distance(magnitude, self.mechanism, self.depth,
+                                                                                 self.azimuth_angle, self.dip_angle,
+                                                                                 local_site, self.source_site)
+
+                z_tor = self.coseismic_rupture_depth
+                down_dip_width = self.get_downdip_rupture_width(magnitude, self.mechanism)
+
+                horizontal_distance = HazardUtil.compute_horizontal_distance(joyner_boore_distance, z_tor,
+                                                                             down_dip_width, self.dip_angle,
+                                                                             self.azimuth_angle, rupture_distance_ergo)
+
+                # TODO - this should come from a shapefile of site conditions maybe?
+                site_condition = self.default_site_condition
+                site_vs30 = 760
+
+                eq_parameters["dist_rup"] = rupture_distance_ergo
+                eq_parameters["dist_jb"] = joyner_boore_distance
+                eq_parameters["dist_x"] = horizontal_distance
+
+                scenario = pygmm.Scenario(**eq_parameters)
+                model = model_class(scenario)
+                if req_demand == "pga":
+                    raw_hazard_value = model.pga
+                elif req_demand == "pgv":
+                    raw_hazard_value = model.pgv
+                elif req_demand == "pgd":
+                    raw_hazard_value = model.pgd
+                elif req_demand == "sa":
+                    sa_values = model.spec_accels
+                    sa_index_array = np.where(model.periods == float(req_demand_period))
+                    if len(sa_index_array[0]) != 0:
+                        sa_index = sa_index_array[0][0]
+                    else:
+                        sa_index = self.find_nearest(model.periods, float(req_demand_period))
+
+                    sa_values = model.spec_accels
+                    raw_hazard_value = sa_values[sa_index]
+
+                if raw_hazard_value is None:
+                    converted_hazard_value = raw_hazard_value
+                else:
+                    converted_hazard_value = Units.convert_hazard(raw_hazard_value, 
+                                                                  original_demand_units=self.default_units[req_demand],
+                                                                  requested_demand_units=req["units"][index])
+
+                hazard_values.append(converted_hazard_value)
+
+            req.update({"hazardValues": hazard_values})
+            response.append(req)
+
+        return response
+
+    def find_nearest(self, period_array, period):
+        array = np.asarray(period_array)
+        idx = (np.abs(period_array - period)).argmin()
+        return idx
+
+    def get_pygmm_parameters(self, metadata):
+
+        eq_parameters = metadata["eqParameters"]
+        magnitude = float(eq_parameters["magnitude"])
+        event_type = eq_parameters["eventType"]
+
+        # TODO - need to handle the multiple attenuation models case
+        fault_type_map = eq_parameters["faultTypeMap"]
+        mechanism = fault_type_map[self.attenuation]
+
+        if mechanism == "Strike-Slip":
+            mechanism_pygmm = "SS"
+        else:
+            # TODO handle other cases
+            mechanism_pygmm = mechanism
+
+        dip_angle = float(eq_parameters["dipAngle"])
+        # self.azimuth_angle = eq_parameters["azimuthAngle"]
+
+        # TODO remove these later
+        self.mechanism = mechanism_pygmm
+        self.dip_angle = dip_angle
+
+        self.azimuth_angle = 130.0
+        self.default_site_condition = "soil"
+        self.default_site_vs30 = 760
+        self.source_site = [float(eq_parameters["srcLatitude"]), float(eq_parameters["srcLongitude"])]
+        self.depth = float(eq_parameters["depth"])
+        self.coseismic_rupture_depth = float(eq_parameters["coseismicRuptureDepth"])
+
+        # TODO all the eq parameters should be mapped into a dictionary that can be passed to the pygmm Scenario so
+        # we can map IN-CORE parameters to pygmm parameters (e.g. Ztor is depth_tor in pygmm)
+
+        eq_parameters = {"mag": magnitude, "site_cond": self.default_site_condition, "v_s30": self.default_site_vs30,
+                         "dip": dip_angle, "mechanism": mechanism_pygmm, "event_type": event_type,
+                         "depth_tor": 3.0, "depth_1_0": 0.0}
+
+        return eq_parameters
+
+    def get_downdip_rupture_width(self, magnitude, mechanism):
+        # TODO this should come from a dict of coefficients for different mechanisms
+        a = -0.76
+        b = 0.27
+
+        return math.pow(10, a + b * magnitude)
diff --git a/pyincore/utils/hazardutil.py b/pyincore/utils/hazardutil.py
@@ -0,0 +1,153 @@
+# Copyright (c) 2024 University of Illinois and others. All rights reserved.
+#
+# This program and the accompanying materials are made available under the
+# terms of the Mozilla Public License v2.0 which accompanies this distribution,
+# and is available at https://www.mozilla.org/en-US/MPL/2.0/
+
+import logging
+import sys
+import math
+
+logging.basicConfig(stream=sys.stderr, level=logging.INFO)
+
+
+class HazardUtil:
+
+    @staticmethod
+    def compute_rupture_distance(magnitude, depth, mechanism, azimuth_angle, dip_angle, local_site, source_site):
+        original_distances = HazardUtil.compute_original_distance(local_site, source_site, depth)
+
+        rupture_length = HazardUtil.get_subsurface_rupture_length(magnitude, mechanism)
+        rupture_width = HazardUtil.get_downdip_rupture_width(magnitude, mechanism)
+        transformed_distances = HazardUtil.compute_transformed_distances(azimuth_angle, dip_angle, original_distances)
+
+        if math.fabs(transformed_distances[0]) <= rupture_width / 2.0 and math.fabs(transformed_distances[1]) <= rupture_length / 2.0:
+            r_rup = math.fabs(transformed_distances[2])
+        elif (math.fabs(transformed_distances[0]) > rupture_width / 2.0 and math.fabs(
+                transformed_distances[1]) <= rupture_length / 2.0):
+            r_rup = math.sqrt(transformed_distances[2] * transformed_distances[2] + math.pow((math.fabs(
+                transformed_distances[0]) - rupture_width / 2.0), 2))
+        elif (math.fabs(transformed_distances[0]) <= rupture_width / 2.0 and math.fabs(
+                transformed_distances[1]) > rupture_length / 2.0):
+            r_rup = math.sqrt(transformed_distances[2] * transformed_distances[2] + math.pow((math.fabs(
+                transformed_distances[1]) - rupture_length / 2.0), 2))
+        else:
+            r_rup = math.sqrt(transformed_distances[2] * transformed_distances[2] + math.pow((math.fabs(
+                transformed_distances[0]) - rupture_width / 2.0), 2) + math.pow((math.fabs(transformed_distances[1]) -
+                                                                                 rupture_length / 2.0), 2))
+
+        if r_rup < 0.0:
+            print("r_rup is less than 0?")
+
+        return r_rup
+
+    @staticmethod
+    def compute_joyner_boore_distance(magnitude, mechanism, depth, azimuth_angle, dip_angle, local_site, source_site):
+        original_distances = HazardUtil.compute_original_distance(local_site, source_site, depth)
+        rupture_length = HazardUtil.get_subsurface_rupture_length(magnitude, mechanism)
+        rupture_width = HazardUtil.get_downdip_rupture_width(magnitude, mechanism)
+
+        transformed_distances = HazardUtil.compute_transformed_distances(azimuth_angle, dip_angle, original_distances)
+
+        dip_angle_radians = math.radians(dip_angle)
+        distance_x = transformed_distances[0]
+        distance_y = transformed_distances[1]
+        distance_z = transformed_distances[2]
+
+        x_proj = distance_x * math.cos(dip_angle_radians) + distance_z * math.sin(dip_angle_radians)
+
+        if (math.fabs(x_proj) <= rupture_width * math.cos(dip_angle_radians) / 2.0 and math.fabs(distance_y) <=
+                rupture_length / 2.0):
+            r_jb = 0.0
+        elif (math.fabs(x_proj) > rupture_width * math.cos(dip_angle_radians) / 2.0 and math.fabs(distance_y) <=
+              rupture_length / 2.0):
+            r_jb = math.fabs(x_proj) - rupture_width * math.cos(dip_angle_radians) / 2.0
+        elif (math.fabs(x_proj) <= rupture_width * math.cos(dip_angle_radians) / 2.0 and math.fabs(distance_y) >
+              rupture_length / 2.0):
+            r_jb = math.fabs(distance_y) - rupture_length / 2.0
+        else:
+            r_jb = math.sqrt(math.pow(math.fabs(x_proj) - rupture_width * math.cos(dip_angle_radians) / 2.0,
+                                      2) + math.pow(math.fabs(distance_y) - rupture_length / 2.0, 2))
+
+        return r_jb
+
+    @staticmethod
+    def compute_horizontal_distance(joyner_boore_distance, z_tor, down_dip_width, dip_angle, azimuth_angle, r_rup):
+        dip_angle_radians = math.radians(dip_angle)
+        azimuth_angle_radians = math.radians(azimuth_angle)
+
+        if dip_angle != 90.0:
+            if azimuth_angle >= 0 and azimuth_angle <= 180 and azimuth_angle != 90.0:
+                if joyner_boore_distance * math.fabs(math.tan(azimuth_angle_radians)) <= down_dip_width * math.cos(
+                        dip_angle_radians):
+                    rx = joyner_boore_distance * math.tan(azimuth_angle_radians)
+                else:
+                    rx = joyner_boore_distance * math.tan(azimuth_angle_radians) * math.cos(azimuth_angle_radians -
+                                                                                            math.asin((down_dip_width *
+                                                                                                       math.cos(dip_angle_radians) *
+                                                                 math.cos(azimuth_angle_radians) / joyner_boore_distance)))
+            elif azimuth_angle == 90.0:
+                if joyner_boore_distance > 0:
+                    rx = joyner_boore_distance + down_dip_width * math.cos(dip_angle_radians)
+                else:
+                    if r_rup < z_tor * (1 / math.cos(dip_angle_radians)):
+                        rx = math.sqrt(math.pow(r_rup, 2) - math.pow(z_tor, 2))
+                    else:
+                        rx = r_rup * (1 / math.sin(azimuth_angle_radians)) - z_tor * (1 / math.tan(dip_angle_radians))
+            else:
+                rx = joyner_boore_distance * math.sin(azimuth_angle_radians)
+        else:
+            rx = joyner_boore_distance * math.sin(azimuth_angle_radians)
+
+        return rx
+
+    @staticmethod
+    def compute_original_distance(local_site, source_site, depth):
+        site_latitude = local_site[0]
+        site_longitude = local_site[1]
+        source_latitude = source_site[0]
+        source_longitude = source_site[1]
+
+        original_distance = []
+        # 6373 is approximate radius of earth
+        original_distance.append((site_longitude - source_longitude) * math.pi * 6373.0 / 180.0)
+        original_distance.append((site_latitude - source_latitude) * math.pi * 6373.0 / 180.0)
+        original_distance.append(depth)
+
+        return original_distance
+
+    @staticmethod
+    def compute_transformed_distances(azimuth_angle, dip_angle, original_distances):
+        azimuth_angle_radians = math.radians(azimuth_angle)
+        dip_angle_radians = math.radians(dip_angle)
+
+        distance_x = (math.cos(dip_angle_radians) * math.cos(azimuth_angle_radians) * original_distances[0] -
+                      math.cos(dip_angle_radians) * math.sin(azimuth_angle_radians) * original_distances[1] -
+                      math.sin(dip_angle_radians) * original_distances[2])
+
+        distance_y = math.sin(azimuth_angle_radians) * original_distances[0] + math.cos(
+            azimuth_angle_radians) * original_distances[1]
+
+        distance_z = (math.sin(dip_angle_radians) * math.cos(azimuth_angle_radians) * original_distances[0] -
+                      math.sin(dip_angle_radians) * math.sin(azimuth_angle_radians) * original_distances[1] + math.cos(
+                    dip_angle_radians) * original_distances[2])
+
+        transformed_distances = [distance_x, distance_y, distance_z]
+
+        return transformed_distances
+
+    @staticmethod
+    def get_subsurface_rupture_length(magnitude, mechanism):
+        # TODO this should come from a dict of coefficients for different mechanisms
+        a = -2.57
+        b = 0.62
+
+        return math.pow(10, a + b * magnitude)
+
+    @staticmethod
+    def get_downdip_rupture_width(magnitude, mechanism):
+        # TODO this should come from a dict of coefficients for different mechanisms
+        a = -0.76
+        b = 0.27
+
+        return math.pow(10, a + b * magnitude)
diff --git a/tests/pyincore/analyses/localearthquake/test_localeq.py b/tests/pyincore/analyses/localearthquake/test_localeq.py
