Skip to content

Commit bfa46ff

Browse files
jinyan1214jinyan1214
committed
Fix a bug in hazard scenario downsampling
1 parent acb159e commit bfa46ff

2 files changed

Lines changed: 159 additions & 2 deletions

File tree

‎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

0 commit comments

Comments
 (0)