Skip to content

Commit fab04f3

Browse files
fmckennafmckenna
authored
Merge pull request #420 from jinyan1214/master
Fix bug reported in https://github.com/orgs/NHERI-SimCenter/discussions/363
2 parents 074838a + bfa46ff commit fab04f3

5 files changed

Lines changed: 219 additions & 9 deletions

File tree

‎modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py‎

Lines changed: 25 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,18 +54,21 @@
5454
'SA': [
5555
'Chiou & Youngs (2014)',
5656
'Abrahamson, Silva & Kamai (2014)',
57+
'Abrahamson, Silva & Kamai (2014) Aftershock',
5758
'Boore, Stewart, Seyhan & Atkinson (2014)',
5859
'Campbell & Bozorgnia (2014)',
5960
],
6061
'PGA': [
6162
'Chiou & Youngs (2014)',
6263
'Abrahamson, Silva & Kamai (2014)',
64+
'Abrahamson, Silva & Kamai (2014) Aftershock',
6365
'Boore, Stewart, Seyhan & Atkinson (2014)',
6466
'Campbell & Bozorgnia (2014)',
6567
],
6668
'PGV': [
6769
'Chiou & Youngs (2014)',
6870
'Abrahamson, Silva & Kamai (2014)',
71+
'Abrahamson, Silva & Kamai (2014) Aftershock',
6972
'Boore, Stewart, Seyhan & Atkinson (2014)',
7073
'Campbell & Bozorgnia (2014)',
7174
],
@@ -136,11 +139,16 @@ def __init__(
136139
gmpe_weights_dict=dict(), # noqa: B006, C408
137140
im_type=None,
138141
site_info=dict(), # noqa: B006, C408
142+
mainshock=None
139143
):
140144
# basic set-ups
141145
self.set_im_gmpe(im_dict, gmpe_dict, gmpe_weights_dict)
142146
self.set_im_type(im_type)
143147
self.set_sites(site_info)
148+
if mainshock is not None:
149+
self.mainshock = mainshock.copy() # single row gdf containing mainshock rupture surface
150+
else:
151+
self.mainshock = None
144152
# self.set_source(source_info)
145153

146154
def set_source(self, source_info): # noqa: D102
@@ -155,6 +163,7 @@ def set_source(self, source_info): # noqa: D102
155163
or 'Abrahamson, Silva & Kamai (2014)' in gmpe_list
156164
or 'Boore, Stewart, Seyhan & Atkinson (2014)' in gmpe_list
157165
or 'Campbell & Bozorgnia (2014)' in gmpe_list
166+
or 'Abrahamson, Silva & Kamai (2014) Aftershock' in gmpe_list
158167
):
159168
source_index = source_info.get('SourceIndex', None)
160169
rupture_index = source_info.get('RuptureIndex', None)
@@ -163,6 +172,14 @@ def set_source(self, source_info): # noqa: D102
163172
self.erf, source_index, rupture_index, self.site_info
164173
)
165174
# self.timeGetRuptureInfo += time.process_time_ns() - start
175+
if 'Abrahamson, Silva & Kamai (2014) Aftershock' in gmpe_list:
176+
source_index = source_info.get('SourceIndex', None)
177+
rupture_index = source_info.get('RuptureIndex', None)
178+
crJB = get_rupture_info_ASK2014_aftershock( # noqa: F405
179+
self.erf, source_index, rupture_index, self.mainshock
180+
)
181+
site_rup_dict['crJB'] = crJB
182+
166183
elif source_info['Type'] == 'PointSource':
167184
if (
168185
'Chiou & Youngs (2014)' in gmpe_list
@@ -456,7 +473,7 @@ def get_im_from_local( # noqa: C901, D102
456473
eq_magnitude, self.site_rup_dict, cur_site, im_info
457474
)
458475
# self.timeGetIM += time.process_time_ns() - start
459-
elif cur_gmpe == 'Abrahamson, Silva & Kamai (2014)':
476+
elif cur_gmpe == 'Abrahamson, Silva & Kamai (2014)' or cur_gmpe == 'Abrahamson, Silva & Kamai (2014) Aftershock':
460477
# start = time.process_time_ns()
461478
tmpResult = self.ASK.get_IM( # noqa: N806
462479
eq_magnitude, self.site_rup_dict, cur_site, im_info
@@ -793,13 +810,17 @@ def compute_im( # noqa: C901, D103
793810
sys.exit(
794811
'SA is used in hazard downsampling but not defined in the intensity measure tab'
795812
)
796-
elif ho_period in im_info['SA'].get('Periods'):
797-
pass
813+
elif im_info.get('Periods', None) is not None:
814+
if ho_period in im_info['Periods']:
815+
pass
816+
elif im_info.get('SA', None) is not None:
817+
if ho_period in im_info['SA'].get('Periods'):
818+
pass
798819
else:
799820
tmp_periods = im_info['SA']['Periods'] + [ho_period]
800821
tmp_periods.sort()
801822
im_info['SA']['Periods'] = tmp_periods
802-
elif ho_period in im_info['SA'].get('Periods'):
823+
elif ho_period in im_info.get('Periods'):
803824
pass
804825
else:
805826
tmp_periods = im_info['SA']['Periods'] + [ho_period]

‎modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -46,6 +46,8 @@
4646
import sys
4747
import psutil
4848
import GlobalVariable
49+
import shapely
50+
import geopandas as gpd
4951

5052
if 'stampede2' not in socket.gethostname():
5153
import GlobalVariable
@@ -395,6 +397,25 @@ def get_rupture_distance(erf, source_index, rupture_index, lat, lon): # noqa: D
395397

396398
return distToRupture
397399

400+
def get_rupture_info_ASK2014_aftershock(
401+
erf, source_index, rupture_indx, mainshock):
402+
rupSource = erf.getSource(source_index) # noqa: N806
403+
rupList = rupSource.getRuptureList() # noqa: N806
404+
rupSurface = rupList.get(rupture_indx).getRuptureSurface()
405+
rupSurface_perimeter = rupSurface.getPerimeter()
406+
coords = []
407+
for i in range(rupSurface_perimeter.size()):
408+
loc = rupSurface_perimeter.get(i)
409+
coords.append((loc.getLongitude(), loc.getLatitude()))
410+
if len(coords) == 1:
411+
rup_polygon = shapely.geometry.Point(coords[0])
412+
else:
413+
rup_polygon = shapely.geometry.Polygon(coords)
414+
rup_gdf = gpd.GeoDataFrame(index=[0], crs='EPSG:4326', geometry=[rup_polygon])
415+
rup_gdf = rup_gdf.to_crs(epsg=6417)
416+
centroid = rup_gdf.geometry.centroid.iloc[0]
417+
return mainshock.geometry.iloc[0].distance(centroid) /1000 # in meters
418+
398419

399420
def get_rupture_info_CY2014(erf, source_index, rupture_index, siteList): # noqa: N802, N803, D103
400421
rupSource = erf.getSource(source_index) # noqa: N806

‎modules/performRegionalEventSimulation/regionalGroundMotion/HazardOccurrence.py‎

Lines changed: 6 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -87,14 +87,18 @@ def configure_hazard_occurrence( # noqa: C901, D103
8787
# return periods
8888
if hc_input is None:
8989
return {}
90-
elif hc_input == 'Inferred_NSHMP': # noqa: RET505
90+
elif hc_input == 'NSHM V1' or hc_input == 'NSHM V2': # noqa: RET505
9191
period = hzo_config.get('Period', 0.0)
9292
if im_type == 'SA':
9393
cur_imt = im_type + f'{period:.1f}'.replace('.', 'P')
9494
else:
9595
cur_imt = im_type
9696
# fetching hazard curve from usgs
97-
cur_edition = hzo_config.get('Edition', 'E2014')
97+
version_str = hc_input.split(' ')[-1]
98+
if version_str == 'V1':
99+
cur_edition = 'E2008'
100+
else:
101+
cur_edition = 'E2014'
98102
hazard_curve_collector = []
99103
for site_id in range(len(site_config)):
100104
cur_site = site_config[site_id]

‎modules/performRegionalEventSimulation/regionalGroundMotion/LoadRupFile.py‎

Lines changed: 153 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,153 @@
1+
import numpy as np
2+
from tqdm import tqdm
3+
import geopandas as gpd
4+
import shapely
5+
import sys
6+
import ujson as json
7+
8+
def get_rups_to_run(scenario_info, num_scenarios): # noqa: C901, D103
9+
# If there is a filter
10+
if scenario_info['Generator'].get('method', None) == 'MonteCarlo':
11+
rup_filter = scenario_info['Generator'].get('RuptureFilter', None)
12+
if rup_filter is None or len(rup_filter) == 0:
13+
rups_to_run = list(range(num_scenarios))
14+
else:
15+
rups_requested = []
16+
for rups in rup_filter.split(','):
17+
if '-' in rups:
18+
asset_low, asset_high = rups.split('-')
19+
rups_requested += list(
20+
range(int(asset_low), int(asset_high) + 1)
21+
)
22+
else:
23+
rups_requested.append(int(rups))
24+
rups_requested = np.array(rups_requested)
25+
rups_requested = (
26+
rups_requested - 1
27+
) # The input index starts from 1, not 0
28+
rups_available = list(range(num_scenarios))
29+
rups_to_run = rups_requested[
30+
np.where(np.isin(rups_requested, rups_available))[0]
31+
]
32+
else:
33+
sys.exit(
34+
f'The scenario selection method {scenario_info["Generator"].get("method", None)} is not available'
35+
)
36+
return rups_to_run
37+
38+
def load_earthquake_rupFile(scenario_info, rupFilePath): # noqa: N802, N803, D103
39+
# Getting earthquake rupture forecast data
40+
source_type = scenario_info['EqRupture']['Type']
41+
try:
42+
with open(rupFilePath) as f: # noqa: PTH123
43+
user_scenarios = json.load(f)
44+
except: # noqa: E722
45+
sys.exit(f'CreateScenario: source file {rupFilePath} not found.')
46+
# number of features (i.e., ruptures)
47+
num_scenarios = len(user_scenarios.get('features', []))
48+
if num_scenarios < 1:
49+
sys.exit('CreateScenario: source file is empty.')
50+
rups_to_run = get_rups_to_run(scenario_info, num_scenarios)
51+
# get rupture and source ids
52+
scenario_data = {}
53+
if source_type == 'ERF':
54+
# source model
55+
source_model = scenario_info['EqRupture']['Model']
56+
for rup_tag in rups_to_run:
57+
cur_rup = user_scenarios.get('features')[rup_tag]
58+
cur_id_source = cur_rup.get('properties').get('Source', None)
59+
cur_id_rupture = cur_rup.get('properties').get('Rupture', None)
60+
scenario_data.update(
61+
{
62+
rup_tag: {
63+
'Type': source_type,
64+
'RuptureForecast': source_model,
65+
'Name': cur_rup.get('properties').get('Name', ''),
66+
'Magnitude': cur_rup.get('properties').get(
67+
'Magnitude', None
68+
),
69+
'MeanAnnualRate': cur_rup.get('properties').get(
70+
'MeanAnnualRate', None
71+
),
72+
'SourceIndex': cur_id_source,
73+
'RuptureIndex': cur_id_rupture,
74+
'SiteSourceDistance': cur_rup.get('properties').get(
75+
'Distance', None
76+
),
77+
'SiteRuptureDistance': cur_rup.get('properties').get(
78+
'DistanceRup', None
79+
),
80+
}
81+
}
82+
)
83+
elif source_type == 'PointSource':
84+
sourceID = 0 # noqa: N806
85+
rupID = 0 # noqa: N806
86+
for rup_tag in rups_to_run:
87+
try:
88+
cur_rup = user_scenarios.get('features')[rup_tag]
89+
magnitude = cur_rup.get('properties')['Magnitude']
90+
location = cur_rup.get('properties')['Location']
91+
average_rake = cur_rup.get('properties')['AverageRake']
92+
average_dip = cur_rup.get('properties')['AverageDip']
93+
scenario_data.update(
94+
{
95+
0: {
96+
'Type': source_type,
97+
'Magnitude': magnitude,
98+
'Location': location,
99+
'AverageRake': average_rake,
100+
'AverageDip': average_dip,
101+
'SourceIndex': sourceID,
102+
'RuptureIndex': rupID,
103+
}
104+
}
105+
)
106+
rupID = rupID + 1 # noqa: N806
107+
except: # noqa: PERF203, E722
108+
print('Please check point-source inputs.') # noqa: T201
109+
# return
110+
return scenario_data
111+
112+
def load_earthquake_rup_scenario(scenario_info, user_scenarios): # noqa: N802, N803, D103
113+
# Getting earthquake rupture forecast data
114+
source_type = scenario_info['EqRupture']['Type']
115+
# number of features (i.e., ruptures)
116+
num_scenarios = len(user_scenarios)
117+
if num_scenarios < 1:
118+
sys.exit('CreateScenario: source file is empty.')
119+
rups_to_run = get_rups_to_run(scenario_info, num_scenarios)
120+
# get rupture and source ids
121+
scenario_data = {}
122+
if source_type == 'ERF':
123+
# source model
124+
source_model = scenario_info['EqRupture']['Model']
125+
for rup_tag in rups_to_run:
126+
cur_rup = user_scenarios.iloc[rup_tag, :]
127+
cur_id_source = cur_rup['Source']
128+
cur_id_rupture = cur_rup['Rupture']
129+
print("DEBUG: cur_id_source, cur_id_rupture", cur_id_source, cur_id_rupture)
130+
scenario_data.update(
131+
{
132+
rup_tag: {
133+
'Type': source_type,
134+
'RuptureForecast': source_model,
135+
'Name': cur_rup.get('Name', ''),
136+
'Magnitude': cur_rup.get(
137+
'Magnitude', None
138+
),
139+
'MeanAnnualRate': cur_rup.get(
140+
'MeanAnnualRate', None
141+
),
142+
'SourceIndex': cur_id_source,
143+
'RuptureIndex': cur_id_rupture,
144+
'SiteSourceDistance': cur_rup.get(
145+
'Distance', None
146+
),
147+
'SiteRuptureDistance': cur_rup.get(
148+
'DistanceRup', None
149+
),
150+
}
151+
}
152+
)
153+
return scenario_data

‎modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py‎

Lines changed: 14 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -317,7 +317,7 @@ class abrahamson_silva_kamai_2014: # noqa: D101
317317
timeCalc = 0 # noqa: N815
318318
supportedImt = None # noqa: N815
319319

320-
def __init__(self):
320+
def __init__(self, aftershock=False):
321321
self.coeff = pd.read_csv(
322322
os.path.join(os.path.dirname(__file__), 'data', 'ASK14.csv') # noqa: PTH118, PTH120
323323
)
@@ -343,6 +343,7 @@ def __init__(self):
343343
self.H2 = 1.5
344344
self.H3 = -0.75
345345
self.PHI_AMP_SQ = 0.16
346+
self.aftershock = aftershock
346347

347348
def setIMT(self, imt): # noqa: N802, D102
348349
if imt not in self.supportedImt:
@@ -442,6 +443,7 @@ def calcValues( # noqa: C901, N802, D102
442443
vsInferred, # noqa: N803
443444
z1p0,
444445
style,
446+
crJB=None
445447
):
446448
if Mw > 5: # noqa: PLR2004
447449
c4mag = self.C4
@@ -537,8 +539,16 @@ def calcValues( # noqa: C901, N802, D102
537539
)
538540
else:
539541
f5 = (self.a10 + self.b * self.N) * np.log(vs30s / self.Vlin)
540-
# total model (no aftershock f11) -- Equation 1
541-
mean = f1 + f78 + f5 + f4 + f6 + f10
542+
543+
# Aftershock term -- Equation 20
544+
if self.aftershock:
545+
if crJB is None:
546+
raise ValueError('crJB must be provided for aftershock calculations')
547+
f11 = self.a14 * np.clip(1 - (crJB - 5)/10, 0, 1)
548+
else:
549+
f11 = 0.0
550+
# total model -- Equation 1
551+
mean = f1 + f78 + f5 + f4 + f6 + f10 + f11
542552

543553
# ****** Aleatory uncertainty model ******
544554
# 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
604614
vsInf,
605615
site_info['z1pt0'] / 1000.0,
606616
style,
617+
site_rup_dict.get('crJB', None)
607618
)
608619
self.timeCalc += time.process_time_ns() - start
609620
meanList.append(mean)

0 commit comments

Comments
 (0)