Add ASK14 aftershock GMPE

Author: jinyan1214

modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py

@@ -54,18 +54,21 @@
     'SA': [
         'Chiou & Youngs (2014)',
         'Abrahamson, Silva & Kamai (2014)',
+        'Abrahamson, Silva & Kamai (2014) Aftershock',
         'Boore, Stewart, Seyhan & Atkinson (2014)',
         'Campbell & Bozorgnia (2014)',
     ],
     'PGA': [
         'Chiou & Youngs (2014)',
         'Abrahamson, Silva & Kamai (2014)',
+        'Abrahamson, Silva & Kamai (2014) Aftershock',
         'Boore, Stewart, Seyhan & Atkinson (2014)',
         'Campbell & Bozorgnia (2014)',
     ],
     'PGV': [
         'Chiou & Youngs (2014)',
         'Abrahamson, Silva & Kamai (2014)',
+        'Abrahamson, Silva & Kamai (2014) Aftershock',
         'Boore, Stewart, Seyhan & Atkinson (2014)',
         'Campbell & Bozorgnia (2014)',
     ],
@@ -136,11 +139,16 @@ def __init__(
         gmpe_weights_dict=dict(),  # noqa: B006, C408
         im_type=None,
         site_info=dict(),  # noqa: B006, C408
+        mainshock=None
     ):
         # basic set-ups
         self.set_im_gmpe(im_dict, gmpe_dict, gmpe_weights_dict)
         self.set_im_type(im_type)
         self.set_sites(site_info)
+        if mainshock is not None:
+            self.mainshock = mainshock.copy() # single row gdf containing mainshock rupture surface
+        else:
+            self.mainshock = None
         # self.set_source(source_info)
 
     def set_source(self, source_info):  # noqa: D102
@@ -155,6 +163,7 @@ def set_source(self, source_info):  # noqa: D102
                 or 'Abrahamson, Silva & Kamai (2014)' in gmpe_list
                 or 'Boore, Stewart, Seyhan & Atkinson (2014)' in gmpe_list
                 or 'Campbell & Bozorgnia (2014)' in gmpe_list
+                or 'Abrahamson, Silva & Kamai (2014) Aftershock' in gmpe_list
             ):
                 source_index = source_info.get('SourceIndex', None)
                 rupture_index = source_info.get('RuptureIndex', None)
@@ -163,6 +172,14 @@ def set_source(self, source_info):  # noqa: D102
                     self.erf, source_index, rupture_index, self.site_info
                 )
                 # self.timeGetRuptureInfo += time.process_time_ns() - start
+            if 'Abrahamson, Silva & Kamai (2014) Aftershock' in gmpe_list:
+                source_index = source_info.get('SourceIndex', None)
+                rupture_index = source_info.get('RuptureIndex', None)
+                crJB = get_rupture_info_ASK2014_aftershock(  # noqa: F405
+                    self.erf, source_index, rupture_index, self.mainshock
+                )
+                site_rup_dict['crJB'] = crJB
+
         elif source_info['Type'] == 'PointSource':
             if (
                 'Chiou & Youngs (2014)' in gmpe_list
@@ -456,7 +473,7 @@ def get_im_from_local(  # noqa: C901, D102
                         eq_magnitude, self.site_rup_dict, cur_site, im_info
                     )
                     # self.timeGetIM += time.process_time_ns() - start
-                elif cur_gmpe == 'Abrahamson, Silva & Kamai (2014)':
+                elif cur_gmpe == 'Abrahamson, Silva & Kamai (2014)' or cur_gmpe == 'Abrahamson, Silva & Kamai (2014) Aftershock':
                     # start = time.process_time_ns()
                     tmpResult = self.ASK.get_IM(  # noqa: N806
                         eq_magnitude, self.site_rup_dict, cur_site, im_info
@@ -793,13 +810,17 @@ def compute_im(  # noqa: C901, D103
                     sys.exit(
                         'SA is used in hazard downsampling but not defined in the intensity measure tab'
                     )
-                elif ho_period in im_info['SA'].get('Periods'):
-                    pass
+                elif im_info.get('Periods', None) is not None:
+                    if ho_period in im_info['Periods']:
+                        pass
+                elif im_info.get('SA', None) is not None:
+                    if ho_period in im_info['SA'].get('Periods'):
+                        pass
                 else:
                     tmp_periods = im_info['SA']['Periods'] + [ho_period]
                     tmp_periods.sort()
                     im_info['SA']['Periods'] = tmp_periods
-            elif ho_period in im_info['SA'].get('Periods'):
+            elif ho_period in im_info.get('Periods'):
                 pass
             else:
                 tmp_periods = im_info['SA']['Periods'] + [ho_period]

modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py

@@ -46,6 +46,8 @@
 import sys
 import psutil
 import GlobalVariable
+import shapely
+import geopandas as gpd
 
 if 'stampede2' not in socket.gethostname():
     import GlobalVariable
@@ -395,6 +397,25 @@ def get_rupture_distance(erf, source_index, rupture_index, lat, lon):  # noqa: D
 
     return distToRupture
 
+def get_rupture_info_ASK2014_aftershock(
+        erf, source_index, rupture_indx, mainshock):
+    rupSource = erf.getSource(source_index)  # noqa: N806
+    rupList = rupSource.getRuptureList()  # noqa: N806
+    rupSurface = rupList.get(rupture_indx).getRuptureSurface()
+    rupSurface_perimeter = rupSurface.getPerimeter()
+    coords = []
+    for i in range(rupSurface_perimeter.size()):
+        loc = rupSurface_perimeter.get(i)
+        coords.append((loc.getLongitude(), loc.getLatitude()))
+    if len(coords) == 1:
+        rup_polygon = shapely.geometry.Point(coords[0])
+    else:
+        rup_polygon = shapely.geometry.Polygon(coords)
+    rup_gdf = gpd.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[rup_polygon])
+    rup_gdf = rup_gdf.to_crs(epsg=6417)
+    centroid = rup_gdf.geometry.centroid.iloc[0]
+    return mainshock.geometry.iloc[0].distance(centroid) /1000 # in meters
+
 
 def get_rupture_info_CY2014(erf, source_index, rupture_index, siteList):  # noqa: N802, N803, D103
     rupSource = erf.getSource(source_index)  # noqa: N806

modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py

@@ -317,7 +317,7 @@ class abrahamson_silva_kamai_2014:  # noqa: D101
     timeCalc = 0  # noqa: N815
     supportedImt = None  # noqa: N815
 
-    def __init__(self):
+    def __init__(self, aftershock=False):
         self.coeff = pd.read_csv(
             os.path.join(os.path.dirname(__file__), 'data', 'ASK14.csv')  # noqa: PTH118, PTH120
         )
@@ -343,6 +343,7 @@ def __init__(self):
         self.H2 = 1.5
         self.H3 = -0.75
         self.PHI_AMP_SQ = 0.16
+        self.aftershock = aftershock
 
     def setIMT(self, imt):  # noqa: N802, D102
         if imt not in self.supportedImt:
@@ -442,6 +443,7 @@ def calcValues(  # noqa: C901, N802, D102
         vsInferred,  # noqa: N803
         z1p0,
         style,
+        crJB=None
     ):
         if Mw > 5:  # noqa: PLR2004
             c4mag = self.C4
@@ -537,8 +539,16 @@ def calcValues(  # noqa: C901, N802, D102
             )
         else:
             f5 = (self.a10 + self.b * self.N) * np.log(vs30s / self.Vlin)
-        # total model (no aftershock f11) -- Equation 1
-        mean = f1 + f78 + f5 + f4 + f6 + f10
+
+        # Aftershock term -- Equation 20
+        if self.aftershock:
+            if crJB is None:
+                raise ValueError('crJB must be provided for aftershock calculations')
+            f11 = self.a14 * np.clip(1 - (crJB - 5)/10, 0, 1)
+        else:
+            f11 = 0.0    
+        # total model -- Equation 1
+        mean = f1 + f78 + f5 + f4 + f6 + f10 + f11
 
         # ****** Aleatory uncertainty model ******
         # Intra-event term -- Equation 24
@@ -604,6 +614,7 @@ def get_IM(self, Mw, site_rup_dict, site_info, im_info):  # noqa: N802, N803, D1
                 vsInf,
                 site_info['z1pt0'] / 1000.0,
                 style,
+                site_rup_dict.get('crJB', None)
             )
             self.timeCalc += time.process_time_ns() - start
             meanList.append(mean)