diff --git a/modules/performREC/CMakeLists.txt b/modules/performREC/CMakeLists.txt
index bf93c9eb8..960dbba6c 100644
--- a/modules/performREC/CMakeLists.txt
+++ b/modules/performREC/CMakeLists.txt
@@ -1,2 +1 @@
-simcenter_add_python_script(SCRIPT transportation.py)
add_subdirectory(pyrecodes)
\ No newline at end of file
diff --git a/modules/performREC/pyrecodes/run_pyrecodes.py b/modules/performREC/pyrecodes/run_pyrecodes.py
index 8a1fc6809..589b3e468 100644
--- a/modules/performREC/pyrecodes/run_pyrecodes.py
+++ b/modules/performREC/pyrecodes/run_pyrecodes.py
@@ -1,3 +1,35 @@
+import ctypes
+import sys
+
+# Programmatic fix for Apple Silicon macOS flat namespace errors
+if sys.platform == "darwin":
+ try:
+ # Force-load libomp globally so its symbols are visible to dlopen()
+ ctypes.CDLL('/opt/homebrew/opt/libomp/lib/libomp.dylib', mode=ctypes.RTLD_GLOBAL)
+ except OSError:
+ # Fallback for systems where Homebrew is not installed or configured differently
+ pass
+
+# Shim pandana.network -> pandarm so pyrecodes' bundled residual_demand_API
+# (which imports pandana at module load) can run in environments where
+# pandana is not installed. Only pdna.Network is used.
+import types
+import pandarm
+_pandana_network = types.ModuleType('pandana.network')
+_pandana_network.Network = pandarm.Network
+_pandana = types.ModuleType('pandana')
+_pandana.network = _pandana_network
+sys.modules['pandana'] = _pandana
+sys.modules['pandana.network'] = _pandana_network
+
+# Silence pandarm's CRS-default warning: residual_demand_API rebuilds the
+# Network every substep without passing a CRS, which floods the log.
+import warnings
+warnings.filterwarnings(
+ 'ignore',
+ message='No CRS was passed to geometry input',
+)
+
import pyrecodes.resilience_calculator
import pyrecodes.resilience_calculator.recodes_calculator
import pyrecodes.resilience_calculator.resilience_calculator
@@ -15,7 +47,7 @@
# sys.path.append('/Users/jinyanzhao/Desktop/SimCenterBuild/r2d_pyrecodes')
import pyrecodes
import pyrecodes.main
-from pyrecodes.geovisualizer.r2d_geovisualizer import R2D_GeoVisualizer
+from pyrecodes.geovisualizer.r2d_geovisualizer import R2DGeoVisualizer
from pyrecodes.plotter.concrete_plotter import ConcretePlotter
@@ -61,7 +93,7 @@ def select_realizations_to_run(damage_input, run_dir):
rlzs_requested.append(int(rlzs))
rlzs_requested = np.array(rlzs_requested)
- rlzs_in_available = np.in1d(rlzs_requested, rlzs_available) # noqa: NPY201
+ rlzs_in_available = np.isin(rlzs_requested, rlzs_available)
if rlzs_in_available.sum() != 0:
rlzs_to_run = rlzs_requested[np.where(rlzs_in_available)[0]]
else:
@@ -79,6 +111,51 @@ def select_realizations_to_run(damage_input, run_dir):
raise ValueError(msg)
return rlzs_to_run
+
+
+def find_display_recodes_calculator(system):
+ """
+ Locate the ReCoDeSCalculator whose output should drive the per-resource visualizations.
+
+ The "All" scope is the convention pyrecodes uses for region-wide
+ accounting, so that's what we match. If a user's system has no
+ such calculator, the caller is expected to handle the None result
+ rather than fall back silently to some other calculator.
+
+ Parameters
+ ----------
+ system : pyrecodes System
+ The system object returned by ``pyrecodes.main.run``.
+
+ Returns
+ -------
+ tuple[int | None, ReCoDeSCalculator | None]
+ ``(index, calculator)`` of the matched calculator, or
+ ``(None, None)`` if no match is configured.
+ """
+ target_cls = pyrecodes.resilience_calculator.recodes_calculator.ReCoDeSCalculator
+ for i, c in enumerate(system.resilience_calculators):
+ if isinstance(c, target_cls) and getattr(c, 'scope', None) == 'All':
+ return i, c
+ return None, None
+
+
+def find_recovery_time_calculator(system):
+ """
+ Return the first calculator that exposes a ``component_recovery_times`` list, or ``None``.
+
+ Both ``ComponentRecoveryTimeCalculator`` and
+ ``R2DComponentRecoveryTimeCalculator`` populate that attribute as a
+ list of length ``len(system.components)``, parallel to
+ ``system.components``. The wrapper just needs the recovery-time
+ value per component, so either subclass is acceptable.
+ """
+ for c in system.resilience_calculators:
+ if hasattr(c, 'component_recovery_times'):
+ return c
+ return None
+
+
def run_one_realization(main_file, rlz, rwhale_run_dir, system_config, save_pickle_file):
"""
Run a single realization of the pyrecodes simulation.
@@ -91,36 +168,44 @@ def run_one_realization(main_file, rlz, rwhale_run_dir, system_config, save_pick
system_config (dict): System configuration dictionary.
"""
# Run the pyrecodes
- system = pyrecodes.main.run(main_file)
+ main_file_path = Path(main_file)
+ system = pyrecodes.main.run(str(main_file_path.parent), main_file_path.name)
system.calculate_resilience(print_output=False)
# Add the component recovery time to the Results_rlz.json file
with Path(rwhale_run_dir / f'Results_{rlz}.json').open() as f:
results_rlz = json.load(f)
- all_recovery_time = system.resilience_calculators[1].component_recovery_times
- for ind, comp in enumerate(system.components):
- if getattr(comp, 'r2d_comp', False) is True:
- asset_type = comp.asset_type
- asset_subtype = comp.asset_subtype
- asset_id = comp.general_information['AIM_id']
- results_rlz[asset_type][asset_subtype][asset_id]['Recovery'] = {
- 'Time': next(iter(all_recovery_time[ind].values()))
- }
+ recovery_calc = find_recovery_time_calculator(system)
+ if recovery_calc is None:
+ print( # noqa: T201
+ f'Warning: no component-recovery-time calculator found in system; '
+ f'skipping recovery-time merge for realization {rlz}.'
+ )
+ else:
+ all_recovery_time = recovery_calc.component_recovery_times
+ for ind, comp in enumerate(system.components):
+ if getattr(comp, 'r2d_comp', False) is True:
+ asset_type = comp.asset_type
+ asset_subtype = comp.asset_subtype
+ asset_id = comp.general_information['AIM_id']
+ results_rlz[asset_type][asset_subtype][asset_id]['Recovery'] = {
+ 'Time': next(iter(all_recovery_time[ind].values()))
+ }
# Write the results to a file in the current realization workdir
with (Path(f'Results_{rlz}.json')).open('w') as f:
json.dump(results_rlz, f)
# Create a gif of the recovery process
- geo_visualizer = R2D_GeoVisualizer(system.components)
+ geo_visualizer = R2DGeoVisualizer(system.components)
# time_step_list = list(range(0, system.time_step, 1))
time_step_list = list(np.linspace(0, system.time_step, min(system.time_step, 20)).astype(int))
for time_step in time_step_list:
fig, ax = plt.subplots(figsize=(10, 10))
- geo_visualizer.create_current_state_figure(time_step, ax=ax)
+ geo_visualizer.create_current_recovery_state_figure(time_step, ax=ax)
legend = ax.get_legend()
legend.set(loc='center left', bbox_to_anchor=(1, 0.5))
plt.savefig(f'status_at_time_step_{time_step}.png', dpi=300, bbox_inches='tight', transparent=False, pad_inches=0)
- geo_visualizer.create_recovery_gif(time_step_list, file_name = \
+ geo_visualizer.create_recovery_gif(time_step_list, load_file_name = \
'status_at_time_step_TIME_STEP.png', fps=2)
# create a plot of the supply and demand recovery of the first resource
@@ -137,33 +222,47 @@ def run_one_realization(main_file, rlz, rwhale_run_dir, system_config, save_pick
units_to_plot.append(system_config['Resources']['TransportationService'].get('Unit', 'unit_TransportationService'))
print(f'Resources to plot {resources_to_plot}')
+ # Write one supply/demand/consumption JSON per (ReCoDeS calculator,
+ # resource) pair, with the calculator's positional index as the
+ # filename suffix.
for calculator_i, calculator in enumerate(system.resilience_calculators):
if isinstance(calculator, pyrecodes.resilience_calculator.recodes_calculator.ReCoDeSCalculator):
- resources_to_plot_i = []
- for resource in resources_to_plot:
- if resource in calculator.resource_names:
- resources_to_plot_i.append(resource)
+ resources_to_plot_i = [
+ resource for resource in resources_to_plot
+ if resource in calculator.resource_names
+ ]
if len(resources_to_plot_i) > 0:
plotter_object.save_supply_demand_consumption(system, resources_to_plot_i, calculator_i)
- # plotter_object.save_component_recovery_progress(system.components[:20])
-
- for resource, unit in zip(resources_to_plot, units_to_plot):
- y_axis_label = f'{resource} {unit} | {system.resilience_calculators[0].scope}'
- x_axis_label = 'Time step [day]'
- axis_object = plotter_object.setup_lor_plot_fig(x_axis_label, y_axis_label)
- time_range = system.time_step+1
- time_steps_before_event = 10
- plotter_object.plot_single_resource(list(range(-time_steps_before_event, time_range)), system.resilience_calculators[0].system_supply[resource][:time_range],
- system.resilience_calculators[0].system_demand[resource][:time_range],
- system.resilience_calculators[0].system_consumption[resource][:time_range], axis_object, warmup=time_steps_before_event,
- show = False
- )
- plotter_object.save_current_figure(savename = f'{resource}_supply_demand_consumption.png')
-
- plotter_object.save_supply_demand_consumption(system, [resource])
-
- plotter_object.save_component_recovery_progress(system.components)
+ # Render the per-resource PNG plots from the "All" scope calculator.
+ display_calc_idx, display_calc = find_display_recodes_calculator(system)
+ if display_calc is None:
+ print( # noqa: T201
+ 'Warning: no ReCoDeSCalculator with scope="All" found in system; '
+ 'skipping per-resource supply/demand PNG plots.'
+ )
+ else:
+ for resource, unit in zip(resources_to_plot, units_to_plot):
+ if resource not in display_calc.resource_names:
+ continue
+ y_axis_label = f'{resource} {unit} | {display_calc.scope}'
+ x_axis_label = 'Time step [day]'
+ axis_object = plotter_object.setup_lor_plot_fig(x_axis_label, y_axis_label)
+ time_range = system.time_step + 1
+ time_steps_before_event = 10
+ plotter_object.plot_single_resource(
+ list(range(-time_steps_before_event, time_range)),
+ display_calc.system_supply[resource][:time_range],
+ display_calc.system_demand[resource][:time_range],
+ display_calc.system_consumption[resource][:time_range],
+ axis_object,
+ warmup=time_steps_before_event,
+ show=False,
+ )
+ plotter_object.save_current_figure(savename=f'{resource}_supply_demand_consumption.png')
+
+ plotter_object.save_component_recovery_progress(system.components)
+
if save_pickle_file:
# Save the system components as a pickle file
with open(f'system_{rlz}.pickle', 'wb') as f:
@@ -258,7 +357,7 @@ def modify_system_config_residual_demand_distribution(system_config, input_data_
def modify_main_file(main_file_dict, component_library, run_dir):
"""
- Modify the main configuration file with updated paths for the component library and system configuration.
+ Copy the component library into run_dir and write main.json there with basename references to the staged inputs.
Parameters
----------
@@ -266,17 +365,18 @@ def modify_main_file(main_file_dict, component_library, run_dir):
The main configuration dictionary.
component_library : str
Path to the component library file.
- rwhale_run_dir : str
- Directory where the results are stored.
+ run_dir : Path
+ Per-realization run directory.
Returns
-------
- bool
- True if the modification is successful.
+ str
+ Path to the written main.json file.
"""
- main_file_dict['System']['SystemConfigurationFile'] = str(
- run_dir / 'SystemConfiguration.json')
- main_file_dict['ComponentLibrary']['ComponentLibraryFile'] = component_library
+ component_library_basename = Path(component_library).name
+ shutil.copy2(component_library, run_dir / component_library_basename)
+ main_file_dict['ComponentLibrary']['ComponentLibraryFile'] = component_library_basename
+ main_file_dict['System']['SystemConfigurationFile'] = 'SystemConfiguration.json'
main_file_path = str(run_dir / 'main.json')
with Path(main_file_path).open('w') as f:
json.dump(main_file_dict, f)
@@ -520,6 +620,11 @@ def run_pyrecodes( # noqa: C901
print(f'Rrealization {rlz} completed') # noqa: T201
# Write the aggregated results to the Results_det.json file
aggregate_results_to_det(results_agg, Path(run_dir / 'Results_det.json'))
+ print( # noqa: T201
+ f'Pyrecodes recovery simulation completed. '
+ f'Aggregated results written to {run_dir / "Results_det.json"}. '
+ f'Per-realization outputs are in {run_dir / "RecoverySimulation"}.'
+ )
else:
msg = 'Damage input type not recognized'
raise ValueError(msg)
diff --git a/modules/performREC/transportation.py b/modules/performREC/transportation.py
deleted file mode 100644
index dc658ec56..000000000
--- a/modules/performREC/transportation.py
+++ /dev/null
@@ -1,2049 +0,0 @@
-"""Methods for performance simulations of transportation networks.""" # noqa: INP001
-
-# -*- coding: utf-8 -*-
-#
-# Copyright (c) 2024 The Regents of the University of California
-#
-# This file is part of SimCenter Backend Applications.
-#
-# Redistribution and use in source and binary forms, with or without
-# modification, are permitted provided that the following conditions are met:
-#
-# 1. Redistributions of source code must retain the above copyright notice,
-# this list of conditions and the following disclaimer.
-#
-# 2. Redistributions in binary form must reproduce the above copyright notice,
-# this list of conditions and the following disclaimer in the documentation
-# and/or other materials provided with the distribution.
-#
-# 3. Neither the name of the copyright holder nor the names of its contributors
-# may be used to endorse or promote products derived from this software without
-# specific prior written permission.
-#
-# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
-# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
-# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
-# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
-# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
-# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
-# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
-# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
-# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
-# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
-# POSSIBILITY OF SUCH DAMAGE.
-#
-# You should have received a copy of the BSD 3-Clause License along with
-# BRAILS. If not, see .
-#
-# Contributors:
-# Barbaros Cetiner
-# Tianyu Han
-#
-# Last updated:
-# 08-14-2024
-
-from __future__ import annotations # noqa: I001
-
-import gc
-import os
-import json
-import logging
-import time
-import copy
-from collections import defaultdict
-from abc import ABC, abstractmethod
-from pathlib import Path
-from typing import List, Set, Tuple, Dict, Literal, Union, Optional
-from shapely.geometry import LineString, Polygon, Point, mapping
-from shapely.ops import split
-from sklearn.metrics.pairwise import cosine_similarity
-
-import geopandas as gpd
-import numpy as np
-import pandana.network as pdna
-import pandas as pd
-from sklearn.feature_extraction.text import TfidfVectorizer
-from scipy.spatial.distance import cdist
-from shapely.wkt import loads
-from brails.utils.geoTools import haversine_dist
-from brails.workflow.TransportationElementHandler import (
- ROADLANES_MAP,
- ROADSPEED_MAP,
- calculate_road_capacity,
-)
-
-
-class TransportationPerformance(ABC): # noqa: B024
- """
- An abstract base class for simulating transportation networks.
-
- This class provides an interface for implementing methods that process
- transportation data (such as system state and origin-destination files) and
- compute network performance metrics.
-
- Subclasses must define how to process these data inputs and update system
- performance in concrete terms.
-
- Attributes__
- assets (list): A list of asset types (e.g., 'Bridge', 'Roadway',
- 'Tunnel') to be analyzed in the transportation network.
- csv_filenames (list): A list of filenames (e.g., 'edges.csv', '
- nodes.csv', 'closed_edges.csv') that are required
- for network simulations.
- capacity_map (dict): A mapping that relates damage states to capacity
- ratios. For example, damage states of 0-2 may
- represent full capacity (1), while higher damage
- states reduce capacity (e.g., 0.5 or 0).
- no_identifier (dict): A mapping of asset types to their unique
- identifiers (e.g., 'StructureNumber' for Bridges,
- 'TigerOID' for Roadways). These are used to track
- and manage assets within the system.
-
- Methods__
- system_state(detfile: str, csv_file_dir: str) -> None:
- Abstract method to process a given det (damage state) file and
- update the system state with capacity ratios. Also checks for
- missing or necessary CSV files for network simulation.
-
- update_od_file(old_nodes: str, old_det: str, new_nodes: str, new_det:
- str, od: str, origin_ids: list) -> pd.DataFrame:
- Abstract method to update the origin-destination (OD) file based
- on population changes between time steps. The method tracks
- population growth in nodes and generates new trips in the OD file
- accordingly.
-
- system_performance(detfile: str, csv_file_dir: str) -> None:
- Abstract method to compute or update the performance of the
- transportation system based on current state and available data.
-
- Notes__
- This is an abstract class. To use it, create a subclass that implements
- the abstract methods for specific behavior related to transportation
- network performance analysis.
- """
-
- def __init__(
- self,
- assets: Union[list[str], None] = None, # noqa: UP007
- capacity_map: Union[dict[int, float], None] = None, # noqa: UP007
- csv_files: Union[dict[str, str], None] = None, # noqa: UP007
- no_identifier: Union[dict[str, str], None] = None, # noqa: UP007
- network_inventory: Union[dict[str, str], None] = None, # noqa: UP007
- ):
- """
- Initialize the TransportationPerformance class with essential data.
-
- Args__
- assets (list): A list of asset types such as 'Bridge', 'Roadway',
- 'Tunnel'.
- capacity_map (dict): A mapping of damage states to capacity ratios.
- csv_files (dict): A dictionary of CSV filenames for network data,
- including 'network_edges', 'network_nodes',
- 'edge_closures', and 'od_pairs'.
- no_identifier (dict): A mapping of asset types to their unique
- identifiers.
- network_inventory (dict): A dictionary of "edges" and "nodes" file.
- """
- if assets is None:
- assets = ['Bridge', 'Roadway', 'Tunnel']
- if capacity_map is None:
- # capacity_map = {0: 1, 1: 1, 2: 1, 3: 0.5, 4: 0}
- capacity_map = {0: 1, 1: 1, 2: 1, 3: 0, 4: 0}
- if csv_files is None:
- csv_files = {
- 'network_edges': 'edges.csv',
- 'network_nodes': 'nodes.csv',
- 'edge_closures': 'closed_edges.csv',
- 'od_pairs': 'od.csv',
- }
- if no_identifier is None:
- no_identifier = {
- 'Bridge': 'StructureNumber',
- 'Roadway': 'TigerOID',
- 'Tunnel': 'TunnelNumber',
- }
-
- self.assets = assets
- self.csv_files = csv_files
- self.capacity_map = capacity_map
- self.no_identifier = no_identifier
- self.attribute_maps = {
- 'lanes': ROADLANES_MAP,
- 'speed': ROADSPEED_MAP,
- 'capcity': calculate_road_capacity,
- }
- self.network_inventory = network_inventory
-
- # @abstractmethod
- def system_state(
- self, initial_state: str, csv_file_dir: str
- ) -> dict: # updated_state
- """
- Process given det and damage results file to get updated system state.
-
- This function reads a JSON file containing undamaged network attributes
- and JSON file containing damage states and updates the state of the
- network (i.e., determines edges experiencing capacity reductions)
- using capacity ratios defined for each damage state. It also checks
- for the existence of required CSV files, created the if missing, and
- generates a file listing closed edges.
-
- Args__
- detfile (str): Path to the JSON file containing the asset data.
- csv_file_dir (str): Directory containing the CSV files needed for
- running network simulations.
-
- Returns__
- The function does not return any value. It creates updated det file
- and CSV file necessary to run network simulations
-
- Raises__
- FileNotFoundError: If the `detfile` or any required CSV files are
- not found.
- json.JSONDecodeError: If the `detfile` cannot be parsed as JSON.
- KeyError: If expected keys are missing in the JSON data.
- ValueError: If there are issues with data conversions, e.g., JSON
- to integer.
-
- Examples__
- >>> system_state('damage_states.json', '/path/to/csv/files')
- Missing files: nodes.csv
- All required files are present.
- # This will process 'damage_states.json', update it, and use CSV
- files in '/path/to/csv/files'
-
- >>> system_state('damage_states.json', '/path/to/nonexistent/dir')
- Missing files: edges.csv, nodes.csv
- # This will attempt to process 'damage_states.json' and handle
- missing files in a non-existent directory
- """
-
- def files_exist(directory, filenames):
- # Convert directory path to a Path object
- dir_path = Path(directory)
-
- # Get a set of files in the directory
- files_in_directory = {f.name for f in dir_path.iterdir() if f.is_file()}
-
- # Check if each file exists in the directory
- missing_files = [
- filename
- for filename in filenames
- if filename not in files_in_directory
- ]
-
- if missing_files:
- print(f"Missing files: {', '.join(missing_files)}") # noqa: T201
- out = False
- else:
- print('All required files are present.') # noqa: T201
- out = True
-
- return out
-
- # Create link closures for network simulations:
- fexists = files_exist(csv_file_dir, [self.csv_files['network_edges']])
-
- if fexists:
- pass
- # graph_edge_file = os.path.join(csv_file_dir, self.csv_files['network_edges'])
- # graph_edge_df = pd.read_csv(graph_edge_file)
- else:
- self.get_graph_network(csv_file_dir)
-
- # Read damage states for det file and determine element capacity ratios
- # 1 denotes fully open and 0 denotes fully closed:
- capacity_dict = {}
- closed_edges = []
- with Path.open(Path(initial_state), encoding='utf-8') as file:
- temp = json.load(file)
- data = temp['TransportationNetwork']
- for asset_type in self.assets:
- datadict = data[asset_type]
- for aim_id in datadict:
- item_id = datadict[aim_id]['GeneralInformation'][
- self.no_identifier[asset_type]
- ]
- damage_state = int(
- datadict[aim_id]['R2Dres'][
- 'R2Dres_MostLikelyCriticalDamageState'
- ]
- )
- capacity_ratio = self.capacity_map[damage_state]
- capacity_dict[item_id] = capacity_ratio
- datadict[aim_id]['GeneralInformation']['Open'] = capacity_ratio
- if capacity_ratio == 0:
- if asset_type == 'Roadway':
- closed_edges.append(int(aim_id))
- elif asset_type == 'Bridge' or asset_type == 'Tunnel': # noqa: PLR1714
- carried_edges = datadict[aim_id]['GeneralInformation'][
- 'RoadID'
- ]
- carried_edges = [int(x) for x in carried_edges]
- closed_edges += carried_edges
-
- # Update det file with closure information:
- temp = os.path.basename(initial_state).split('.') # noqa: PTH119
- detfile_updated = temp[0] + '_updated.' + temp[1]
- detfile_updated = os.path.join(csv_file_dir, detfile_updated) # noqa: PTH118
-
- with Path.open(Path(detfile_updated), 'w', encoding='utf-8') as file:
- json.dump(data, file, indent=2)
-
- # Write closed edges:
- edge_closure_file = os.path.join( # noqa: PTH118
- csv_file_dir, self.csv_files['edge_closures']
- )
- with open(edge_closure_file, 'w', encoding='utf-8') as file: # noqa: PTH123
- # Write each item on a new line
- file.write('uniqueid\n')
- for item in closed_edges:
- file.write(str(item) + '\n')
-
- return
-
-
- def get_graph_network(self, csv_file_dir) -> None: # noqa: D102
- # Get edges and nodes from the network inventory
- edges_file = self.network_inventory['edges']
- nodes_file = self.network_inventory['nodes']
- edges_gpd = gpd.read_file(edges_file)
- nodes_gpd = gpd.read_file(nodes_file)
- # Format edges and nodes for Pandana
- ## Edges
- edges_gpd['id'] = edges_gpd['id'].astype(int)
- edges_gpd = edges_gpd.rename(
- columns={
- 'id': 'uniqueid',
- 'StartNode': 'start_nid',
- 'EndNode': 'end_nid',
- 'RoadType': 'type',
- 'NumOfLanes': 'lanes',
- 'MaxMPH': 'maxspeed',
- }
- )
- edges_gpd['capacity'] = edges_gpd['lanes'].map(calculate_road_capacity)
- edges_gpd = edges_gpd.to_crs('EPSG:6500')
- edges_gpd['length'] = edges_gpd['geometry'].apply(lambda x: x.length)
- edges_gpd = edges_gpd.to_crs('EPSG:4326')
- edges_fils = os.path.join(csv_file_dir, self.csv_files['network_edges']) # noqa: PTH118
- edges_gpd.to_csv(edges_fils, index=False)
- ## Nodes
- nodes_gpd = nodes_gpd.rename(columns={'nodeID': 'node_id'})
- nodes_gpd['lon'] = nodes_gpd['geometry'].apply(lambda x: x.x)
- nodes_gpd['lat'] = nodes_gpd['geometry'].apply(lambda x: x.y)
- nodes_file = os.path.join(csv_file_dir, self.csv_files['network_nodes']) # noqa: PTH118
- nodes_gpd.to_csv(nodes_file, index=False)
- return # noqa: PLR1711
-
- # @abstractmethod
- def system_performance( # noqa: C901, D102, PLR0915
- self, state # noqa: ARG002
- ) -> None: # Move the CSV creation here
- def substep_assignment(
- nodes_df=None,
- weighted_edges_df=None,
- od_ss=None,
- quarter_demand=None,
- assigned_demand=None,
- quarter_counts=4,
- trip_info=None,
- agent_time_limit=0,
- sample_interval=1,
- agents_path=None,
- hour=None,
- quarter=None,
- ss_id=None,
- alpha_f=0.3,
- beta_f=3,
- ):
- open_edges_df = weighted_edges_df.loc[weighted_edges_df['fft'] < 36000] # noqa: PLR2004
-
- net = pdna.Network(
- nodes_df['x'],
- nodes_df['y'],
- open_edges_df['start_nid'],
- open_edges_df['end_nid'],
- open_edges_df[['weight']],
- twoway=False,
- )
-
- print('network') # noqa: T201
- net.set(pd.Series(net.node_ids))
- print('net') # noqa: T201
-
- nodes_origin = od_ss['origin_nid'].to_numpy()
- nodes_destin = od_ss['destin_nid'].to_numpy()
- nodes_current = od_ss['current_nid'].to_numpy()
- agent_ids = od_ss['agent_id'].to_numpy()
- agent_current_links = od_ss['current_link'].to_numpy()
- agent_current_link_times = od_ss['current_link_time'].to_numpy()
- paths = net.shortest_paths(nodes_current, nodes_destin)
-
- # check agent time limit
- path_lengths = net.shortest_path_lengths(nodes_current, nodes_destin)
- remove_agent_list = []
- if agent_time_limit is None:
- pass
- else:
- for agent_idx, agent_id in enumerate(agent_ids):
- planned_trip_length = path_lengths[agent_idx]
- # agent_time_limit[agent_id]
- trip_length_limit = agent_time_limit
- if planned_trip_length > trip_length_limit + 0:
- remove_agent_list.append(agent_id)
-
- edge_travel_time_dict = weighted_edges_df['t_avg'].T.to_dict()
- edge_current_vehicles = weighted_edges_df['veh_current'].T.to_dict()
- edge_quarter_vol = weighted_edges_df['vol_true'].T.to_dict()
- # edge_length_dict = weighted_edges_df['length'].T.to_dict()
- od_residual_ss_list = []
- # all_paths = []
- path_i = 0
- for path in paths:
- trip_origin = nodes_origin[path_i]
- trip_destin = nodes_destin[path_i]
- agent_id = agent_ids[path_i]
- # remove some agent (path too long)
- if agent_id in remove_agent_list:
- path_i += 1
- # no need to update trip info
- continue
- remaining_time = (
- 3600 / quarter_counts + agent_current_link_times[path_i]
- )
- used_time = 0
- for edge_s, edge_e in zip(path, path[1:]):
- edge_str = f'{edge_s}-{edge_e}'
- edge_travel_time = edge_travel_time_dict[edge_str]
-
- if (remaining_time > edge_travel_time) and (
- edge_travel_time < 36000 # noqa: PLR2004
- ):
- # all_paths.append(edge_str)
- # p_dist += edge_travel_time
- remaining_time -= edge_travel_time
- used_time += edge_travel_time
- edge_quarter_vol[edge_str] += 1 * sample_interval
- trip_stop = edge_e
-
- if edge_str == agent_current_links[path_i]:
- edge_current_vehicles[edge_str] -= 1 * sample_interval
- else:
- if edge_str != agent_current_links[path_i]:
- edge_current_vehicles[edge_str] += 1 * sample_interval
- new_current_link = edge_str
- new_current_link_time = remaining_time
- trip_stop = edge_s
- od_residual_ss_list.append(
- [
- agent_id,
- trip_origin,
- trip_destin,
- trip_stop,
- new_current_link,
- new_current_link_time,
- ]
- )
- break
- trip_info[(agent_id, trip_origin, trip_destin)][0] += (
- 3600 / quarter_counts
- )
- trip_info[(agent_id, trip_origin, trip_destin)][1] += used_time
- trip_info[(agent_id, trip_origin, trip_destin)][2] = trip_stop
- trip_info[(agent_id, trip_origin, trip_destin)][3] = hour
- trip_info[(agent_id, trip_origin, trip_destin)][4] = quarter
- trip_info[(agent_id, trip_origin, trip_destin)][5] = ss_id
- path_i += 1
-
- new_edges_df = weighted_edges_df[
- [
- 'uniqueid',
- 'u',
- 'v',
- 'start_nid',
- 'end_nid',
- 'fft',
- 'capacity',
- 'normal_fft',
- 'normal_capacity',
- 'length',
- 'vol_true',
- 'vol_tot',
- 'veh_current',
- 'geometry',
- ]
- ].copy()
- # new_edges_df = new_edges_df.join(edge_volume, how='left')
- # new_edges_df['vol_ss'] = new_edges_df['vol_ss'].fillna(0)
- # new_edges_df['vol_true'] += new_edges_df['vol_ss']
- new_edges_df['vol_true'] = new_edges_df.index.map(edge_quarter_vol)
- new_edges_df['veh_current'] = new_edges_df.index.map(
- edge_current_vehicles
- )
- # new_edges_df['vol_tot'] += new_edges_df['vol_ss']
- new_edges_df['flow'] = (
- new_edges_df['vol_true'] * quarter_demand / assigned_demand
- ) * quarter_counts
- new_edges_df['t_avg'] = new_edges_df['fft'] * (
- 1
- + alpha_f
- * (new_edges_df['flow'] / new_edges_df['capacity']) ** beta_f
- )
- new_edges_df['t_avg'] = np.where(
- new_edges_df['t_avg'] > 36000, 36000, new_edges_df['t_avg'] # noqa: PLR2004
- )
- new_edges_df['t_avg'] = new_edges_df['t_avg'].round(2)
-
- return new_edges_df, od_residual_ss_list, trip_info, agents_path
-
- def write_edge_vol(
- edges_df=None,
- simulation_outputs=None,
- quarter=None,
- hour=None,
- scen_nm=None,
- ):
- if 'flow' in edges_df.columns:
- if edges_df.shape[0] < 10: # noqa: PLR2004
- edges_df[
- [
- 'uniqueid',
- 'start_nid',
- 'end_nid',
- 'capacity',
- 'veh_current',
- 'vol_true',
- 'vol_tot',
- 'flow',
- 't_avg',
- 'geometry',
- ]
- ].to_csv(
- f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_'
- f'qt{quarter}_{scen_nm}.csv',
- index=False,
- )
-
- else:
- edges_df.loc[
- edges_df['vol_true'] > 0,
- [
- 'uniqueid',
- 'start_nid',
- 'end_nid',
- 'capacity',
- 'veh_current',
- 'vol_true',
- 'vol_tot',
- 'flow',
- 't_avg',
- 'geometry',
- ],
- ].to_csv(
- f'{simulation_outputs}/edge_vol/edge_vol_hr{hour}_'
- f'qt{quarter}_{scen_nm}.csv',
- index=False,
- )
-
- def write_final_vol(
- edges_df=None,
- simulation_outputs=None,
- quarter=None,
- hour=None,
- scen_nm=None,
- ):
- edges_df.loc[
- edges_df['vol_tot'] > 0,
- ['uniqueid', 'start_nid', 'end_nid', 'vol_tot', 'geometry'],
- ].to_csv(
- f'{simulation_outputs}/edge_vol/final_edge_vol_hr{hour}_qt'
- f'{quarter}_{scen_nm}.csv',
- index=False,
- )
-
- def assignment( # noqa: C901, PLR0913
- quarter_counts=6,
- substep_counts=15,
- substep_size=30000,
- edges_df=None,
- nodes_df=None,
- od_all=None,
- simulation_outputs=None,
- scen_nm=None,
- hour_list=None,
- quarter_list=None,
- cost_factor=None, # noqa: ARG001
- closure_hours=None,
- closed_links=None,
- agent_time_limit=None,
- sample_interval=1,
- agents_path=None,
- alpha_f=0.3,
- beta_f=4,
- ):
- if closure_hours is None:
- closure_hours = []
-
- # Check if all od pairs has path. If not,
- orig = od_all['origin_nid'].to_numpy()
- dest = od_all['destin_nid'].to_numpy()
- open_edges_df = edges_df[
- ~edges_df['uniqueid'].isin(closed_links['uniqueid'].to_numpy())
- ]
- net = pdna.Network(
- nodes_df['x'],
- nodes_df['y'],
- open_edges_df['start_nid'],
- open_edges_df['end_nid'],
- open_edges_df[['fft']],
- twoway=False,
- )
- paths = net.shortest_paths(orig, dest)
- no_path_ind = [i for i in range(len(paths)) if len(paths[i]) == 0]
- od_no_path = od_all.iloc[no_path_ind].copy()
- od_all = od_all.drop(no_path_ind)
-
- od_all['current_nid'] = od_all['origin_nid']
- trip_info = {
- (od.agent_id, od.origin_nid, od.destin_nid): [
- 0,
- 0,
- od.origin_nid,
- 0,
- od.hour,
- od.quarter,
- 0,
- 0,
- ]
- for od in od_all.itertuples()
- }
-
- # Quarters and substeps
- # probability of being in each division of hour
- if quarter_list is None:
- quarter_counts = 4
- else:
- quarter_counts = len(quarter_list)
- quarter_ps = [1 / quarter_counts for i in range(quarter_counts)]
- quarter_ids = list(range(quarter_counts))
-
- # initial setup
- od_residual_list = []
- # accumulator
- edges_df['vol_tot'] = 0
- edges_df['veh_current'] = 0
-
- # Loop through days and hours
- for _day in ['na']:
- for hour in hour_list:
- gc.collect()
- if hour in closure_hours:
- for row in closed_links.itertuples():
- edges_df.loc[
- (edges_df['uniqueid'] == row.uniqueid), 'capacity'
- ] = 1
- edges_df.loc[
- (edges_df['uniqueid'] == row.uniqueid), 'fft'
- ] = 36000
- else:
- edges_df['capacity'] = edges_df['normal_capacity']
- edges_df['fft'] = edges_df['normal_fft']
-
- # Read OD
- od_hour = od_all[od_all['hour'] == hour].copy()
- if od_hour.shape[0] == 0:
- od_hour = pd.DataFrame([], columns=od_all.columns)
- od_hour['current_link'] = None
- od_hour['current_link_time'] = 0
-
- # Divide into quarters
- if 'quarter' in od_all.columns:
- pass
- else:
- od_quarter_msk = np.random.choice(
- quarter_ids, size=od_hour.shape[0], p=quarter_ps
- )
- od_hour['quarter'] = od_quarter_msk
-
- if quarter_list is None:
- quarter_list = quarter_ids
- for quarter in quarter_list:
- # New OD in assignment period
- od_quarter = od_hour.loc[
- od_hour['quarter'] == quarter,
- [
- 'agent_id',
- 'origin_nid',
- 'destin_nid',
- 'current_nid',
- 'current_link',
- 'current_link_time',
- ],
- ]
- # Add resudal OD
- od_residual = pd.DataFrame(
- od_residual_list,
- columns=[
- 'agent_id',
- 'origin_nid',
- 'destin_nid',
- 'current_nid',
- 'current_link',
- 'current_link_time',
- ],
- )
- od_residual['quarter'] = quarter
- # Total OD in each assignment period is the combined
- # of new and residual OD:
- od_quarter = pd.concat(
- [od_quarter, od_residual], sort=False, ignore_index=True
- )
- # Residual OD is no longer residual after it has been
- # merged to the quarterly OD:
- od_residual_list = []
- od_quarter = od_quarter[
- od_quarter['current_nid'] != od_quarter['destin_nid']
- ]
-
- # total demand for this quarter, including total and
- # residual demand:
- quarter_demand = od_quarter.shape[0]
- # how many among the OD pairs to be assigned in this
- # quarter are actually residual from previous quarters
- residual_demand = od_residual.shape[0]
- assigned_demand = 0
-
- substep_counts = max(1, (quarter_demand // substep_size) + 1)
- logging.info(
- f'HR {hour} QT {quarter} has '
- f'{quarter_demand}/{residual_demand} od'
- f's/residuals {substep_counts} substeps'
- )
- substep_ps = [
- 1 / substep_counts for i in range(substep_counts)
- ]
- substep_ids = list(range(substep_counts))
- od_substep_msk = np.random.choice(
- substep_ids, size=quarter_demand, p=substep_ps
- )
- od_quarter['ss_id'] = od_substep_msk
-
- # reset volume at each quarter
- edges_df['vol_true'] = 0
-
- for ss_id in substep_ids:
- gc.collect()
-
- time_ss_0 = time.time()
- logging.info(
- f'Hour: {hour}, Quarter: {quarter}, '
- 'SS ID: {ss_id}'
- )
- od_ss = od_quarter[od_quarter['ss_id'] == ss_id]
- assigned_demand += od_ss.shape[0]
- if assigned_demand == 0:
- continue
- # calculate weight
- weighted_edges_df = edges_df.copy()
- # weight by travel distance
- # weighted_edges_df['weight'] = edges_df['length']
- # weight by travel time
- # weighted_edges_df['weight'] = edges_df['t_avg']
- # weight by travel time
- # weighted_edges_df['weight'] = (edges_df['t_avg']
- # - edges_df['fft']) * 0.5 + edges_df['length']*0.1
- # + cost_factor*edges_df['length']*0.1*(
- # edges_df['is_highway'])
- # 10 yen per 100 m --> 0.1 yen per m
- weighted_edges_df['weight'] = edges_df['t_avg']
- # weighted_edges_df['weight'] = np.where(
- # weighted_edges_df['weight']<0.1, 0.1,
- # weighted_edges_df['weight'])
-
- # traffic assignment with truncated path
- (
- edges_df,
- od_residual_ss_list,
- trip_info,
- agents_path,
- ) = substep_assignment(
- nodes_df=nodes_df,
- weighted_edges_df=weighted_edges_df,
- od_ss=od_ss,
- quarter_demand=quarter_demand,
- assigned_demand=assigned_demand,
- quarter_counts=quarter_counts,
- trip_info=trip_info,
- agent_time_limit=agent_time_limit,
- sample_interval=sample_interval,
- agents_path=agents_path,
- hour=hour,
- quarter=quarter,
- ss_id=ss_id,
- alpha_f=alpha_f,
- beta_f=beta_f,
- )
-
- od_residual_list += od_residual_ss_list
- # write_edge_vol(edges_df=edges_df,
- # simulation_outputs=simulation_outputs,
- # quarter=quarter,
- # hour=hour,
- # scen_nm='ss{}_{}'.format(ss_id, scen_nm))
- logging.info(
- f'HR {hour} QT {quarter} SS {ss_id}'
- ' finished, max vol '
- f'{np.max(edges_df["vol_true"])}, '
- f'time {time.time() - time_ss_0}'
- )
-
- # write quarterly results
- edges_df['vol_tot'] += edges_df['vol_true']
- if True: # hour >=16 or (hour==15 and quarter==3):
- write_edge_vol(
- edges_df=edges_df,
- simulation_outputs=simulation_outputs,
- quarter=quarter,
- hour=hour,
- scen_nm=scen_nm,
- )
-
- if hour % 3 == 0:
- trip_info_df = pd.DataFrame(
- [
- [
- trip_key[0],
- trip_key[1],
- trip_key[2],
- trip_value[0],
- trip_value[1],
- trip_value[2],
- trip_value[3],
- trip_value[4],
- trip_value[5],
- ]
- for trip_key, trip_value in trip_info.items()
- ],
- columns=[
- 'agent_id',
- 'origin_nid',
- 'destin_nid',
- 'travel_time',
- 'travel_time_used',
- 'stop_nid',
- 'stop_hour',
- 'stop_quarter',
- 'stop_ssid',
- ],
- )
- trip_info_df.to_csv(
- simulation_outputs + '/trip_info'
- f'/trip_info_{scen_nm}_hr{hour}'
- '.csv',
- index=False,
- )
-
- # output individual trip travel time and stop location
-
- trip_info_df = pd.DataFrame(
- [
- [
- trip_key[0],
- trip_key[1],
- trip_key[2],
- trip_value[0],
- trip_value[1],
- trip_value[2],
- trip_value[3],
- trip_value[4],
- trip_value[5],
- ]
- for trip_key, trip_value in trip_info.items()
- ],
- columns=[
- 'agent_id',
- 'origin_nid',
- 'destin_nid',
- 'travel_time',
- 'travel_time_used',
- 'stop_nid',
- 'stop_hour',
- 'stop_quarter',
- 'stop_ssid',
- ],
- )
- # Add the no path OD to the trip info
- trip_info_no_path = od_no_path.drop(
- columns=[
- col
- for col in od_no_path.columns
- if col not in ['agent_id', 'origin_nid', 'destin_nid']
- ]
- )
- trip_info_no_path['travel_time'] = 360000
- trip_info_no_path['travel_time_used'] = np.nan
- trip_info_no_path['stop_nid'] = np.nan
- trip_info_no_path['stop_hour'] = np.nan
- trip_info_no_path['stop_quarter'] = np.nan
- trip_info_no_path['stop_ssid'] = np.nan
- trip_info_df = pd.concat(
- [trip_info_df, trip_info_no_path], ignore_index=True
- )
-
- trip_info_df.to_csv(
- simulation_outputs + f'/trip_info/trip_info_{scen_nm}.csv',
- index=False,
- )
-
- write_final_vol(
- edges_df=edges_df,
- simulation_outputs=simulation_outputs,
- quarter=quarter,
- hour=hour,
- scen_nm=scen_nm,
- )
-
- network_edges = self.csv_filenames[0]
- network_nodes = self.csv_filenames[1]
- closed_edges_file = self.csv_filenames[2]
- demand_file = self.csv_filenames[3]
- simulation_outputs = self.simulation_outputs
- scen_nm = 'simulation_out'
-
- hour_list = list(range(6, 9))
- quarter_list = [0, 1, 2, 3, 4, 5]
- closure_hours = hour_list
-
- edges_df = pd.read_csv(network_edges)
- # edges_df = edges_df[["uniqueid", "geometry", "osmid", "length", "type",
- # "lanes", "maxspeed", "fft", "capacity",
- # "start_nid", "end_nid"]]
- edges_df = edges_df[
- [
- 'uniqueid',
- 'geometry',
- 'length',
- 'type',
- 'lanes',
- 'maxspeed',
- 'capacity',
- 'start_nid',
- 'end_nid',
- ]
- ]
- edges_df = gpd.GeoDataFrame(
- edges_df, crs='epsg:4326', geometry=edges_df['geometry'].map(loads)
- )
- # pay attention to the unit conversion, length is in meters, maxspeed is mph
- # fft is in seconds
- edges_df['fft'] = edges_df['length'] / edges_df['maxspeed'] * 2.23694
- edges_df = edges_df.sort_values(by='fft', ascending=False).drop_duplicates(
- subset=['start_nid', 'end_nid'], keep='first'
- )
- edges_df['edge_str'] = (
- edges_df['start_nid'].astype('str')
- + '-'
- + edges_df['end_nid'].astype('str')
- )
- edges_df['capacity'] = np.where(
- edges_df['capacity'] < 1, 950, edges_df['capacity']
- )
- edges_df['normal_capacity'] = edges_df['capacity']
- edges_df['normal_fft'] = edges_df['fft']
- edges_df['t_avg'] = edges_df['fft']
- edges_df['u'] = edges_df['start_nid']
- edges_df['v'] = edges_df['end_nid']
- edges_df = edges_df.set_index('edge_str')
- # closure locations
- closed_links = pd.read_csv(closed_edges_file)
- for row in closed_links.itertuples():
- edges_df.loc[(edges_df['uniqueid'] == row.uniqueid), 'capacity'] = 1
- edges_df.loc[(edges_df['uniqueid'] == row.uniqueid), 'fft'] = 36000
- # output closed file for visualization
- # edges_df.loc[edges_df['fft'] == 36000, ['uniqueid',
- # 'start_nid',
- # 'end_nid',
- # 'capacity',
- # 'fft',
- # 'geometry']].to_csv(
- # simulation_outputs +
- # '/closed_links_'
- # f'{scen_nm}.csv')
-
- # nodes processing
- nodes_df = pd.read_csv(network_nodes)
-
- nodes_df['x'] = nodes_df['lon']
- nodes_df['y'] = nodes_df['lat']
- nodes_df = nodes_df.set_index('node_id')
-
- # demand processing
- t_od_0 = time.time()
- od_all = pd.read_csv(demand_file)
- t_od_1 = time.time()
- logging.info('%d sec to read %d OD pairs', t_od_1 - t_od_0, od_all.shape[0])
- # run residual_demand_assignment
- assignment(
- edges_df=edges_df,
- nodes_df=nodes_df,
- od_all=od_all,
- simulation_outputs=simulation_outputs,
- scen_nm=scen_nm,
- hour_list=hour_list,
- quarter_list=quarter_list,
- closure_hours=closure_hours,
- closed_links=closed_links,
- )
-
- # @abstractmethod
- def update_od_file(
- self,
- old_nodes: str,
- old_det: str,
- new_nodes: str,
- new_det: str,
- od: str,
- origin_ids: list[int],
- ) -> pd.DataFrame:
- """
- Update origin-destination (OD) file from changes in population data.
-
- This function updates an OD file by calculating the population changes
- at each node between two time steps and generates trips originating
- from specified origins and ending at nodes where the population has
- increased. The updated OD data is saved to a new file.
-
- Args__
- old_nodes (str): Path to the CSV file containing the node
- information at the previous time step.
- old_det (str): Path to the JSON file containing the building
- information at the previous time step.
- new_nodes (str): Path to the CSV file containing the node
- information at the current time step.
- new_det (str): Path to the JSON file containing the building
- information at the current time step.
- od (str): Path to the existing OD file to be updated.
- origin_ids (List[int]): List of IDs representing possible origins
- for generating trips.
-
- Returns__
- pd.DataFrame: The updated OD DataFrame with new trips based on
- population changes.
-
- Raises__
- FileNotFoundError: If any of the provided file paths are incorrect
- or the files do not exist.
- json.JSONDecodeError: If the JSON files cannot be read or parsed
- correctly.
- KeyError: If expected keys are missing in the JSON data.
- ValueError: If there are issues with data conversions or
- calculations.
-
- Examples__
- >>> update_od_file(
- 'old_nodes.csv',
- 'old_building_data.json',
- 'new_nodes.csv',
- 'new_building_data.json',
- 'existing_od.csv',
- [1, 2, 3, 4]
- )
- The function will process the provided files, calculate population
- changes, and update the OD file with new trips. The result will be
- saved to 'updated_od.csv'.
-
- Notes__
- - Ensure that the columns `lat`, `lon`, and `node_id` are present
- in the nodes CSV files.
- - The `det` files should contain the `Buildings` and
- `GeneralInformation` sections with appropriate keys.
- - The OD file should have the columns `agent_id`, `origin_nid`,
- `destin_nid`, `hour`, and `quarter`.
- """
-
- # Extract the building information from the det file and convert it to
- # a pandas dataframe
- def extract_building_from_det(det):
- # Open the det file
- with Path.open(det, encoding='utf-8') as file:
- # Return the JSON object as a dictionary
- json_data = json.load(file)
-
- # Extract the required information and convert it to a pandas
- # dataframe
- extracted_data = []
-
- for aim_id, info in json_data['Buildings']['Building'].items():
- general_info = info.get('GeneralInformation', {})
- extracted_data.append(
- {
- 'AIM_id': aim_id,
- 'Latitude': general_info.get('Latitude'),
- 'Longitude': general_info.get('Longitude'),
- 'Population': general_info.get('Population'),
- }
- )
-
- return pd.DataFrame(extracted_data)
-
- # Aggregate the population in buildings to the closest road network
- # node
- def closest_neighbour(building_df, nodes_df):
- # Find the nearest road network node to each building
- nodes_xy = np.array(
- [nodes_df['lat'].to_numpy(), nodes_df['lon'].to_numpy()]
- ).transpose()
- building_df['closest_node'] = building_df.apply(
- lambda x: cdist(
- [(x['Latitude'], x['Longitude'])], nodes_xy
- ).argmin(),
- axis=1,
- )
-
- # Merge the road network and building dataframes
- merged_df = nodes_df.merge(
- building_df, left_on='node_id', right_on='closest_node', how='left'
- )
- merged_df = merged_df.drop(
- columns=['AIM_id', 'Latitude', 'Longitude', 'closest_node']
- )
- merged_df = merged_df.fillna(0)
-
- # Aggregate population of neareast buildings to the road network
- # node
- updated_nodes_df = (
- merged_df.groupby('node_id')
- .agg(
- {
- 'lon': 'first',
- 'lat': 'first',
- 'geometry': 'first',
- 'Population': 'sum',
- }
- )
- .reset_index()
- )
-
- return updated_nodes_df # noqa: RET504
-
- # Function to add the population information to the nodes file
- def update_nodes_file(nodes, det):
- # Read the nodes file
- nodes_df = pd.read_csv(nodes)
- # Extract the building information from the det file and convert it
- # to a pandas dataframe
- building_df = extract_building_from_det(det)
- # Aggregate the population in buildings to the closest road network
- # node
- updated_nodes_df = closest_neighbour(building_df, nodes_df)
-
- return updated_nodes_df # noqa: RET504
-
- # Read the od file
- od_df = pd.read_csv(od)
- # Add the population information to nodes dataframes at the last and
- # current time steps
- old_nodes_df = update_nodes_file(old_nodes, old_det)
- new_nodes_df = update_nodes_file(new_nodes, new_det)
- # Calculate the population changes at each node (assuming that
- # population will only increase at each node)
- population_change_df = old_nodes_df.copy()
- population_change_df['Population Change'] = (
- new_nodes_df['Population'] - old_nodes_df['Population']
- )
- population_change_df['Population Change'].astype(int)
- # Randomly generate the trips that start at one of the connections
- # between Alameda Island and Oakland, and end at the nodes where
- # population increases and append them to the od file
- # Generate OD data for each node with increased population and append
- # it to the od file
- for _, row in population_change_df.iterrows():
- # Only process rows with positive population difference
- if row['Population_Difference'] > 0:
- for _ in range(row['Population_Difference']):
- # Generate random origin_nid
- origin_nid = np.random.choice(origin_ids)
- # Set the destin_nid to the node_id of the increased
- # population
- destin_nid = row['node_id']
- # Generate random hour and quarter
- hour = np.random.randint(5, 24)
- quarter = np.random.randint(0, 5)
- # Append to od dataframe
- od_df = od_df.append(
- {
- 'agent_id': 0,
- 'origin_nid': origin_nid,
- 'destin_nid': destin_nid,
- 'hour': hour,
- 'quarter': quarter,
- },
- ignore_index=True,
- )
- od_df.to_csv('updated_od.csv')
- return od_df
-
- # def get_graph_network(self,
- # det_file_path: str,
- # output_files: Optional[List[str]] = None,
- # save_additional_attributes:
- # Literal['', 'residual_demand'] = ''
- # ) -> Tuple[Dict, Dict]:
- # """
- # Create a graph network from inventory data .
-
- # This function processes inventory files containing road and structural
- # data, constructs a graph network with nodes and edges, and optionally
- # saves additional attributes such as residual demand. The node and edge
- # features are saved to specified output files.
-
- # Args__
- # det_file_path (str): The path to the deterministic result JSON
- # output_files (Optional[List[str]]): A list of file paths where the
- # node and edge features will be saved. The first file in the
- # list is used for edges and the second for nodes.
- # save_additional_attributes (Literal['', 'residual_demand']):
- # A flag indicating whether to save additional attributes.
- # The only supported additional attribute is 'residual_demand'.
-
- # Returns__
- # Tuple[Dict, Dict]: A tuple containing two dictionaries:
- # - The first dictionary contains the edge features.
- # - The second dictionary contains the node features.
-
- # Example__
- # >>> det_file_path = 'path/to/Results_det.json'
- # >>> output_files = ['edges.csv', 'nodes.csv']
- # >>> save_additional_attributes = 'residual_demand'
- # >>> edges, nodes = get_graph_network(det_file_path,
- # output_files,
- # save_additional_attributes)
- # >>> print(edges)
- # >>> print(nodes)
- # """
- # # if output_files is None:
- # # self.graph_network['output_files'] = ['edges.csv', 'nodes.csv']
-
- # def create_circle_from_lat_lon(lat, lon, radius_ft, num_points=100):
- # """
- # Create a circle polygon from latitude and longitude.
-
- # Args__
- # lat (float): Latitude of the center.
- # lon (float): Longitude of the center.
- # radius_ft (float): Radius of the circle in feet.
- # num_points (int): Number of points to approximate the circle.
-
- # Returns__
- # Polygon: A Shapely polygon representing the circle.
- # """
- # # Earth's radius in kilometers
- # R_EARTH_FT = 20925721.784777
-
- # # Convert radius from kilometers to radians
- # radius_rad = radius_ft / R_EARTH_FT
-
- # # Convert latitude and longitude to radians
- # lat_rad = np.radians(lat)
- # lon_rad = np.radians(lon)
-
- # # Generate points around the circle
- # angle = np.linspace(0, 2 * np.pi, num_points)
- # lat_points = lat_rad + radius_rad * np.cos(angle)
- # lon_points = lon_rad + radius_rad * np.sin(angle)
-
- # # Convert radians back to degrees
- # lat_points_deg = np.degrees(lat_points)
- # lon_points_deg = np.degrees(lon_points)
-
- # # Create a polygon from the points
- # points = list(zip(lon_points_deg, lat_points_deg))
- # return Polygon(points)
-
- # def parse_element_geometries_from_geojson(
- # det_file: str,
- # save_additional_attributes: str
- # ) -> Tuple[List[Dict], Tuple[List[LineString], List[Dict]]]:
- # """
- # Parse geometries from R2D deterministic result json.
-
- # This function reads json files specified in `output_files`. It
- # extracts and parses `LineString` geometries from files that contain
- # "road" in their name. For other files, it accumulates point
- # features.
-
- # Args__
- # det_file (str): Path to the R2D deterministic result Json.
-
- # Returns__
- # tuple: A tuple containing two elements:
- # - ptdata (list of dict): List of point features parsed from
- # GeoJSON files for bridges and tunnels
- # - road_data (tuple): A tuple containing:
- # - road_polys (list of LineString): List of `LineString`
- # objects created from the road geometries in GeoJSON
- # files that contain "road" in their name.
- # - road_features (list of dict): List of road features
- # as dictionaries from GeoJSON files that contain
- # "roads" in their name.
-
- # Raises__
- # FileNotFoundError: If any of the specified files in
- # `output_files` do not exist.
- # json.JSONDecodeError: If a file cannot be parsed as valid JSON.
- # """
- # ptdata = {}
- # with open(det_file, "r", encoding="utf-8") as file:
- # # Read the JSON file:
- # temp = json.load(file)
- # if 'TransportationNetwork' not in temp:
- # raise KeyError(
- # 'The deterministic result JSON file does not contain TransportationNetwork')
- # temp = temp['TransportationNetwork']
- # # If the file contains road information:
- # if 'Roadway' in temp:
- # # Extract road features:
- # road_features = temp['Roadway']
- # # Read road polygons, add asset type information to
- # # road features and, if required, parse existing road
- # # attributes to infer information required for network
- # # analysis:
- # road_polys = {}
- # for AIM_ID, item in road_features.items():
- # road_polys.update({AIM_ID: loads(item["GeneralInformation"]["geometry"])})
- # # road_features[ind]['properties']['asset_type'] = \
- # # 'road'
- # if save_additional_attributes:
- # # Access properties for element at index ind:
- # properties = item['GeneralInformation']
- # # Get MTFCC value:
- # mtfcc = properties['MTFCC']
- # # Update road features for the indexed element
- # # with number of lanes, maximum speed, and
- # # road capacity:
- # properties.update({
- # 'lanes':
- # self.attribute_maps['lanes'][mtfcc],
- # 'maxspeed':
- # self.attribute_maps['speed'][mtfcc],
- # 'capacity':
- # self.attribute_maps['capacity'][
- # properties['lanes']]
- # })
- # if 'Bridge' in temp:
- # ptdata['Bridge'] = temp['Bridge']
- # # for feature in temp["features"]:
- # # feature['properties']['asset_type'] = asset_type
- # # ptdata += temp["features"]
- # if 'Tunnel' in temp:
- # ptdata['Tunnel'] = temp['Tunnel']
-
- # return ptdata, (road_polys, road_features)
-
- # def find_intersections(lines: List[LineString]) -> Set[Point]:
- # """
- # Find intersection points between pairs of LineString geometries.
-
- # This function takes a list of `LineString` objects and identifies
- # points where any two lines intersect. The intersections are
- # returned as a set of `Point` objects.
-
- # Args__
- # lines (List[LineString]): A list of `LineString` objects. The
- # function computes intersections between each pair of
- # `LineString` objects in this list.
-
- # Returns__
- # Set[Point]: A set of `Point` objects representing the
- # intersection points between the `LineString` objects. If
- # multiple intersections occur at the same location, it will
- # only be included once in the set.
-
- # Example__
- # >>> from shapely.geometry import LineString
- # >>> line1 = LineString([(0, 0), (1, 1)])
- # >>> line2 = LineString([(0, 1), (1, 0)])
- # >>> find_intersections([line1, line2])
- # {}
-
- # Notes__
- # - The function assumes that all input geometries are valid
- # `LineString`
- # objects.
- # - The resulting set may be empty if no intersections are found.
-
- # Raises__
- # TypeError: If any element in `lines` is not a `LineString`
- # object.
- # """
- # intersections = set()
- # for i, line1 in enumerate(lines):
- # for line2 in lines[i + 1:]:
- # if line1.intersects(line2):
- # inter_points = line1.intersection(line2)
- # if inter_points.geom_type == "Point":
- # intersections.add(inter_points)
- # elif inter_points.geom_type == "MultiPoint":
- # intersections.update(inter_points.geoms)
- # return intersections
-
- # def cut_lines_at_intersections(lines: List[LineString],
- # line_features: List[Dict],
- # intersections: List[Point]
- # ) -> List[LineString]:
- # """
- # Cut LineStrings at intersection points & return resulting segments.
-
- # This function takes a list of `LineString` objects and a list of
- # `Point` objects representing intersection points. For each
- # `LineString`, it splits the line at each intersection point. The
- # resulting segments are collected and returned.
-
- # Args__
- # lines (List[LineString]): A list of `LineString` objects to be
- # cut at the intersection points.
- # line_features (List[Dict]): List of features for the
- # `LineString` objects in lines.
- # intersections (List[Point]): A list of `Point` objects where
- # the lines are intersected and split.
-
- # Returns__
- # List[LineString]: A list of `LineString` objects resulting from
- # cutting the original lines at the
- # intersection points.
-
- # Notes__
- # - The `split` function from `shapely.ops` is used to perform
- # the cutting of lines at intersection points.
- # - The function handles cases where splitting results in a
- # `GeometryCollection` by extracting only `LineString`
- # geometries.
-
- # Example__
- # >>> from shapely.geometry import LineString, Point
- # >>> from shapely.ops import split
- # >>> lines = [
- # ... LineString([(0, 0), (2, 2)]),
- # ... LineString([(2, 0), (0, 2)])
- # ... ]
- # >>> intersections = [
- # ... Point(1, 1)
- # ... ]
- # >>> cut_lines_at_intersections(lines, intersections)
- # [,
- # ]
- # """
- # new_lines = []
- # new_line_features = []
-
- # for ind_line, line in enumerate(lines):
- # segments = [line] # Start with the original line
- # for point in intersections:
- # new_segments = []
- # features = []
- # for segment in segments:
- # if segment.intersects(point):
- # split_result = split(segment, point)
- # if split_result.geom_type == "GeometryCollection":
- # new_segments.extend(
- # geom
- # for geom in split_result.geoms
- # if geom.geom_type == "LineString"
- # )
- # features.extend([copy.deepcopy(
- # line_features[ind_line])
- # for _ in range(
- # len(
- # split_result.geoms
- # ))])
- # elif split_result.geom_type == "LineString":
- # segments.append(split_result)
- # features.append(line_features[ind_line])
- # else:
- # new_segments.append(segment)
- # features.append(line_features[ind_line])
- # segments = new_segments
-
- # # Add remaining segments that were not intersected by any
- # # points:
- # new_lines.extend(segments)
- # new_line_features.extend(features)
-
- # return (new_lines, new_line_features)
-
- # def save_cut_lines_and_intersections(lines: List[LineString],
- # points: List[Point]) -> None:
- # """
- # Save LineString and Point objects to separate GeoJSON files.
-
- # This function converts lists of `LineString` and `Point` objects to
- # GeoJSON format and saves them to separate files. The `LineString`
- # objects are saved to "lines.geojson", and the `Point` objects are
- # saved to "points.geojson".
-
- # Args__
- # lines (List[LineString]): A list of `LineString` objects to be
- # saved in GeoJSON format.
- # intersections (List[Point]): A list of `Point` objects to be
- # saved in GeoJSON format.
-
- # Returns__
- # None: This function does not return any value. It writes
- # GeoJSON data to files.
-
- # Notes__
- # - The function uses the `mapping` function from
- # `shapely.geometry` to convert geometries to GeoJSON format.
- # - Two separate GeoJSON files are created: one for lines and one
- # for points.
- # - The output files are named "lines.geojson" and
- # "points.geojson" respectively.
-
- # Example__
- # >>> from shapely.geometry import LineString, Point
- # >>> lines = [
- # ... LineString([(0, 0), (1, 1)]),
- # ... LineString([(1, 1), (2, 2)])
- # ... ]
- # >>> points = [
- # ... Point(0.5, 0.5),
- # ... Point(1.5, 1.5)
- # ... ]
- # >>> save_cut_lines_and_intersections(lines, points)
- # # This will create "lines.geojson" and "points.geojson" files
- # with the corresponding GeoJSON data.
- # """
- # # Convert LineString objects to GeoJSON format
- # features = []
- # for line in lines:
- # features.append(
- # {"type": "Feature", "geometry": mapping(
- # line), "properties": {}}
- # )
-
- # # Create a GeoJSON FeatureCollection
- # geojson = {"type": "FeatureCollection", "features": features}
-
- # # Save the GeoJSON to a file
- # with open("lines.geojson", "w", encoding="utf-8") as file:
- # json.dump(geojson, file, indent=2)
-
- # # Convert Point objects to GeoJSON format
- # features = []
- # for point in points:
- # features.append(
- # {"type": "Feature", "geometry": mapping(
- # point), "properties": {}}
- # )
-
- # # Create a GeoJSON FeatureCollection
- # geojson = {"type": "FeatureCollection", "features": features}
-
- # # Save the GeoJSON to a file
- # with open("points.geojson", "w", encoding="utf-8") as file:
- # json.dump(geojson, file, indent=2)
-
- # def find_repeated_line_pairs(lines: List[LineString]) -> Set[Tuple]:
- # """
- # Find and groups indices of repeated LineString objects from a list.
-
- # This function processes a list of `LineString` objects to identify
- # and group all unique index pairs where LineString objects are
- # repeated. The function converts each `LineString` to its
- # Well-Known Text (WKT) representation to identify duplicates.
-
- # Args__
- # lines (List[LineString]): A list of `LineString` objects to be
- # analyzed for duplicates.
-
- # Returns__
- # Set[Tuple]: A set of tuples, where each tuple contains indices
- # for LineString objects that are repeated.
-
- # Raises__
- # TypeError: If any element in the `lines` list is not an
- # instance of `LineString`.
-
- # Example__
- # >>> from shapely.geometry import LineString
- # >>> lines = [
- # ... LineString([(0, 0), (1, 1)]),
- # ... LineString([(0, 0), (1, 1)]),
- # ... LineString([(1, 1), (2, 2)]),
- # ... LineString([(0, 0), (1, 1)]),
- # ... LineString([(1, 1), (2, 2)])
- # ... ]
- # >>> find_repeated_line_pairs(lines)
- # [{0, 1, 3}, {2, 4}]
- # """
- # line_indices = defaultdict(list)
-
- # for id, line in lines.items():
- # if not isinstance(line, LineString):
- # raise TypeError(
- # 'All elements in the input list must be LineString'
- # ' objects.')
-
- # # Convert LineString to its WKT representation to use as a
- # # unique identifier:
- # line_wkt = line.wkt
- # line_indices[line_wkt].append(id)
-
- # repeated_pairs = set()
- # for indices in line_indices.values():
- # if len(indices) > 1:
- # # Create pairs of all indices for the repeated LineString
- # for i, _ in enumerate(indices):
- # for j in range(i + 1, len(indices)):
- # repeated_pairs.add((indices[i], indices[j]))
-
- # repeated_pairs = list(repeated_pairs)
- # ind_matched = []
- # repeated_polys = []
- # for index_p1, pair1 in enumerate(repeated_pairs):
- # pt1 = set(pair1)
- # for index_p2, pair2 in enumerate(repeated_pairs[index_p1+1:]):
- # if (index_p1 + index_p2 + 1) not in ind_matched:
- # pt2 = set(pair2)
- # if bool(pt1 & pt2):
- # pt1 |= pt2
- # ind_matched.append(index_p1 + index_p2 + 1)
- # if pt1 not in repeated_polys and index_p1 not in ind_matched:
- # repeated_polys.append(pt1)
-
- # return repeated_polys
-
- # def match_edges_to_points(ptdata: List[Dict],
- # road_polys: List[LineString],
- # road_features: List[Dict]
- # ) -> List[List[int]]:
- # """
- # Match points to closest road polylines based on name similarity.
-
- # This function takes a list of points and a list of road polylines.
- # For each point, it searches for intersecting road polylines within
- # a specified radius and calculates the similarity between the
- # point's associated facility name and the road's name. It returns a
- # list of lists where each sublist contains indices of the road
- # polylines that best match the point based on the similarity score.
-
- # Args__
- # ptdata (List[Dict[str, Any]]): A list of dictionaries where
- # each dictionary represents a point with its geometry and
- # properties. The 'geometry' key should contain 'coordinates'
- # , and the 'properties' key should contain 'tunnel_name' or
- # 'facility_carried'.
- # road_polys (List[LineString]): A list of `LineString` objects
- # representing road polylines.
-
- # Returns__
- # List[List[int]]: A list of lists where each sublist contains
- # the indices of road polylines that have the highest textual
- # similarity to the point's facility name. If no similarity
- # is found, the sublist is empty.
-
- # Notes__
- # - The function uses a search radius of 100 feet to find
- # intersecting road polylines.
- # - TF-IDF vectors are used to compute the textual similarity
- # between the facility names and road names.
- # - Cosine similarity is used to determine the best matches.
-
- # Example__
- # >>> from shapely.geometry import Point, LineString
- # >>> ptdata = [
- # ... {"geometry": {"coordinates": [1.0, 1.0]},
- # "properties": {"tunnel_name": "Tunnel A"}},
- # ... {"geometry": {"coordinates": [2.0, 2.0]},
- # "properties": {"facility_carried": "Road B"}}
- # ... ]
- # >>> road_polys = [
- # ... LineString([(0, 0), (2, 2)]),
- # ... LineString([(1, 1), (3, 3)])
- # ... ]
- # >>> match_edges_to_points(ptdata, road_polys)
- # [[0], [1]]
- # """
- # edge_matches = []
- # for ptdata_type, type_data in ptdata.items():
- # for AIM_ID, point in type_data.items():
- # lon = point["GeneralInformation"]["location"]["longitude"]
- # lat = point["GeneralInformation"]["location"]["latitude"]
- # search_radius = create_circle_from_lat_lon(lat, lon, 100)
- # # Check for intersections:
- # intersecting_polylines = [
- # AIM_ID
- # for (AIM_ID, poly) in road_polys.items()
- # if poly.intersects(search_radius)
- # ]
- # properties = point["GeneralInformation"]
- # if ptdata_type == "Tunnel":
- # facility = properties.get("tunnel_name", "")
- # elif ptdata_type == "Bridge":
- # # facility = properties["facility_carried"]
- # facility = properties.get("facility_carried", "")
- # max_similarity = 0
- # similarities = {}
- # for AIM_ID in intersecting_polylines:
- # roadway = road_features[AIM_ID]["GeneralInformation"].get("NAME", None)
- # if roadway:
- # # Create TF-IDF vectors:
- # vectorizer = TfidfVectorizer()
- # tfidf_matrix = vectorizer.fit_transform(
- # [facility, roadway.lower()]
- # )
-
- # # Calculate cosine similarity:
- # similarity = cosine_similarity(
- # tfidf_matrix[0:1], tfidf_matrix[1:2]
- # )
- # else:
- # similarity = -1
- # similarities.update({AIM_ID:similarity})
- # if similarity > max_similarity:
- # max_similarity = similarity
- # if max_similarity == 0:
- # max_similarity = -1
- # indices_of_max = [
- # intersecting_polylines[index]
- # for index, value in (similarities).items()
- # if value == max_similarity
- # ]
- # edge_matches.append(indices_of_max)
-
- # return edge_matches
-
- # def merge_brtn_features(ptdata: List[Dict],
- # road_features_brtn: List[Dict],
- # edge_matches: List[List],
- # save_additional_attributes: str) -> List[Dict]:
- # """
- # Merge bridge/tunnel features to road features using edge matches.
-
- # This function updates road features with additional attributes
- # derived from bridge or tunnel point data. It uses edge matches to
- # determine how to distribute lane and capacity information among
- # road features.
-
- # Args__
- # ptdata (List[Dict]): A list of dictionaries where each
- # dictionary contains properties of bridge or tunnel
- # features. Each dictionary should have 'asset_type',
- # 'traffic_lanes_on', 'structure_number',
- # 'total_number_of_lanes', and 'tunnel_number' as keys.
- # road_features_brtn (List[Dict]): A list of dictionaries
- # representing road features that will be updated. Each
- # dictionary should have a 'properties' key where attributes
- # are stored.
- # edge_matches (List[List[int]]): A list of lists, where each
- # sublist contains indices that correspond to
- # `road_features_brtn` and represent which features should be
- # updated together.
- # save_additional_attributes (str): A flag indicating whether to
- # save additional attributes like 'lanes' and 'capacity'. If
- # non-empty, additional attributes will be saved.
-
- # Returns__
- # List[Dict]: The updated list of road features with merged
- # attributes.
-
- # Example__
- # >>> ptdata = [
- # ... {'properties': {'asset_type': 'bridge',
- # 'traffic_lanes_on': 4,
- # 'structure_number': '12345'}},
- # ... {'properties': {'asset_type': 'tunnel',
- # 'total_number_of_lanes': 6,
- # 'tunnel_number': '67890'}}
- # ... ]
- # # List of empty
- # >>> road_features_brtn = [{} for _ in range(4)]
- # dictionaries for demonstration
- # >>> edge_matches = [[0, 1], [2, 3]]
- # >>> save_additional_attributes = 'yes'
- # >>> updated_features = merge_brtn_features(
- # ptdata,
- # road_features_brtn,
- # edge_matches,
- # save_additional_attributes)
- # >>> print(updated_features)
- # """
- # poly_index = 0
- # for item_index, edge_indices in enumerate(edge_matches):
- # nitems = len(edge_indices)
- # features = ptdata[item_index]['properties']
- # asset_type = features['asset_type']
- # if asset_type == 'bridge':
- # total_nlanes = features['traffic_lanes_on']
- # struct_no = features['structure_number']
- # else:
- # total_nlanes = features['total_number_of_lanes']
- # struct_no = features['tunnel_number']
-
- # lanes_per_item = round(int(total_nlanes)/nitems)
-
- # for index in range(poly_index, poly_index + nitems):
- # properties = road_features_brtn[index]['properties']
- # properties['asset_type'] = asset_type
- # properties['OID'] = struct_no
- # if save_additional_attributes:
- # properties['lanes'] = lanes_per_item
- # properties['capacity'] = calculate_road_capacity(
- # lanes_per_item)
-
- # poly_index += nitems
- # return road_features_brtn
-
- # def get_nodes_edges(lines: List[LineString],
- # length_unit: Literal['ft', 'm'] = 'ft'
- # ) -> Tuple[Dict, Dict]:
- # """
- # Extract nodes and edges from a list of LineString objects.
-
- # This function processes a list of `LineString` geometries to
- # generate nodes and edges. Nodes are defined by their unique
- # coordinates, and edges are represented by their start and end
- # nodes, length, and geometry.
-
- # Args__
- # lines (List[LineString]): A list of `LineString` objects
- # representing road segments.
- # length_unit (Literal['ft', 'm']): The unit of length for the
- # edge distances. Defaults to 'ft'. Acceptable values are
- # 'ft' for feet and 'm' for meters.
-
- # Returns__
- # Tuple[Dict[int, Dict[str, float]], Dict[int, Dict[str, Any]]]:
- # - A dictionary where keys are node IDs and values are
- # dictionaries with node attributes:
- # - 'lon': Longitude of the node.
- # - 'lat': Latitude of the node.
- # - 'geometry': WKT representation of the node.
- # - A dictionary where keys are edge IDs and values are
- # dictionaries with edge attributes:
- # - 'start_nid': ID of the start node.
- # - 'end_nid': ID of the end node.
- # - 'length': Length of the edge in the specified unit.
- # - 'geometry': WKT representation of the edge.
-
- # Raises__
- # TypeError: If any element in the `lines` list is not an
- # instance of `LineString`.
-
- # Example__
- # >>> from shapely.geometry import LineString
- # >>> lines = [
- # ... LineString([(0, 0), (1, 1)]),
- # ... LineString([(1, 1), (2, 2)])
- # ... ]
- # >>> nodes, edges = get_nodes_edges(lines, length_unit='m')
- # >>> print(nodes)
- # >>> print(edges)
- # """
- # node_list = []
- # edges = {}
- # node_counter = 0
- # for line_counter, line in enumerate(lines):
- # # Ensure the object is a LineString
- # if not isinstance(line, LineString):
- # raise TypeError(
- # "All elements in the list must be LineString objects.")
-
- # # Extract coordinates
- # coords = list(line.coords)
- # ncoord_pairs = len(coords)
- # if ncoord_pairs > 0:
- # start_node_coord = coords[0]
- # end_node_coord = coords[-1]
-
- # if start_node_coord not in node_list:
- # node_list.append(start_node_coord)
- # start_nid = node_counter
- # node_counter += 1
- # else:
- # start_nid = node_list.index(start_node_coord)
-
- # if end_node_coord not in node_list:
- # node_list.append(end_node_coord)
- # end_nid = node_counter
- # node_counter += 1
- # else:
- # end_nid = node_list.index(end_node_coord)
-
- # length = 0
- # (lon, lat) = line.coords.xy
- # for pair_no in range(ncoord_pairs - 1):
- # length += haversine_dist([lat[pair_no], lon[pair_no]],
- # [lat[pair_no+1],
- # lon[pair_no+1]])
- # if length_unit == 'm':
- # length = 0.3048*length
-
- # edges[line_counter] = {'start_nid': start_nid,
- # 'end_nid': end_nid,
- # 'length': length,
- # 'geometry': line.wkt}
-
- # nodes = {}
- # for node_id, node_coords in enumerate(node_list):
- # nodes[node_id] = {'lon': node_coords[0],
- # 'lat': node_coords[1],
- # 'geometry': 'POINT ('
- # f'{node_coords[0]:.7f} '
- # f'{node_coords[1]:.7f})'}
-
- # return (nodes, edges)
-
- # def get_node_edge_features(updated_road_polys: List[LineString],
- # updated_road_features: List[Dict],
- # output_files: List[str]
- # ) -> Tuple[Dict, Dict]:
- # """
- # Extract/write node and edge features from updated road polygons.
-
- # This function processes road polygon data to generate nodes and
- # edges, then writes the extracted features to specified output
- # files.
-
- # Args__
- # updated_road_polys (List[LineString]): A list of LineString
- # objects representing updated road polygons.
- # updated_road_features (List[Dict]): A list of dictionaries
- # containing feature properties for each road segment.
- # output_files (List[str]): A list of two file paths where the
- # first path is for edge data and the second for node data.
-
- # Returns__
- # Tuple[Dict, Dict]: A tuple containing two dictionaries:
- # - The first dictionary contains edge data.
- # - The second dictionary contains node data.
-
- # Raises__
- # TypeError: If any input is not of the expected type.
-
- # Example__
- # >>> from shapely.geometry import LineString
- # >>> road_polys = [LineString([(0, 0), (1, 1)]),
- # LineString([(1, 1), (2, 2)])]
- # >>> road_features = [{'properties': {'OID': 1,
- # 'asset_type': 'road',
- # 'lanes': 2,
- # 'capacity': 2000,
- # 'maxspeed': 30}}]
- # >>> output_files = ['edges.csv', 'nodes.csv']
- # >>> get_node_edge_features(road_polys,
- # road_features,
- # output_files)
- # """
- # (nodes, edges) = get_nodes_edges(
- # updated_road_polys, length_unit='m')
-
- # with open(output_files[1], 'w', encoding="utf-8") as nodes_file:
- # nodes_file.write('node_id, lon, lat, geometry\n')
- # for key in nodes:
- # nodes_file.write(f'{key}, {nodes[key]["lon"]}, '
- # f'{nodes[key]["lat"]}, '
- # f'{nodes[key]["geometry"]}\n')
-
- # with open(output_files[0], 'w', encoding="utf-8") as edge_file:
- # edge_file.write('uniqueid, start_nid, end_nid, osmid, length, '
- # 'type, lanes, maxspeed, capacity, fft, '
- # 'geometry\n')
-
- # for key in edges:
- # features = updated_road_features[key]['properties']
- # edges[key]['osmid'] = features['OID']
- # edges[key]['type'] = features['asset_type']
- # edges[key]['lanes'] = features['lanes']
- # edges[key]['capacity'] = features['capacity']
- # maxspeed = features['maxspeed']
- # edges[key]['maxspeed'] = maxspeed
- # free_flow_time = edges[key]['length'] / \
- # (maxspeed*1609.34/3600)
- # edges[key]['fft'] = free_flow_time
- # edge_file.write(f'{key}, {edges[key]["start_nid"]}, '
- # f'{edges[key]["end_nid"]}, '
- # f'{edges[key]["osmid"]}, '
- # f'{edges[key]["length"]}, '
- # f'{edges[key]["type"]}, '
- # f'{edges[key]["lanes"]}, {maxspeed}, '
- # f'{edges[key]["capacity"]}, '
- # f'{free_flow_time}, '
- # f'{edges[key]["geometry"]}\n')
-
- # return (edges, nodes)
-
- # print('Getting graph network for elements in '
- # f'{det_file_path}...')
-
- # # Read inventory data:
- # ptdata, (road_polys, road_features) = \
- # parse_element_geometries_from_geojson(det_file_path,
- # save_additional_attributes=save_additional_attributes)
-
- # # Find edges that match bridges and tunnels:
- # edge_matches = match_edges_to_points(ptdata, road_polys, road_features)
-
- # # Get the indices for bridges and tunnels:
- # brtn_poly_idx = [item for sublist in edge_matches for item in sublist]
-
- # # Detect repeated edges and save their indices:
- # repeated_edges = find_repeated_line_pairs(road_polys)
-
- # edges_remove = []
- # for edge_set in repeated_edges:
- # bridge_poly = set(brtn_poly_idx)
- # if edge_set & bridge_poly:
- # remove = list(edge_set - bridge_poly)
- # else:
- # temp = list(edge_set)
- # remove = temp[1:].copy()
- # edges_remove.extend(remove)
-
- # # Save polygons that are not bridge or tunnel edges or marked for
- # # removal in road polygons:
- # road_polys_local = [poly for (ind, poly) in road_polys.items() if
- # ind not in brtn_poly_idx + edges_remove]
- # road_features_local = [feature for (ind, feature) in
- # road_features.items()
- # if ind not in brtn_poly_idx + edges_remove]
-
- # # Save polygons that are not bridge or tunnel edges:
- # road_polys_brtn = [poly for (ind, poly) in enumerate(road_polys)
- # if ind in brtn_poly_idx]
- # road_features_brtn = [feature for (ind, feature)
- # in enumerate(road_features)
- # if ind in brtn_poly_idx]
- # road_features_brtn = merge_brtn_features(ptdata,
- # road_features_brtn,
- # edge_matches,
- # save_additional_attributes)
-
- # # Compute the intersections of road polygons:
- # intersections = find_intersections(road_polys_local)
-
- # # Cut road polygons at intersection points:
- # cut_lines, cut_features = cut_lines_at_intersections(
- # road_polys_local,
- # road_features_local,
- # intersections)
- # # Come back and update cut_lines_at_intersections to not intersect
- # # lines within a certain diameter of a bridge point
-
- # # Combine all polygons and their features:
- # updated_road_polys = cut_lines + road_polys_brtn
- # updated_road_features = cut_features + road_features_brtn
-
- # # Save created polygons (for debugging only)
- # # save_cut_lines_and_intersections(updated_road_polys, intersections)
-
- # # Get nodes and edges of the final set of road polygons:
- # (edges, nodes) = get_node_edge_features(updated_road_polys,
- # updated_road_features,
- # output_files)
- # self.graph_network['edges'] = edges
- # self.graph_network['nodes'] = nodes
-
- # print('Edges and nodes of the graph network are saved in '
- # f'{", ".join(output_files)}')
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/CMakeLists.txt b/modules/performRegionalEventSimulation/regionalGroundMotion/CMakeLists.txt
index 120def497..3c2d9e20a 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/CMakeLists.txt
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/CMakeLists.txt
@@ -5,7 +5,8 @@ simcenter_add_python_script(SCRIPT ComputeIntensityMeasure.py)
simcenter_add_python_script(SCRIPT CreateScenario.py)
simcenter_add_python_script(SCRIPT CreateStation.py)
simcenter_add_python_script(SCRIPT FetchOpenQuake.py)
-simcenter_add_python_script(SCRIPT FetchOpenSHA_25_4.py)
+simcenter_add_python_script(SCRIPT FetchOpenSHA.py)
+simcenter_add_python_script(SCRIPT FetchOpenSHA_1.5.2.py)
simcenter_add_python_script(SCRIPT HazardSimulation.py DEPENDS database)
simcenter_add_python_script(SCRIPT HazardSimulationEQ.py DEPENDS database)
simcenter_add_python_script(SCRIPT SelectGroundMotion.py)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py b/modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py
index 0664c625b..4d0e9c881 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/ComputeIntensityMeasure.py
@@ -101,7 +101,6 @@
import json # noqa: E402
import os # noqa: E402
import socket # noqa: E402
-import sys # noqa: E402
import time # noqa: E402
from pathlib import Path # noqa: E402
@@ -111,8 +110,14 @@
if 'stampede2' not in socket.gethostname():
from FetchOpenQuake import get_site_rup_info_oq
- #from FetchOpenSHA import * # noqa: F403
- from FetchOpenSHA_25_4 import * # noqa: F403
+ from FetchOpenSHA import (
+ getERF,
+ get_IM,
+ get_PointSource_info_CY2014,
+ get_rupture_info_ASK2014_aftershock,
+ get_rupture_info_CY2014,
+ get_site_prop,
+ )
import threading # noqa: E402
import ujson # noqa: E402
@@ -808,7 +813,7 @@ def compute_im( # noqa: C901, D103
ho_period = generator_info['Parameters'].get('Period')
if im_info['Type'] == 'Vector':
if im_info.get('SA') is None:
- sys.exit(
+ raise ValueError(
'SA is used in hazard downsampling but not defined in the intensity measure tab'
)
elif im_info.get('Periods', None) is not None:
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/CreateScenario.py b/modules/performRegionalEventSimulation/regionalGroundMotion/CreateScenario.py
index a7264b4a8..b0106de74 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/CreateScenario.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/CreateScenario.py
@@ -37,19 +37,23 @@
# Kuanshi Zhong
#
-import ujson as python_json
+import ujson as json
import os
import random
import socket
-import sys
import time
import numpy as np
import pandas as pd
if 'stampede2' not in socket.gethostname():
- #from FetchOpenSHA import * # noqa: F403
- from FetchOpenSHA_25_4 import * # noqa: F403
+ from FetchOpenSHA import (
+ export_to_json,
+ getERF,
+ get_rupture_distance,
+ get_source_distance,
+ get_source_rupture,
+ )
def get_rups_to_run(scenario_info, user_scenarios, num_scenarios): # noqa: C901, D103
@@ -99,7 +103,7 @@ def get_rups_to_run(scenario_info, user_scenarios, num_scenarios): # noqa: C901
elif scenario_info['Generator'].get('method', None) == 'Subsampling':
rups_to_run = list(range(num_scenarios))
else:
- sys.exit(
+ raise NotImplementedError(
f'The scenario selection method {scenario_info["Generator"].get("method", None)} is not available'
)
return rups_to_run
@@ -110,13 +114,15 @@ def load_earthquake_rupFile(scenario_info, rupFilePath): # noqa: N802, N803, D1
source_type = scenario_info['EqRupture']['Type']
try:
with open(rupFilePath) as f: # noqa: PTH123
- user_scenarios = python_json.load(f)
+ user_scenarios = json.load(f)
except: # noqa: E722
- sys.exit(f'CreateScenario: source file {rupFilePath} not found.')
+ raise FileNotFoundError( # noqa: B904
+ f'CreateScenario: source file {rupFilePath} not found.'
+ )
# number of features (i.e., ruptures)
num_scenarios = len(user_scenarios.get('features', []))
if num_scenarios < 1:
- sys.exit('CreateScenario: source file is empty.')
+ raise ValueError('CreateScenario: source file is empty.')
rups_to_run = get_rups_to_run(scenario_info, user_scenarios, num_scenarios)
# get rupture and source ids
scenario_data = {}
@@ -203,13 +209,15 @@ def load_ruptures_openquake(scenario_info, stations, work_dir, siteFile, rupFile
try:
with open(rupFile) as f: # noqa: PTH123
- user_scenarios = python_json.load(f)
+ user_scenarios = json.load(f)
except: # noqa: E722
- sys.exit(f'CreateScenario: source file {rupFile} not found.')
+ raise FileNotFoundError( # noqa: B904
+ f'CreateScenario: source file {rupFile} not found.'
+ )
# number of features (i.e., ruptures)
num_scenarios = len(user_scenarios.get('features', []))
if num_scenarios < 1:
- sys.exit('CreateScenario: source file is empty.')
+ raise ValueError('CreateScenario: source file is empty.')
rups_to_run = get_rups_to_run(scenario_info, user_scenarios, num_scenarios)
in_dir = os.path.join(work_dir, 'Input') # noqa: PTH118
oq = readinput.get_oqparam(
@@ -374,7 +382,7 @@ def load_earthquake_scenarios(scenario_info, stations, dir_info): # noqa: D103
)
try:
with open(user_scenario_file) as f: # noqa: PTH123
- user_scenarios = python_json.load(f)
+ user_scenarios = json.load(f)
except: # noqa: E722
print(f'CreateScenario: source file {user_scenario_file} not found.') # noqa: T201
return {}
@@ -589,7 +597,7 @@ def create_earthquake_scenarios( # noqa: C901, D103
if outfile is not None:
print(f'The collected point source ruptures are saved in {outfile}') # noqa: T201
with open(outfile, 'w') as f: # noqa: PTH123
- python_json.dump(pointSource_data, f, indent=2)
+ json.dump(pointSource_data, f, indent=2)
elif source_type == 'oqSourceXML':
import FetchOpenQuake
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/CreateStation.py b/modules/performRegionalEventSimulation/regionalGroundMotion/CreateStation.py
index 3b6e6499e..c505a25c3 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/CreateStation.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/CreateStation.py
@@ -37,26 +37,14 @@
# Kuanshi Zhong
#
-import socket # noqa: I001
-import sys
+import contextlib
+import socket
+import joblib
import numpy as np
import pandas as pd
-from tqdm import tqdm
-import importlib
-import subprocess
-
-
-if importlib.util.find_spec('joblib') is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'joblib']) # noqa: S603
-
-if importlib.util.find_spec('contextlib') is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'contextlib']) # noqa: S603
-
-import joblib # noqa: I001
-import contextlib
from joblib import Parallel, delayed
-import multiprocessing
+from tqdm import tqdm
@contextlib.contextmanager
@@ -78,8 +66,7 @@ def __call__(self, *args, **kwargs):
if 'stampede2' not in socket.gethostname():
- #from FetchOpenSHA import (
- from FetchOpenSHA_25_4 import (
+ from FetchOpenSHA import (
get_site_vs30_from_opensha,
get_site_z1pt0_from_opensha,
get_site_z2pt5_from_opensha,
@@ -242,7 +229,7 @@ def create_stations( # noqa: C901, PLR0912, PLR0915
]
# Check if any duplicated points
if selected_stn.duplicated(subset=[lon_label, lat_label]).any():
- sys.exit(
+ raise ValueError(
'Error: Duplicated lat and lon in the Site File (.csv), '
f'please check site \n{selected_stn[selected_stn.duplicated(subset=[lon_label, lat_label], keep = False)].index.tolist()}'
)
@@ -252,7 +239,7 @@ def create_stations( # noqa: C901, PLR0912, PLR0915
# Get Vs30
if vs30Config['Type'] == 'User-specified':
if vs30_label not in selected_stn.keys(): # noqa: SIM118
- sys.exit(
+ raise ValueError(
'ERROR: User-specified option is selected for Vs30 model but the provided.' # noqa: ISC003
+ "but the provided Site File doesn't contain a column named 'Vs30'."
+ '\nNote: the User-specified Vs30 model is only supported for Scattering Locations site definition.'
@@ -612,7 +599,7 @@ def create_stations( # noqa: C901, PLR0912, PLR0915
if stn.get('vsInferred'):
if stn.get('vsInferred') not in [0, 1]:
- sys.exit(
+ raise ValueError(
"CreateStation: Only '0' or '1' can be assigned to the" # noqa: ISC003
+ " 'vsInferred' column in the Site File (.csv), where 0 stands for false and 1 stands for true."
)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenQuake.py b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenQuake.py
index 2a7cf2574..70e480a2b 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenQuake.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenQuake.py
@@ -358,135 +358,48 @@ def openquake_config(site_info, scen_info, event_info, workDir): # noqa: C901,
with open(filename_ini, 'w') as configfile: # noqa: PTH123
cfg.write(configfile)
- # openquake module
oq_ver_loaded = None
try:
from importlib_metadata import version
except: # noqa: E722
from importlib.metadata import version
- if scen_info['EqRupture'].get('OQLocal', None):
- # using user-specific local OQ
- # first to validate the path
- if not os.path.isdir(scen_info['EqRupture'].get('OQLocal')): # noqa: PTH112
+
+ oqlocal = scen_info['EqRupture'].get('OQLocal', None)
+ if oqlocal:
+ if not os.path.isdir(oqlocal): # noqa: PTH112
print( # noqa: T201
- 'FetchOpenQuake: Local OpenQuake instance {} not found.'.format(
- scen_info['EqRupture'].get('OQLocal')
- )
+ f'FetchOpenQuake: Local OpenQuake instance {oqlocal} not found.'
)
return 0
- else: # noqa: RET505
- # getting version
- try:
- oq_ver = version('openquake.engine')
- if oq_ver:
- print( # noqa: T201
- f'FetchOpenQuake: Removing previous installation of OpenQuake {oq_ver}.'
- )
- sys.modules.pop('openquake')
- subprocess.check_call( # noqa: S603
- [
- sys.executable,
- '-m',
- 'pip',
- 'uninstall',
- '-y',
- 'openquake.engine',
- ]
- )
- except: # noqa: E722
- # no installed OQ python package
- # do nothing
- print( # noqa: T201
- 'FetchOpenQuake: No previous installation of OpenQuake python package found.'
- )
- # load the local OQ
- try:
- print('FetchOpenQuake: Setting up the user-specified local OQ.') # noqa: T201
- sys.path.insert(
- 0,
- os.path.dirname(scen_info['EqRupture'].get('OQLocal')), # noqa: PTH120
- )
- # owd = os.getcwd()
- # os.chdir(os.path.dirname(scen_info['EqRupture'].get('OQLocal')))
- if 'openquake' in list(sys.modules.keys()):
- sys.modules.pop('openquake')
- from openquake import baselib
-
- oq_ver_loaded = baselib.__version__
- # sys.modules.pop('openquake')
- # os.chdir(owd)
- except: # noqa: E722
- print( # noqa: T201
- 'FetchOpenQuake: {} cannot be loaded.'.format(
- scen_info['EqRupture'].get('OQLocal')
- )
- )
+ try:
+ print('FetchOpenQuake: Setting up the user-specified local OQ.') # noqa: T201
+ sys.path.insert(0, os.path.dirname(oqlocal)) # noqa: PTH120
+ if 'openquake' in list(sys.modules.keys()):
+ sys.modules.pop('openquake')
+ from openquake import baselib
+ oq_ver_loaded = baselib.__version__
+ except: # noqa: E722
+ print( # noqa: T201
+ f'FetchOpenQuake: {oqlocal} cannot be loaded.'
+ )
else:
- # using the official released OQ
+ required = scen_info['EqRupture'].get('OQVersion', default_oq_version)
try:
- oq_ver = version('openquake.engine')
- if oq_ver != scen_info['EqRupture'].get('OQVersion', default_oq_version):
- print( # noqa: T201
- 'FetchOpenQuake: Required OpenQuake version is not found and being installed now.'
- )
- if oq_ver:
- # pop the old version first
- sys.modules.pop('openquake')
- subprocess.check_call( # noqa: S603
- [
- sys.executable,
- '-m',
- 'pip',
- 'uninstall',
- '-y',
- 'openquake.engine',
- ]
- )
-
- # install the required version
- subprocess.check_call( # noqa: S603
- [
- sys.executable,
- '-m',
- 'pip',
- 'install',
- 'openquake.engine=='
- + scen_info['EqRupture'].get(
- 'OQVersion', default_oq_version
- ),
- '--user',
- ]
- )
- oq_ver_loaded = version('openquake.engine')
-
- else:
- oq_ver_loaded = oq_ver
-
+ oq_ver_loaded = version('openquake.engine')
except: # noqa: E722
print( # noqa: T201
- 'FetchOpenQuake: No OpenQuake is not found and being installed now.'
+ 'FetchOpenQuake: openquake.engine is not installed in this '
+ 'Python environment. Install it (e.g., '
+ '`pip install --upgrade "nheri_simcenter[r2d]"`) before '
+ 'running OpenQuake-based scenarios.'
)
- try:
- subprocess.check_call( # noqa: S603
- [
- sys.executable,
- '-m',
- 'pip',
- 'install',
- 'openquake.engine=='
- + scen_info['EqRupture'].get(
- 'OQVersion', default_oq_version
- ),
- '--user',
- ]
- )
- oq_ver_loaded = version('openquake.engine')
- except: # noqa: E722
+ else:
+ if oq_ver_loaded != required:
print( # noqa: T201
- 'FetchOpenQuake: Install of OpenQuake {} failed - please check the version.'.format(
- scen_info['EqRupture'].get('OQVersion', default_oq_version)
- )
+ f'FetchOpenQuake: this scenario expects OpenQuake '
+ f'{required}, but the installed engine is '
+ f'{oq_ver_loaded}. Proceeding with the installed version.'
)
print('FetchOpenQuake: OpenQuake configured.') # noqa: T201
@@ -891,7 +804,7 @@ def __init__( # noqa: C901
# dbserver.ensure_on()
if dbserver.get_status() == 'not-running':
if config.dbserver.multi_user:
- sys.exit(
+ raise RuntimeError(
'Please start the DbServer: ' 'see the documentation for details'
)
# otherwise start the DbServer automatically; NB: I tried to use
@@ -912,7 +825,7 @@ def __init__( # noqa: C901
waiting_seconds = 30
while dbserver.get_status() == 'not-running':
if waiting_seconds == 0:
- sys.exit(
+ raise RuntimeError(
'The DbServer cannot be started after 30 seconds. '
'Please check the configuration'
)
@@ -924,7 +837,7 @@ def __init__( # noqa: C901
# check if we are talking to the right server
err = dbserver.check_foreign()
if err:
- sys.exit(err)
+ raise RuntimeError(err)
# Copy the event_info
self.event_info = event_info
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py
index 93a530f4f..dcffbe103 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA.py
@@ -41,9 +41,6 @@
import pandas as pd
import ujson
import socket
-import subprocess
-import importlib
-import sys
import psutil
import GlobalVariable
import shapely
@@ -54,32 +51,33 @@
if GlobalVariable.JVM_started is False:
GlobalVariable.JVM_started = True
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'JPype1']) # noqa: S603
import jpype # noqa: I001, RUF100
-
- # from jpype import imports
import jpype.imports
- from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
- jpype.addClassPath('./lib/OpenSHA-1.5.2.jar')
- jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
-from java.io import * # noqa: F403
-from java.lang import * # noqa: F403
-from java.lang.reflect import * # noqa: F403
-from java.util import * # noqa: F403
-from org.opensha.commons.data import * # noqa: F403
-from org.opensha.commons.data.function import * # noqa: F403
-from org.opensha.commons.data.siteData import * # noqa: F403
-from org.opensha.commons.geo import * # noqa: F403
-from org.opensha.commons.param import * # noqa: F403
-from org.opensha.commons.param.constraint import * # noqa: F403
-from org.opensha.commons.param.event import * # noqa: F403
-from org.opensha.sha.calc import * # noqa: F403
-from org.opensha.sha.earthquake import * # noqa: F403
-from org.opensha.sha.earthquake.param import * # noqa: F403
+ jpype.addClassPath('./lib/opensha-all.jar')
+ jpype.startJVM(
+ f'-Xmx{memory_request}G',
+ convertStrings=False,
+ jvmpath=GlobalVariable.find_compatible_jvm_path(),
+ )
+
+from java.lang import Class, Double
+from java.util import ArrayList, Arrays
+from org.opensha.commons.data import Site, TimeSpan
+from org.opensha.commons.data.siteData import OrderedSiteDataProviderList, SiteData
+from org.opensha.commons.geo import Location
+from org.opensha.commons.param import Parameter
+from org.opensha.commons.param.event import ParameterChangeWarningListener
+from org.opensha.sha.earthquake import EqkRupture
+from org.opensha.sha.earthquake.param import (
+ BPTAveragingTypeOptions,
+ BackgroundRupType,
+ IncludeBackgroundOption,
+ MagDependentAperiodicityOptions,
+ ProbabilityModelOptions,
+)
from org.opensha.sha.earthquake.rupForecastImpl.Frankel02 import (
Frankel02_AdjustableEqkRupForecast,
)
@@ -90,34 +88,43 @@
from org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.MeanUCERF2 import (
MeanUCERF2,
)
-from org.opensha.sha.faultSurface import * # noqa: F403
-from org.opensha.sha.faultSurface.utils import PtSrcDistCorr
-from org.opensha.sha.imr import * # noqa: F403
-from org.opensha.sha.imr.attenRelImpl import * # noqa: F403
-from org.opensha.sha.imr.attenRelImpl.ngaw2 import * # noqa: F403
-from org.opensha.sha.imr.attenRelImpl.ngaw2.NGAW2_Wrappers import * # noqa: F403
-from org.opensha.sha.imr.param.IntensityMeasureParams import * # noqa: F403
-from org.opensha.sha.imr.param.OtherParams import * # noqa: F403
-from org.opensha.sha.util import * # noqa: F403
-from tqdm import tqdm
+from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrections
+from org.opensha.sha.gcim.imr.attenRelImpl import (
+ BommerEtAl_2009_AttenRel,
+ KS_2006_AttenRel,
+)
+from org.opensha.sha.imr.attenRelImpl import AfshariStewart_2016_AttenRel
+from org.opensha.sha.imr.attenRelImpl.ngaw2 import (
+ ASK_2014,
+ BSSA_2014,
+ CB_2014,
+ CY_2014,
+ NGAW2_Wrappers,
+)
+from org.opensha.sha.imr.param.IntensityMeasureParams import PeriodParam, SA_Param
+from org.opensha.sha.imr.param.OtherParams import StdDevTypeParam
+from org.opensha.sha.util import SiteTranslator
+# NGAW2_Wrappers is a Java class; expose its nested wrapper classes as
+# top-level names so callers can use them directly.
+ASK_2014_Wrapper = NGAW2_Wrappers.ASK_2014_Wrapper
+BSSA_2014_Wrapper = NGAW2_Wrappers.BSSA_2014_Wrapper
+CB_2014_Wrapper = NGAW2_Wrappers.CB_2014_Wrapper
+CY_2014_Wrapper = NGAW2_Wrappers.CY_2014_Wrapper
+
+# UCERF3 lives under the upstream `scratch.*` namespace.
try:
from scratch.UCERF3.erf.mean import MeanUCERF3
except ModuleNotFoundError:
- MeanUCERF3 = jpype.JClass('scratch.UCERF3.erf.mean.MeanUCERF3') # noqa: F405, RUF100
+ MeanUCERF3 = jpype.JClass('scratch.UCERF3.erf.mean.MeanUCERF3')
-from org.opensha.sha.gcim.calc import * # noqa: F403
-from org.opensha.sha.gcim.imr.attenRelImpl import * # noqa: F403
-from org.opensha.sha.gcim.imr.param.EqkRuptureParams import * # noqa: F403
-from org.opensha.sha.gcim.imr.param.IntensityMeasureParams import * # noqa: F403
+from tqdm import tqdm
def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
- # Initialization
erf = None
erf_name = scenario_info['EqRupture']['Model']
erf_selection = scenario_info['EqRupture']['ModelParameters']
- # ERF model options
if erf_name == 'WGCEP (2007) UCERF2 - Single Branch':
erf = MeanUCERF2()
if (erf_selection.get('Background Seismicity', None) == 'Exclude') and (
@@ -131,7 +138,6 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
if type(value) is int:
value = float(value) # noqa: PLW2901
erf.setParameter(key, value)
- # erf.getParameter(key).setValue(value)
elif erf_name == 'USGS/CGS 2002 Adj. Cal. ERF':
erf = Frankel02_AdjustableEqkRupForecast()
elif erf_name == 'WGCEP UCERF 1.0 (2005)':
@@ -150,22 +156,17 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
print( # noqa: T201
f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
- # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
- # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
- # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
- # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
- # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
elif erf_selection.get('preset', None) == 'FM3.1 Branch Averaged':
tmp.setPreset(MeanUCERF3.Presets.FM3_1_BRANCH_AVG)
@@ -176,24 +177,18 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
print( # noqa: T201
f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
- # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
- # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
- # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
- # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
- # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
- # Set Probability Model Option
setERFProbabilityModelOptions(tmp, erf_selection)
elif erf_selection.get('preset', None) == 'FM3.2 Branch Averaged':
tmp.setPreset(MeanUCERF3.Presets.FM3_2_BRANCH_AVG)
@@ -204,24 +199,18 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
print( # noqa: T201
f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
- # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
- # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
- # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
- # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
- # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
- # Set Probability Model Option
setERFProbabilityModelOptions(tmp, erf_selection)
else:
print( # noqa: T201
@@ -234,12 +223,29 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
else:
print('Please check the ERF model name.') # noqa: T201
- if erf_name and update_flag:
- erf.updateForecast()
- # return
+ if erf is not None:
+ _set_point_source_distance_correction(erf)
+ if update_flag:
+ erf.updateForecast()
return erf
+def _set_point_source_distance_correction(erf):
+ """Pin the ERF's point-source distance correction to NSHM_2013.
+
+ NSHM_2013 is the current USGS standard. Frankel02 and WGCEP_UCERF1 do
+ not expose this parameter, and on MeanUCERF3 it is nested under
+ "Gridded Seismicity Settings"; the try/except skips those cases.
+ """
+ try:
+ erf.setParameter(
+ 'Point Source Distance Correction',
+ PointSourceDistanceCorrections.NSHM_2013,
+ )
+ except: # noqa: E722
+ pass
+
+
def setERFbackgroundOptions(erf, selection): # noqa: N802, D103
option = selection.get('Background Seismicity', None)
if option == 'Include':
@@ -375,8 +381,6 @@ def get_source_rupture(erf, source_index, rupture_index): # noqa: D103
def get_source_distance(erf, source_index, lat, lon): # noqa: D103
rupSource = erf.getSource(source_index) # noqa: N806
sourceSurface = rupSource.getSourceSurface() # noqa: N806
- # print(lon)
- # print(lat)
distToSource = [] # noqa: N806
for i in range(len(lat)):
distToSource.append( # noqa: PERF401
@@ -414,7 +418,7 @@ def get_rupture_info_ASK2014_aftershock(
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
+ return mainshock.geometry.iloc[0].distance(centroid) / 1000
def get_rupture_info_CY2014(erf, source_index, rupture_index, siteList): # noqa: N802, N803, D103
@@ -477,7 +481,6 @@ def horzDistanceFast(lat1, lon1, lat2, lon2): # noqa: N802, D103
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
c = 2 * np.arcsin(np.sqrt(a))
EARTH_RADIUS_MEAN = 6371.0072 # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/commons/geo/GeoTools.java#L22 # noqa: N806
- # return EARTH_RADIUS_MEAN * np.sqrt((dLat * dLat) + (dLon * dLon))
return EARTH_RADIUS_MEAN * c
@@ -492,7 +495,6 @@ def getPtSrcDistCorr(horzDist, mag, type): # noqa: A002, N802, N803, D103
print( # noqa: T201
'The NSHMP08 rJB correction has not been implemented. corr=1.0 is used instead'
)
- # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/utils/PtSrcDistCorr.java#L20
return 1.0
else:
return 1.0
@@ -517,6 +519,7 @@ def get_PointSource_info_CY2014(source_info, siteList): # noqa: N802, N803, D10
'dip': float(source_info['AverageDip']),
'width': 0.0, # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/PointSurface.java#L68
'zTop': sourceDepth,
+ 'zHyp': sourceDepth,
'aveRake': float(source_info['AverageRake']),
}
return site_rup_info, siteList
@@ -532,28 +535,22 @@ def export_to_json( # noqa: C901, D103
maxDistance=1000.0, # noqa: N803
use_hdf5=False, # noqa: FBT002
):
- # Initializing
erf_data = {'type': 'FeatureCollection'}
site_loc = Location(site_loc[0], site_loc[1]) # type: ignore # noqa: F405
site = Site(site_loc) # type: ignore # noqa: F405
- # Total source number
num_sources = erf.getNumSources()
source_tag = []
source_dist = []
for i in tqdm(range(num_sources), desc=f'Find sources with in {maxDistance} km'):
rup_source = erf.getSource(i)
distance_to_source = rup_source.getMinDistance(site)
- # sourceSurface = rupSource.getSourceSurface()
- # distanceToSource = sourceSurface.getDistanceRup(site_loc)
source_tag.append(i)
source_dist.append(distance_to_source)
df = pd.DataFrame.from_dict({'sourceID': source_tag, 'sourceDist': source_dist}) # noqa: PD901
- # Sorting sources
source_collection = df.sort_values(['sourceDist'], ascending=(True))
source_collection = source_collection[
source_collection['sourceDist'] < maxDistance
]
- # Collecting source features
if not use_hdf5:
feature_collection = []
for i in tqdm(
@@ -561,7 +558,6 @@ def export_to_json( # noqa: C901, D103
desc=f'Find ruptures with in {maxDistance} km',
):
source_index = source_collection.iloc[i, 0]
- # Getting rupture distances
rupSource = erf.getSource(source_index) # noqa: N806
try:
rupList = rupSource.getRuptureList() # noqa: N806
@@ -575,38 +571,26 @@ def export_to_json( # noqa: C901, D103
rup_dist = []
for j in range(rupList.size()):
ruptureSurface = rupList.get(j).getRuptureSurface() # noqa: N806
- # If pointsource rupture distance correction
- if isinstance(ruptureSurface, PointSurface): # noqa: F405
- # or 'FIELD' or 'NSHMP08'
- distCorrType = PtSrcDistCorr.Type.NONE # noqa: N806
- (PointSurface @ ruptureSurface).setDistCorrMagAndType( # noqa: F405
- rupList.get(j).getMag(), distCorrType
- )
cur_dist = ruptureSurface.getDistanceRup(site_loc)
rup_tag.append(j)
if cur_dist < maxDistance:
rup_dist.append(cur_dist)
else:
- # exceeding the maxDistance requirement
rup_dist.append(-1.0)
df = pd.DataFrame.from_dict({'rupID': rup_tag, 'rupDist': rup_dist}) # noqa: PD901
- # Sorting
rup_collection = df.sort_values(['rupDist'], ascending=(True))
- # Preparing the dict of ruptures
for j in range(rupList.size()):
cur_dict = dict() # noqa: C408
cur_dict.update({'type': 'Feature'})
rup_index = rup_collection.iloc[j, 0]
cur_dist = rup_collection.iloc[j, 1]
if cur_dist <= 0.0:
- # skipping ruptures with distance exceeding the maxDistance
continue
rupture = rupList.get(rup_index)
maf = rupture.getMeanAnnualRate(erf.getTimeSpan().getDuration())
if maf <= 0.0:
continue
ruptureSurface = rupture.getRuptureSurface() # noqa: N806
- # Properties
cur_dict['properties'] = dict() # noqa: C408
name = str(rupSource.getName())
if EqName is not None:
@@ -645,10 +629,8 @@ def export_to_json( # noqa: C901, D103
cur_dict['properties'].update(
{'MeanAnnualRate': abs(float(maf))}
)
- # Geometry
cur_dict['geometry'] = dict() # noqa: C408
if ruptureSurface.isPointSurface():
- # Point source
pointSurface = ruptureSurface # noqa: N806
location = pointSurface.getLocation()
cur_dict['geometry'].update({'type': 'Point'})
@@ -661,7 +643,6 @@ def export_to_json( # noqa: C901, D103
}
)
else:
- # Line source
try:
trace = ruptureSurface.getUpperEdge()
except: # noqa: E722
@@ -673,9 +654,7 @@ def export_to_json( # noqa: C901, D103
)
cur_dict['geometry'].update({'type': 'LineString'})
cur_dict['geometry'].update({'coordinates': coordinates})
- # Appending
feature_collection.append(cur_dict)
- # sort the list
maf_list_n = [-x['properties']['MeanAnnualRate'] for x in feature_collection]
sort_ids = np.argsort(maf_list_n)
feature_collection_sorted = [feature_collection[i] for i in sort_ids]
@@ -694,17 +673,14 @@ def export_to_json( # noqa: C901, D103
import h5py
with h5py.File(outfile, 'w') as h5file:
- # Store the geometry as a string array
h5file.create_dataset(
'geometry',
data=gdf.geometry.astype(str).values.astype('S'), # noqa: PD011, F405
) # noqa: F405, PD011, RUF100
- # return
return erf_data
def CreateIMRInstance(gmpe_name): # noqa: N802, D103
- # GMPE name map
gmpe_map = {
str(ASK_2014.NAME): ASK_2014_Wrapper.class_.getName(), # noqa: F405
str(BSSA_2014.NAME): BSSA_2014_Wrapper.class_.getName(), # noqa: F405
@@ -718,17 +694,28 @@ def CreateIMRInstance(gmpe_name): # noqa: N802, D103
AfshariStewart_2016_AttenRel.NAME # noqa: F405
): AfshariStewart_2016_AttenRel.class_.getName(), # noqa: F405
}
- # Mapping GMPE name
imrClassName = gmpe_map.get(gmpe_name) # noqa: N806
if imrClassName is None:
return imrClassName
- # Getting the java class
imrClass = Class.forName(imrClassName) # noqa: N806, F405
- ctor = imrClass.getConstructor()
- imr = ctor.newInstance()
- # Setting default parameters
+ # Some GMMs (KS_2006, BommerEtAl_2009, ...) take a
+ # ParameterChangeWarningListener instead of a no-arg constructor.
+ try:
+ ctor = imrClass.getConstructor()
+ imr = ctor.newInstance()
+ except: # noqa: E722
+ try:
+ @jpype.JImplements(ParameterChangeWarningListener)
+ class _NoopListener:
+ @jpype.JOverride
+ def parameterChangeWarning(self, event): # noqa: ARG002
+ pass
+
+ ctor = imrClass.getConstructor(ParameterChangeWarningListener)
+ imr = ctor.newInstance(_NoopListener())
+ except: # noqa: E722
+ return None
imr.setParamDefaults()
- # return
return imr
@@ -755,13 +742,11 @@ def get_DataSource(paramName, siteData): # noqa: N802, N803, D103
def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
- # GMPE
try:
imr = CreateIMRInstance(gmpe_name)
except: # noqa: E722
print('Please check GMPE name.') # noqa: T201
return 1
- # Site data
sites = ArrayList() # noqa: F405
for cur_site in siteSpec:
cur_loc = Location( # noqa: F405
@@ -776,15 +761,11 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
print( # noqa: T201
'remote getAllAvailableData is not available temporarily, will use site Vs30 in the site csv file.'
)
- # return 1
siteTrans = SiteTranslator() # noqa: N806, F405
- # Looping over all sites
site_prop = []
for i in range(len(siteSpec)):
site_tmp = dict() # noqa: C408
- # Current site
site = sites.get(i)
- # Location
cur_site = siteSpec[i]
locResults = { # noqa: N806
'Latitude': cur_site['Location']['Latitude'],
@@ -798,7 +779,6 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
siteDataValues.add(availableSiteData.get(j).getValue(i))
imrSiteParams = imr.getSiteParams() # noqa: N806
siteDataResults = [] # noqa: N806
- # Setting site parameters
for j in range(imrSiteParams.size()):
siteParam = imrSiteParams.getByIndex(j) # noqa: N806
newParam = Parameter.clone(siteParam) # noqa: N806, F405
@@ -856,13 +836,10 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
}
)
site.addParameter(newParam)
- # End for j
- # Updating site specifications
siteSpec[i] = cur_site
site_tmp.update({'Location': locResults, 'SiteData': siteDataResults})
site_prop.append(site_tmp)
- # Return
return siteSpec, sites, site_prop
@@ -876,20 +853,16 @@ def get_IM( # noqa: C901, N802, D103
station_info,
im_info,
):
- # GMPE name
gmpe_name = gmpe_info['Type']
- # Creating intensity measure relationship instance
try:
imr = CreateIMRInstance(gmpe_name)
except: # noqa: E722
print('Please check GMPE name.') # noqa: T201
return 1, station_info
- # Getting supported intensity measure types
ims = imr.getSupportedIntensityMeasures()
saParam = ims.getParameter(SA_Param.NAME) # noqa: N806, F405
supportedPeriods = saParam.getPeriodParam().getPeriods() # noqa: N806
Arrays.sort(supportedPeriods) # noqa: F405
- # Rupture
eqRup = EqkRupture() # noqa: N806, F405
if source_info['Type'] == 'PointSource':
eqRup.setMag(source_info['Magnitude'])
@@ -905,38 +878,27 @@ def get_IM( # noqa: C901, N802, D103
elif source_info['Type'] == 'ERF':
timeSpan = TimeSpan(TimeSpan.NONE, TimeSpan.YEARS) # noqa: N806, F405
erfParams = source_info.get('Parameters', None) # noqa: N806
- # Additional parameters (if any)
if erfParams is not None:
for k in erfParams.keys:
erf.setParameter(k, erfParams[k])
- # Time span
timeSpan = erf.getTimeSpan() # noqa: N806
- # Source
eqSource = erf.getSource(source_info['SourceIndex']) # noqa: N806
eqSource.getName()
- # Rupture
eqRup = eqSource.getRupture(source_info['RuptureIndex']) # noqa: N806
- # Properties
magnitude = eqRup.getMag()
averageDip = eqRup.getRuptureSurface().getAveDip() # noqa: N806, F841
averageRake = eqRup.getAveRake() # noqa: N806, F841
- # Probability
probEqRup = eqRup # noqa: N806
probability = probEqRup.getProbability() # noqa: F841
- # MAF
meanAnnualRate = probEqRup.getMeanAnnualRate(timeSpan.getDuration()) # noqa: N806
- # Rupture surface
surface = eqRup.getRuptureSurface() # noqa: F841
- # Setting up imr
imr.setEqkRupture(eqRup)
imrParams = gmpe_info['Parameters'] # noqa: N806
if bool(imrParams):
for k in imrParams.keys(): # noqa: SIM118
imr.getParameter(k).setValue(imrParams[k])
- # Station
if station_info['Type'] == 'SiteList':
siteSpec = station_info['SiteList'] # noqa: N806
- # Intensity measure
periods = im_info.get('Periods', None)
if periods is not None:
periods = supportedPeriods
@@ -951,15 +913,11 @@ def get_IM( # noqa: C901, N802, D103
tag_PGA = True # noqa: N806
if 'PGV' in im_info['Type']:
tag_PGV = True # noqa: N806
- # Looping over sites
gm_collector = []
for i in range(len(siteSpec)):
gmResults = site_prop[i] # noqa: N806
- # Current site
site = sites.get(i)
- # Location
cur_site = siteSpec[i] # noqa: F841
- # Set up the site in the imr
imr.setSite(site)
try:
stdDevParam = imr.getParameter(StdDevTypeParam.NAME) # noqa: N806, F405
@@ -967,7 +925,7 @@ def get_IM( # noqa: C901, N802, D103
StdDevTypeParam.STD_DEV_TYPE_INTER # noqa: F405
) and stdDevParam.isAllowed(StdDevTypeParam.STD_DEV_TYPE_INTRA) # noqa: F405
except: # noqa: E722
- stdDevParaam = None # noqa: N806, F841
+ stdDevParam = None # noqa: N806, F841
hasIEStats = False # noqa: N806
cur_T = im_info.get('Periods', None) # noqa: N806
if tag_SA:
@@ -996,7 +954,6 @@ def get_IM( # noqa: C901, N802, D103
saResult['IntraEvStdDev'].append(float(intraEvStdDev))
gmResults.update({'lnSA': saResult})
if tag_PGA:
- # for PGV current T = 0
cur_T = [0.00] # noqa: N806
pgaResult = {'Mean': [], 'TotalStdDev': []} # noqa: N806
if hasIEStats:
@@ -1016,7 +973,6 @@ def get_IM( # noqa: C901, N802, D103
pgaResult['IntraEvStdDev'].append(float(intraEvStdDev))
gmResults.update({'lnPGA': pgaResult})
if tag_PGV:
- # for PGV current T = 0
cur_T = [0.00] # noqa: N806
pgvResult = {'Mean': [], 'TotalStdDev': []} # noqa: N806
if hasIEStats:
@@ -1037,10 +993,8 @@ def get_IM( # noqa: C901, N802, D103
gmResults.update({'lnPGV': pgvResult})
gm_collector.append(gmResults)
- # Updating station information
if station_info['Type'] == 'SiteList':
station_info.update({'SiteList': siteSpec})
- # Final results
res = {
'Magnitude': magnitude,
'MeanAnnualRate': meanAnnualRate,
@@ -1049,22 +1003,18 @@ def get_IM( # noqa: C901, N802, D103
'Periods': cur_T,
'GroundMotions': gm_collector,
}
- # return
return res, station_info
def get_site_vs30_from_opensha(lat, lon, vs30model='CGS/Wills VS30 Map (2015)'): # noqa: D103
- # set up site java object
sites = ArrayList() # noqa: F405
num_sites = len(lat)
for i in range(num_sites):
sites.add(Site(Location(lat[i], lon[i]))) # noqa: F405
- # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
- # search name
vs30 = []
for i in range(int(siteData.size())):
cur_siteData = siteData.get(i) # noqa: N806
@@ -1076,21 +1026,33 @@ def get_site_vs30_from_opensha(lat, lon, vs30model='CGS/Wills VS30 Map (2015)'):
else: # noqa: RET508
continue
- # check if any nan (Wills Map return nan for offshore sites)
- # Using global vs30 as default patch - 'Global Vs30 from Topographic Slope (Wald & Allen 2008)'
+ # The Wills 2015 map returns NaN for offshore sites; fall back to the
+ # global topographic-slope Vs30 model. Look up the fallback provider by
+ # name so an upstream reorder of the provider list does not silently
+ # route this code to the wrong dataset.
if any([np.isnan(x) for x in vs30]): # noqa: C419
- non_list = np.where(np.isnan(vs30))[0].tolist()
- for i in non_list:
- vs30[i] = float(siteData.get(3).getValue(i).getValue())
+ fallback_name = 'Global Vs30 from Topographic Slope (Wald & Allen 2008)'
+ fallback_provider_data = None
+ for i in range(int(siteData.size())):
+ if str(siteData.get(i).getSourceName()) == fallback_name:
+ fallback_provider_data = siteData.get(i)
+ break
+ if fallback_provider_data is None:
+ print( # noqa: T201
+ f'FetchOpenSHA: fallback Vs30 provider "{fallback_name}" '
+ f'not found; leaving NaN Vs30 values unchanged.'
+ )
+ else:
+ non_list = np.where(np.isnan(vs30))[0].tolist()
+ for i in non_list:
+ vs30[i] = float(fallback_provider_data.getValue(i).getValue())
- # return
return vs30
def get_site_z1pt0_from_opensha(lat, lon): # noqa: D103
sites = ArrayList() # noqa: F405
sites.add(Site(Location(lat, lon))) # noqa: F405
- # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
for data in siteData:
@@ -1104,7 +1066,6 @@ def get_site_z1pt0_from_opensha(lat, lon): # noqa: D103
def get_site_z2pt5_from_opensha(lat, lon): # noqa: D103
sites = ArrayList() # noqa: F405
sites.add(Site(Location(lat, lon))) # noqa: F405
- # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
for data in siteData:
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_25_4.py b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_1.5.2.py
similarity index 82%
rename from modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_25_4.py
rename to modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_1.5.2.py
index 813784da5..453e9296c 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_25_4.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/FetchOpenSHA_1.5.2.py
@@ -1,3 +1,46 @@
+# # noqa: INP001, D100
+# Copyright (c) 2018 Leland Stanford Junior University
+# Copyright (c) 2018 The Regents of the University of California
+#
+# This file is part of the SimCenter Backend Applications
+#
+# Redistribution and use in source and binary forms, with or without
+# modification, are permitted provided that the following conditions are met:
+#
+# 1. Redistributions of source code must retain the above copyright notice,
+# this list of conditions and the following disclaimer.
+#
+# 2. Redistributions in binary form must reproduce the above copyright notice,
+# this list of conditions and the following disclaimer in the documentation
+# and/or other materials provided with the distribution.
+#
+# 3. Neither the name of the copyright holder nor the names of its contributors
+# may be used to endorse or promote products derived from this software without
+# specific prior written permission.
+#
+# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
+# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
+# ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
+# LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
+# CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
+# SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
+# INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
+# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
+# ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
+# POSSIBILITY OF SUCH DAMAGE.
+#
+# You should have received a copy of the BSD 3-Clause License along with
+# this file. If not, see .
+#
+# Contributors:
+# Kuanshi Zhong
+#
+
+# Legacy module targeting OpenSHA v1.5.2. The current module is
+# `FetchOpenSHA.py`. To use this file, place a compatible OpenSHA-1.5.2.jar
+# at `lib/OpenSHA-1.5.2.jar` (the jar is not shipped with the backend).
+
import numpy as np # noqa: I001
import pandas as pd
import ujson
@@ -19,14 +62,14 @@
subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'JPype1']) # noqa: S603
import jpype # noqa: I001, RUF100
+ # from jpype import imports
import jpype.imports
from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
- jpype.addClassPath('./lib/opensha-all.jar')
+ jpype.addClassPath('./lib/OpenSHA-1.5.2.jar')
jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
-
from java.io import * # noqa: F403
from java.lang import * # noqa: F403
from java.lang.reflect import * # noqa: F403
@@ -34,8 +77,7 @@
from org.opensha.commons.data import * # noqa: F403
from org.opensha.commons.data.function import * # noqa: F403
from org.opensha.commons.data.siteData import * # noqa: F403
-#from org.opensha.commons.geo import * # noqa: F403
-from org.opensha.commons.geo import Location # noqa: F403
+from org.opensha.commons.geo import * # noqa: F403
from org.opensha.commons.param import * # noqa: F403
from org.opensha.commons.param.constraint import * # noqa: F403
from org.opensha.commons.param.event import * # noqa: F403
@@ -53,22 +95,7 @@
MeanUCERF2,
)
from org.opensha.sha.faultSurface import * # noqa: F403
-
-try:
- from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrection
- from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrections
-except ImportError:
- try:
- import jpype
- PointSourceDistanceCorrection = jpype.JClass('org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrection')
- PointSourceDistanceCorrections = jpype.JClass('org.opensha.sha.faultSurface.utils.PointSourceDistanceCorrections')
- except:
- PointSourceDistanceCorrection = None
- PointSourceDistanceCorrections = None
-
-#from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrection
-#from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrections
-
+from org.opensha.sha.faultSurface.utils import PtSrcDistCorr
from org.opensha.sha.imr import * # noqa: F403
from org.opensha.sha.imr.attenRelImpl import * # noqa: F403
from org.opensha.sha.imr.attenRelImpl.ngaw2 import * # noqa: F403
@@ -90,22 +117,25 @@
def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
+ # Initialization
erf = None
erf_name = scenario_info['EqRupture']['Model']
erf_selection = scenario_info['EqRupture']['ModelParameters']
+ # ERF model options
if erf_name == 'WGCEP (2007) UCERF2 - Single Branch':
erf = MeanUCERF2()
if (erf_selection.get('Background Seismicity', None) == 'Exclude') and (
- 'Treat Background Seismicity As' in erf_selection.keys()
+ 'Treat Background Seismicity As' in erf_selection.keys() # noqa: SIM118
):
value = erf_selection.pop('Treat Background Seismicity As')
- print(
- f'Background Seismicity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
+ print( # noqa: T201
+ f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
for key, value in erf_selection.items():
if type(value) is int:
- value = float(value)
+ value = float(value) # noqa: PLW2901
erf.setParameter(key, value)
+ # erf.getParameter(key).setValue(value)
elif erf_name == 'USGS/CGS 2002 Adj. Cal. ERF':
erf = Frankel02_AdjustableEqkRupForecast()
elif erf_name == 'WGCEP UCERF 1.0 (2005)':
@@ -118,70 +148,87 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
):
tmp.setPreset(MeanUCERF3.Presets.BOTH_FM_BRANCH_AVG)
if (erf_selection.get('Background Seismicity', None) == 'Exclude') and (
- 'Treat Background Seismicity As' in erf_selection.keys()
+ 'Treat Background Seismicity As' in erf_selection.keys() # noqa: SIM118
):
value = erf_selection.pop('Treat Background Seismicity As')
- print(
- f'Background Seismicity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
+ print( # noqa: T201
+ f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
+ # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
+ # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
+ # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
+ # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
+ # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
elif erf_selection.get('preset', None) == 'FM3.1 Branch Averaged':
tmp.setPreset(MeanUCERF3.Presets.FM3_1_BRANCH_AVG)
if (erf_selection.get('Background Seismicity', None) == 'Exclude') and (
- 'Treat Background Seismicity As' in erf_selection.keys()
+ 'Treat Background Seismicity As' in erf_selection.keys() # noqa: SIM118
):
value = erf_selection.pop('Treat Background Seismicity As')
- print(
- f'Background Seismicity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
+ print( # noqa: T201
+ f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
+ # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
+ # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
+ # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
+ # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
+ # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
+ # Set Probability Model Option
setERFProbabilityModelOptions(tmp, erf_selection)
elif erf_selection.get('preset', None) == 'FM3.2 Branch Averaged':
tmp.setPreset(MeanUCERF3.Presets.FM3_2_BRANCH_AVG)
if (erf_selection.get('Background Seismicity', None) == 'Exclude') and (
- 'Treat Background Seismicity As' in erf_selection.keys()
+ 'Treat Background Seismicity As' in erf_selection.keys() # noqa: SIM118
):
value = erf_selection.pop('Treat Background Seismicity As')
- print(
- f'Background Seismicity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
+ print( # noqa: T201
+ f'Background Seismicvity is set as Excluded, Treat Background Seismicity As: {value} is ignored'
)
+ # Some parameters in MeanUCERF3 have overloaded setValue() Need to set one by one
+ # Set Apply Aftershock Filter
if erf_selection.get('Apply Aftershock Filter', None):
tmp.setParameter(
'Apply Aftershock Filter',
erf_selection['Apply Aftershock Filter'],
)
+ # Set Aleatoiry mag-area stdDev
if erf_selection.get('Aleatory Mag-Area StdDev', None):
tmp.setParameter(
'Aleatory Mag-Area StdDev',
erf_selection['Aleatory Mag-Area StdDev'],
)
+ # Set IncludeBackgroundOpetion
setERFbackgroundOptions(tmp, erf_selection)
+ # Set Treat Background Seismicity As Option
setERFtreatBackgroundOptions(tmp, erf_selection)
+ # Set Probability Model Option
setERFProbabilityModelOptions(tmp, erf_selection)
else:
- print(
+ print( # noqa: T201
f"""The specified Mean UCERF3 preset {erf_selection.get("preset", None)} is not implemented"""
)
erf = tmp
@@ -189,10 +236,11 @@ def getERF(scenario_info, update_flag=True): # noqa: FBT002, C901, N802, D103
elif erf_name == 'WGCEP Eqk Rate Model 2 ERF':
erf = UCERF2()
else:
- print('Please check the ERF model name.')
+ print('Please check the ERF model name.') # noqa: T201
if erf_name and update_flag:
erf.updateForecast()
+ # return
return erf
@@ -331,11 +379,14 @@ def get_source_rupture(erf, source_index, rupture_index): # noqa: D103
def get_source_distance(erf, source_index, lat, lon): # noqa: D103
rupSource = erf.getSource(source_index) # noqa: N806
sourceSurface = rupSource.getSourceSurface() # noqa: N806
+ # print(lon)
+ # print(lat)
distToSource = [] # noqa: N806
for i in range(len(lat)):
- distToSource.append(
+ distToSource.append( # noqa: PERF401
float(sourceSurface.getDistanceRup(Location(lat[i], lon[i]))) # noqa: F405
)
+
return distToSource
@@ -344,11 +395,11 @@ def get_rupture_distance(erf, source_index, rupture_index, lat, lon): # noqa: D
rupSurface = rupSource.getRupture(rupture_index).getRuptureSurface() # noqa: N806
distToRupture = [] # noqa: N806
for i in range(len(lat)):
- distToRupture.append(
+ distToRupture.append( # noqa: PERF401
float(rupSurface.getDistanceRup(Location(lat[i], lon[i]))) # noqa: F405
)
- return distToRupture
+ return distToRupture
def get_rupture_info_ASK2014_aftershock(
erf, source_index, rupture_indx, mainshock):
@@ -367,7 +418,7 @@ def get_rupture_info_ASK2014_aftershock(
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
+ 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
@@ -375,6 +426,7 @@ def get_rupture_info_CY2014(erf, source_index, rupture_index, siteList): # noqa
rupList = rupSource.getRuptureList() # noqa: N806
rupSurface = rupList.get(rupture_index).getRuptureSurface() # noqa: N806
if rupList.get(rupture_index).getHypocenterLocation() is None:
+ # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/nshmp2/imr/ngaw2/NSHMP14_WUS_CB.java#L242
dip = float(rupSurface.getAveDip())
width = float(rupSurface.getAveWidth())
zTop = float(rupSurface.getAveRupTopDepth()) # noqa: N806
@@ -428,26 +480,30 @@ def horzDistanceFast(lat1, lon1, lat2, lon2): # noqa: N802, D103
dlat = np.abs(lat2 - lat1)
a = np.sin(dlat / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2) ** 2
c = 2 * np.arcsin(np.sqrt(a))
- EARTH_RADIUS_MEAN = 6371.0072 # noqa: N806
+ EARTH_RADIUS_MEAN = 6371.0072 # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/commons/geo/GeoTools.java#L22 # noqa: N806
+ # return EARTH_RADIUS_MEAN * np.sqrt((dLat * dLat) + (dLon * dLon))
return EARTH_RADIUS_MEAN * c
def getPtSrcDistCorr(horzDist, mag, type): # noqa: A002, N802, N803, D103
+ # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/utils/PtSrcDistCorr.java#L20
if type == 'FIELD':
rupLen = np.power(10.0, -3.22 + 0.69 * mag) # noqa: N806
return 0.7071 + (1.0 - 0.7071) / (
1 + np.power(rupLen / (horzDist * 0.87), 1.1)
)
- elif type == 'NSHMP08':
- print(
+ elif type == 'NSHMP08': # noqa: RET505
+ print( # noqa: T201
'The NSHMP08 rJB correction has not been implemented. corr=1.0 is used instead'
)
+ # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/utils/PtSrcDistCorr.java#L20
return 1.0
else:
return 1.0
def get_PointSource_info_CY2014(source_info, siteList): # noqa: N802, N803, D103
+ # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/PointSurface.java#L118
sourceLat = source_info['Location']['Latitude'] # noqa: N806
sourceLon = source_info['Location']['Longitude'] # noqa: N806
sourceDepth = source_info['Location']['Depth'] # noqa: N806
@@ -463,9 +519,8 @@ def get_PointSource_info_CY2014(source_info, siteList): # noqa: N802, N803, D10
siteList[i].update({'rX': rX})
site_rup_info = {
'dip': float(source_info['AverageDip']),
- 'width': 0.0,
+ 'width': 0.0, # https://github.com/opensha/opensha/blob/master/src/main/java/org/opensha/sha/faultSurface/PointSurface.java#L68
'zTop': sourceDepth,
- 'zHyp': sourceDepth,
'aveRake': float(source_info['AverageRake']),
}
return site_rup_info, siteList
@@ -481,22 +536,28 @@ def export_to_json( # noqa: C901, D103
maxDistance=1000.0, # noqa: N803
use_hdf5=False, # noqa: FBT002
):
+ # Initializing
erf_data = {'type': 'FeatureCollection'}
- site_loc = Location(site_loc[0], site_loc[1]) # noqa: F405
- site = Site(site_loc) # noqa: F405
+ site_loc = Location(site_loc[0], site_loc[1]) # type: ignore # noqa: F405
+ site = Site(site_loc) # type: ignore # noqa: F405
+ # Total source number
num_sources = erf.getNumSources()
source_tag = []
source_dist = []
for i in tqdm(range(num_sources), desc=f'Find sources with in {maxDistance} km'):
rup_source = erf.getSource(i)
distance_to_source = rup_source.getMinDistance(site)
+ # sourceSurface = rupSource.getSourceSurface()
+ # distanceToSource = sourceSurface.getDistanceRup(site_loc)
source_tag.append(i)
source_dist.append(distance_to_source)
- df = pd.DataFrame.from_dict({'sourceID': source_tag, 'sourceDist': source_dist})
+ df = pd.DataFrame.from_dict({'sourceID': source_tag, 'sourceDist': source_dist}) # noqa: PD901
+ # Sorting sources
source_collection = df.sort_values(['sourceDist'], ascending=(True))
source_collection = source_collection[
source_collection['sourceDist'] < maxDistance
]
+ # Collecting source features
if not use_hdf5:
feature_collection = []
for i in tqdm(
@@ -504,6 +565,7 @@ def export_to_json( # noqa: C901, D103
desc=f'Find ruptures with in {maxDistance} km',
):
source_index = source_collection.iloc[i, 0]
+ # Getting rupture distances
rupSource = erf.getSource(source_index) # noqa: N806
try:
rupList = rupSource.getRuptureList() # noqa: N806
@@ -517,29 +579,38 @@ def export_to_json( # noqa: C901, D103
rup_dist = []
for j in range(rupList.size()):
ruptureSurface = rupList.get(j).getRuptureSurface() # noqa: N806
- # v25.4: PtSrcDistCorr removed. Original code used NONE
- # (no correction). In v25.4, PointSurface applies no
- # correction by default, so no action needed.
+ # If pointsource rupture distance correction
+ if isinstance(ruptureSurface, PointSurface): # noqa: F405
+ # or 'FIELD' or 'NSHMP08'
+ distCorrType = PtSrcDistCorr.Type.NONE # noqa: N806
+ (PointSurface @ ruptureSurface).setDistCorrMagAndType( # noqa: F405
+ rupList.get(j).getMag(), distCorrType
+ )
cur_dist = ruptureSurface.getDistanceRup(site_loc)
rup_tag.append(j)
if cur_dist < maxDistance:
rup_dist.append(cur_dist)
else:
+ # exceeding the maxDistance requirement
rup_dist.append(-1.0)
- df = pd.DataFrame.from_dict({'rupID': rup_tag, 'rupDist': rup_dist})
+ df = pd.DataFrame.from_dict({'rupID': rup_tag, 'rupDist': rup_dist}) # noqa: PD901
+ # Sorting
rup_collection = df.sort_values(['rupDist'], ascending=(True))
+ # Preparing the dict of ruptures
for j in range(rupList.size()):
cur_dict = dict() # noqa: C408
cur_dict.update({'type': 'Feature'})
rup_index = rup_collection.iloc[j, 0]
cur_dist = rup_collection.iloc[j, 1]
if cur_dist <= 0.0:
+ # skipping ruptures with distance exceeding the maxDistance
continue
rupture = rupList.get(rup_index)
maf = rupture.getMeanAnnualRate(erf.getTimeSpan().getDuration())
if maf <= 0.0:
continue
ruptureSurface = rupture.getRuptureSurface() # noqa: N806
+ # Properties
cur_dict['properties'] = dict() # noqa: C408
name = str(rupSource.getName())
if EqName is not None:
@@ -553,16 +624,18 @@ def export_to_json( # noqa: C901, D103
cur_dict['properties'].update({'Rupture': int(rup_index)})
cur_dict['properties'].update({'Source': int(source_index)})
if outfile is not None:
+ # these calls are time-consuming, so only run them if one needs
+ # detailed outputs of the sources
cur_dict['properties'].update({'Distance': float(cur_dist)})
distanceRup = rupture.getRuptureSurface().getDistanceRup( # noqa: N806
site_loc
- )
+ ) # noqa: N806, RUF100
cur_dict['properties'].update(
{'DistanceRup': float(distanceRup)}
)
distanceSeis = rupture.getRuptureSurface().getDistanceSeis( # noqa: N806
site_loc
- )
+ ) # noqa: N806, RUF100
cur_dict['properties'].update(
{'DistanceSeis': float(distanceSeis)}
)
@@ -576,8 +649,10 @@ def export_to_json( # noqa: C901, D103
cur_dict['properties'].update(
{'MeanAnnualRate': abs(float(maf))}
)
+ # Geometry
cur_dict['geometry'] = dict() # noqa: C408
if ruptureSurface.isPointSurface():
+ # Point source
pointSurface = ruptureSurface # noqa: N806
location = pointSurface.getLocation()
cur_dict['geometry'].update({'type': 'Point'})
@@ -590,45 +665,50 @@ def export_to_json( # noqa: C901, D103
}
)
else:
+ # Line source
try:
trace = ruptureSurface.getUpperEdge()
except: # noqa: E722
trace = ruptureSurface.getEvenlyDiscritizedUpperEdge()
coordinates = []
for k in trace:
- coordinates.append(
+ coordinates.append( # noqa: PERF401
[float(k.getLongitude()), float(k.getLatitude())]
)
cur_dict['geometry'].update({'type': 'LineString'})
cur_dict['geometry'].update({'coordinates': coordinates})
+ # Appending
feature_collection.append(cur_dict)
+ # sort the list
maf_list_n = [-x['properties']['MeanAnnualRate'] for x in feature_collection]
sort_ids = np.argsort(maf_list_n)
feature_collection_sorted = [feature_collection[i] for i in sort_ids]
del feature_collection
erf_data.update({'features': feature_collection_sorted})
- print(
+ print( # noqa: T201
f'FetchOpenSHA: total {len(feature_collection_sorted)} ruptures are collected.'
)
if outfile is not None:
- print(
+ print( # noqa: T201
f'The collected ruptures are sorted by MeanAnnualRate and saved in {outfile}'
)
- print(f"FetchOpenSRA_24 {outfile}")
- with open(outfile, 'w') as f:
+ with open(outfile, 'w') as f: # noqa: PTH123
ujson.dump(erf_data, f, indent=2)
else:
import h5py
with h5py.File(outfile, 'w') as h5file:
+ # Store the geometry as a string array
h5file.create_dataset(
'geometry',
- data=gdf.geometry.astype(str).values.astype('S'),
- )
+ data=gdf.geometry.astype(str).values.astype('S'), # noqa: PD011, F405
+ ) # noqa: F405, PD011, RUF100
+ # return
return erf_data
def CreateIMRInstance(gmpe_name): # noqa: N802, D103
+ # GMPE name map
gmpe_map = {
str(ASK_2014.NAME): ASK_2014_Wrapper.class_.getName(), # noqa: F405
str(BSSA_2014.NAME): BSSA_2014_Wrapper.class_.getName(), # noqa: F405
@@ -642,33 +722,20 @@ def CreateIMRInstance(gmpe_name): # noqa: N802, D103
AfshariStewart_2016_AttenRel.NAME # noqa: F405
): AfshariStewart_2016_AttenRel.class_.getName(), # noqa: F405
}
+ # Mapping GMPE name
imrClassName = gmpe_map.get(gmpe_name) # noqa: N806
if imrClassName is None:
return imrClassName
+ # Getting the java class
imrClass = Class.forName(imrClassName) # noqa: N806, F405
- # Try no-arg constructor first
- try:
- ctor = imrClass.getConstructor()
- imr = ctor.newInstance()
- except: # noqa: E722
- # Some GMMs (KS_2006, BommerEtAl_2009, etc.) require a
- # ParameterChangeWarningListener. Create a dummy one.
- try:
- from org.opensha.commons.param.event import ParameterChangeWarningListener
-
- @jpype.JImplements(ParameterChangeWarningListener)
- class _DummyListener:
- @jpype.JOverride
- def parameterChangeWarning(self, event):
- pass
-
- ctor = imrClass.getConstructor(ParameterChangeWarningListener)
- imr = ctor.newInstance(_DummyListener())
- except: # noqa: E722
- return None
+ ctor = imrClass.getConstructor()
+ imr = ctor.newInstance()
+ # Setting default parameters
imr.setParamDefaults()
+ # return
return imr
+
def get_DataSource(paramName, siteData): # noqa: N802, N803, D103
typeMap = SiteTranslator.DATA_TYPE_PARAM_NAME_MAP # noqa: N806, F405
for dataType in typeMap.getTypesForParameterName(paramName): # noqa: N806
@@ -692,11 +759,13 @@ def get_DataSource(paramName, siteData): # noqa: N802, N803, D103
def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
+ # GMPE
try:
imr = CreateIMRInstance(gmpe_name)
except: # noqa: E722
- print('Please check GMPE name.')
+ print('Please check GMPE name.') # noqa: T201
return 1
+ # Site data
sites = ArrayList() # noqa: F405
for cur_site in siteSpec:
cur_loc = Location( # noqa: F405
@@ -708,14 +777,18 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
availableSiteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
except: # noqa: E722
availableSiteData = [] # noqa: N806
- print(
+ print( # noqa: T201
'remote getAllAvailableData is not available temporarily, will use site Vs30 in the site csv file.'
)
+ # return 1
siteTrans = SiteTranslator() # noqa: N806, F405
+ # Looping over all sites
site_prop = []
for i in range(len(siteSpec)):
site_tmp = dict() # noqa: C408
+ # Current site
site = sites.get(i)
+ # Location
cur_site = siteSpec[i]
locResults = { # noqa: N806
'Latitude': cur_site['Location']['Latitude'],
@@ -729,6 +802,7 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
siteDataValues.add(availableSiteData.get(j).getValue(i))
imrSiteParams = imr.getSiteParams() # noqa: N806
siteDataResults = [] # noqa: N806
+ # Setting site parameters
for j in range(imrSiteParams.size()):
siteParam = imrSiteParams.getByIndex(j) # noqa: N806
newParam = Parameter.clone(siteParam) # noqa: N806, F405
@@ -786,10 +860,13 @@ def get_site_prop(gmpe_name, siteSpec): # noqa: C901, N803, D103
}
)
site.addParameter(newParam)
+ # End for j
+ # Updating site specifications
siteSpec[i] = cur_site
site_tmp.update({'Location': locResults, 'SiteData': siteDataResults})
site_prop.append(site_tmp)
+ # Return
return siteSpec, sites, site_prop
@@ -803,16 +880,20 @@ def get_IM( # noqa: C901, N802, D103
station_info,
im_info,
):
+ # GMPE name
gmpe_name = gmpe_info['Type']
+ # Creating intensity measure relationship instance
try:
imr = CreateIMRInstance(gmpe_name)
except: # noqa: E722
- print('Please check GMPE name.')
+ print('Please check GMPE name.') # noqa: T201
return 1, station_info
+ # Getting supported intensity measure types
ims = imr.getSupportedIntensityMeasures()
saParam = ims.getParameter(SA_Param.NAME) # noqa: N806, F405
supportedPeriods = saParam.getPeriodParam().getPeriods() # noqa: N806
Arrays.sort(supportedPeriods) # noqa: F405
+ # Rupture
eqRup = EqkRupture() # noqa: N806, F405
if source_info['Type'] == 'PointSource':
eqRup.setMag(source_info['Magnitude'])
@@ -828,27 +909,38 @@ def get_IM( # noqa: C901, N802, D103
elif source_info['Type'] == 'ERF':
timeSpan = TimeSpan(TimeSpan.NONE, TimeSpan.YEARS) # noqa: N806, F405
erfParams = source_info.get('Parameters', None) # noqa: N806
+ # Additional parameters (if any)
if erfParams is not None:
for k in erfParams.keys:
erf.setParameter(k, erfParams[k])
+ # Time span
timeSpan = erf.getTimeSpan() # noqa: N806
+ # Source
eqSource = erf.getSource(source_info['SourceIndex']) # noqa: N806
eqSource.getName()
+ # Rupture
eqRup = eqSource.getRupture(source_info['RuptureIndex']) # noqa: N806
+ # Properties
magnitude = eqRup.getMag()
averageDip = eqRup.getRuptureSurface().getAveDip() # noqa: N806, F841
averageRake = eqRup.getAveRake() # noqa: N806, F841
+ # Probability
probEqRup = eqRup # noqa: N806
probability = probEqRup.getProbability() # noqa: F841
+ # MAF
meanAnnualRate = probEqRup.getMeanAnnualRate(timeSpan.getDuration()) # noqa: N806
+ # Rupture surface
surface = eqRup.getRuptureSurface() # noqa: F841
+ # Setting up imr
imr.setEqkRupture(eqRup)
imrParams = gmpe_info['Parameters'] # noqa: N806
if bool(imrParams):
- for k in imrParams.keys():
+ for k in imrParams.keys(): # noqa: SIM118
imr.getParameter(k).setValue(imrParams[k])
+ # Station
if station_info['Type'] == 'SiteList':
siteSpec = station_info['SiteList'] # noqa: N806
+ # Intensity measure
periods = im_info.get('Periods', None)
if periods is not None:
periods = supportedPeriods
@@ -863,11 +955,15 @@ def get_IM( # noqa: C901, N802, D103
tag_PGA = True # noqa: N806
if 'PGV' in im_info['Type']:
tag_PGV = True # noqa: N806
+ # Looping over sites
gm_collector = []
for i in range(len(siteSpec)):
gmResults = site_prop[i] # noqa: N806
+ # Current site
site = sites.get(i)
+ # Location
cur_site = siteSpec[i] # noqa: F841
+ # Set up the site in the imr
imr.setSite(site)
try:
stdDevParam = imr.getParameter(StdDevTypeParam.NAME) # noqa: N806, F405
@@ -875,7 +971,7 @@ def get_IM( # noqa: C901, N802, D103
StdDevTypeParam.STD_DEV_TYPE_INTER # noqa: F405
) and stdDevParam.isAllowed(StdDevTypeParam.STD_DEV_TYPE_INTRA) # noqa: F405
except: # noqa: E722
- stdDevParam = None # noqa: N806
+ stdDevParaam = None # noqa: N806, F841
hasIEStats = False # noqa: N806
cur_T = im_info.get('Periods', None) # noqa: N806
if tag_SA:
@@ -904,6 +1000,7 @@ def get_IM( # noqa: C901, N802, D103
saResult['IntraEvStdDev'].append(float(intraEvStdDev))
gmResults.update({'lnSA': saResult})
if tag_PGA:
+ # for PGV current T = 0
cur_T = [0.00] # noqa: N806
pgaResult = {'Mean': [], 'TotalStdDev': []} # noqa: N806
if hasIEStats:
@@ -923,6 +1020,7 @@ def get_IM( # noqa: C901, N802, D103
pgaResult['IntraEvStdDev'].append(float(intraEvStdDev))
gmResults.update({'lnPGA': pgaResult})
if tag_PGV:
+ # for PGV current T = 0
cur_T = [0.00] # noqa: N806
pgvResult = {'Mean': [], 'TotalStdDev': []} # noqa: N806
if hasIEStats:
@@ -943,8 +1041,10 @@ def get_IM( # noqa: C901, N802, D103
gmResults.update({'lnPGV': pgvResult})
gm_collector.append(gmResults)
+ # Updating station information
if station_info['Type'] == 'SiteList':
station_info.update({'SiteList': siteSpec})
+ # Final results
res = {
'Magnitude': magnitude,
'MeanAnnualRate': meanAnnualRate,
@@ -953,16 +1053,22 @@ def get_IM( # noqa: C901, N802, D103
'Periods': cur_T,
'GroundMotions': gm_collector,
}
+ # return
return res, station_info
def get_site_vs30_from_opensha(lat, lon, vs30model='CGS/Wills VS30 Map (2015)'): # noqa: D103
+ # set up site java object
sites = ArrayList() # noqa: F405
num_sites = len(lat)
for i in range(num_sites):
sites.add(Site(Location(lat[i], lon[i]))) # noqa: F405
+
+ # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
+
+ # search name
vs30 = []
for i in range(int(siteData.size())):
cur_siteData = siteData.get(i) # noqa: N806
@@ -971,18 +1077,24 @@ def get_site_vs30_from_opensha(lat, lon, vs30model='CGS/Wills VS30 Map (2015)'):
float(cur_siteData.getValue(x).getValue()) for x in range(num_sites)
]
break
- else:
+ else: # noqa: RET508
continue
- if any([np.isnan(x) for x in vs30]):
+
+ # check if any nan (Wills Map return nan for offshore sites)
+ # Using global vs30 as default patch - 'Global Vs30 from Topographic Slope (Wald & Allen 2008)'
+ if any([np.isnan(x) for x in vs30]): # noqa: C419
non_list = np.where(np.isnan(vs30))[0].tolist()
for i in non_list:
vs30[i] = float(siteData.get(3).getValue(i).getValue())
+
+ # return
return vs30
def get_site_z1pt0_from_opensha(lat, lon): # noqa: D103
sites = ArrayList() # noqa: F405
sites.add(Site(Location(lat, lon))) # noqa: F405
+ # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
for data in siteData:
@@ -996,6 +1108,7 @@ def get_site_z1pt0_from_opensha(lat, lon): # noqa: D103
def get_site_z2pt5_from_opensha(lat, lon): # noqa: D103
sites = ArrayList() # noqa: F405
sites.add(Site(Location(lat, lon))) # noqa: F405
+ # prepare site data java object
siteDataProviders = OrderedSiteDataProviderList.createSiteDataProviderDefaults() # noqa: N806, F405
siteData = siteDataProviders.getAllAvailableData(sites) # noqa: N806
for data in siteData:
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/GMSimulators.py b/modules/performRegionalEventSimulation/regionalGroundMotion/GMSimulators.py
index 2bca9adb4..d2caf8946 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/GMSimulators.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/GMSimulators.py
@@ -41,8 +41,7 @@
# Anne Husley
# Kuanshi Zhong
# Jinyan Zhao
-import sys # noqa: I001
-import time
+import time # noqa: I001
import warnings
import h5py
@@ -384,7 +383,7 @@ def compute_inter_event_residual_ij(self, cm, im_name_list_1, im_name_list_2):
).reshape([len(im_name_list_1), len(im_name_list_2)])
else:
# TODO: extending this to more inter-event correlation models # noqa: TD002
- sys.exit(
+ raise NotImplementedError(
'GM_Simulator.compute_inter_event_residual: currently supporting Baker & Jayaram (2008), Baker & Bradley (2017)'
)
return rho
@@ -495,7 +494,7 @@ def compute_intra_event_residual_i(self, cm, im_name_list, num_simu): # noqa: D
)
else:
# TODO: extending this to more inter-event correlation models # noqa: TD002
- sys.exit(
+ raise NotImplementedError(
'GM_Simulator.compute_intra_event_residual: currently supporting Jayaram & Baker (2009), Loth & Baker (2013),Markhvida et al. (2017), Du & Ning (2021)'
)
return residuals
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/GlobalVariable.py b/modules/performRegionalEventSimulation/regionalGroundMotion/GlobalVariable.py
index cbc4addbe..56be2e714 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/GlobalVariable.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/GlobalVariable.py
@@ -1 +1,107 @@
-JVM_started = False # noqa: INP001, D100
+"""Shared JVM state and helpers for the regional ground-motion module.""" # noqa: INP001
+
+import os
+import platform
+import subprocess
+import sys
+
+# JPype allows only one JVM per Python process; entry-point scripts check
+# this flag before calling jpype.startJVM().
+JVM_started = False
+
+
+def _find_jvm_macos():
+ """Locate an arch-matched JVM on macOS via /usr/libexec/java_home -a."""
+ py_arch = platform.machine()
+ if py_arch not in ('arm64', 'x86_64'):
+ return None
+
+ try:
+ result = subprocess.run( # noqa: S603
+ ['/usr/libexec/java_home', '-a', py_arch],
+ capture_output=True,
+ text=True,
+ timeout=10,
+ check=False,
+ )
+ except (FileNotFoundError, subprocess.TimeoutExpired, OSError):
+ return None
+
+ if result.returncode != 0:
+ return None
+
+ jvm_home = result.stdout.strip()
+ if not jvm_home:
+ return None
+
+ libjli = os.path.join(jvm_home, 'lib', 'libjli.dylib') # noqa: PTH118
+ return libjli if os.path.exists(libjli) else None # noqa: PTH110
+
+
+def _find_jvm_windows():
+ """Locate a JPype-compatible JVM on Windows.
+
+ Honors ``JAVA_HOME`` first, then scans the standard install folders
+ under ``Program Files``. When several JDKs are installed, the one with
+ the highest-sorting directory name wins (so JDK 17 beats JDK 16).
+ Returns None if nothing usable is found.
+ """
+
+ def _jvm_dll(jdk_home):
+ candidate = os.path.join(jdk_home, 'bin', 'server', 'jvm.dll') # noqa: PTH118
+ return candidate if os.path.exists(candidate) else None # noqa: PTH110
+
+ java_home = os.environ.get('JAVA_HOME')
+ if java_home:
+ dll = _jvm_dll(java_home)
+ if dll:
+ return dll
+
+ is_64bit = sys.maxsize > 2**32
+ program_files = os.environ.get(
+ 'ProgramFiles' if is_64bit else 'ProgramFiles(x86)',
+ r'C:\Program Files',
+ )
+
+ found = []
+ for parent in (
+ os.path.join(program_files, 'Eclipse Adoptium'), # noqa: PTH118
+ os.path.join(program_files, 'Java'), # noqa: PTH118
+ os.path.join(program_files, 'Microsoft', 'jdk'), # noqa: PTH118
+ ):
+ if not os.path.isdir(parent): # noqa: PTH112
+ continue
+ try:
+ entries = os.listdir(parent)
+ except OSError:
+ continue
+ for entry in entries:
+ dll = _jvm_dll(os.path.join(parent, entry)) # noqa: PTH118
+ if dll:
+ found.append((entry, dll))
+
+ if not found:
+ return None
+
+ # Modern Adoptium directory names embed the version
+ # (e.g. "jdk-17.0.10.7-hotspot"), so reverse-sorted alphabetically
+ # is also reverse-sorted by version.
+ found.sort(reverse=True)
+ return found[0][1]
+
+
+def find_compatible_jvm_path():
+ """Return a JVM library path matching the running Python interpreter.
+
+ Guards against JPype's default lookup picking a JVM whose architecture
+ does not match Python (common on macOS when both arm64 and x86_64 JDKs
+ are installed) or a stale JAVA_HOME on Windows. Returns None on Linux
+ and other platforms, and as a safe fallback when no usable JVM is
+ found; callers should pass the result through to
+ ``jpype.startJVM(jvmpath=...)`` unconditionally.
+ """
+ if sys.platform == 'darwin':
+ return _find_jvm_macos()
+ if sys.platform == 'win32':
+ return _find_jvm_windows()
+ return None
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardOccurrence.py b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardOccurrence.py
index 597644943..bc512f734 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardOccurrence.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardOccurrence.py
@@ -39,9 +39,8 @@
import collections
import itertools
-import json as python_json
+import json
import os
-import sys
import threading
import time
@@ -52,7 +51,7 @@
from scipy.stats import norm
from sklearn.linear_model import lasso_path
from tqdm import tqdm
-from USGS_API import * # noqa: F403
+from USGS_API import USGS_HazardCurve
def configure_hazard_occurrence( # noqa: C901, D103
@@ -166,7 +165,7 @@ def configure_hazard_occurrence( # noqa: C901, D103
cur_imt = im_type
if IMfile.lower().endswith('.json'):
with open(IMfile) as f: # noqa: PTH123
- IMdata = python_json.load(f) # noqa: N806
+ IMdata = json.load(f) # noqa: N806
hc_data = calc_hazard_curves(IMdata, site_config, cur_imt)
elif IMfile.lower().endswith('.hdf5'):
hc_data = calc_hazard_curves_hdf5(
@@ -208,7 +207,7 @@ def configure_hazard_occurrence( # noqa: C901, D103
}
# output the hazard occurrence information file
with open(os.path.join(output_dir, 'HazardCurves.json'), 'w') as f: # noqa: PTH118, PTH123
- python_json.dump(occ_dict, f, indent=2)
+ json.dump(occ_dict, f, indent=2)
occ_dict = {
'Model': model_type,
'NumTargetEQs': num_target_eqs,
@@ -383,7 +382,7 @@ def get_hazard_curves(input_dir=None, input_csv=None, input_json=None): # noqa:
if input_json is not None: # noqa: RET503
with open(input_json) as f: # noqa: PTH123
- hc_data = python_json.load(f)
+ hc_data = json.load(f)
return hc_data # noqa: RET504
@@ -405,7 +404,7 @@ def get_im_exceedance_probility( # noqa: C901, D103
# initialize output
if IMfile.lower().endswith('.json'):
with open(IMfile) as f: # noqa: PTH123
- im_raw = python_json.load(f)
+ im_raw = json.load(f)
num_sites = len(im_raw[scenario_idx[0]].get('GroundMotions'))
elif IMfile.lower().endswith('.hdf5'):
with h5py.File(IMfile, 'r') as f:
@@ -571,7 +570,7 @@ def export_sampled_earthquakes(error, id_selected_eqs, eqdata, P, output_dir=Non
if output_dir is not None:
with open(os.path.join(output_dir, 'RupSampled.json'), 'w') as f: # noqa: PTH118, PTH123
- python_json.dump(dict_selected_eqs, f, indent=2)
+ json.dump(dict_selected_eqs, f, indent=2)
# def export_sampled_earthquakes(occ_dict, im_raw, site_config, id_selected_eqs, eqdata, P, output_dir=None):
@@ -613,7 +612,7 @@ def export_sampled_earthquakes(error, id_selected_eqs, eqdata, P, output_dir=Non
# if output_dir is not None:
# with open(os.path.join(output_dir,'RupSampled.json'), 'w') as f:
-# python_json.dump(dict_selected_eqs, f, indent=2)
+# json.dump(dict_selected_eqs, f, indent=2)
class OccurrenceModel_ManzourDavidson2016: # noqa: D101
@@ -807,7 +806,7 @@ def export_sampled_gmms( # noqa: D102
if output_dir is not None:
with open(os.path.join(output_dir, 'InfoSampledGM.json'), 'w') as f: # noqa: PTH118, PTH123
- python_json.dump(dict_selected_gmms, f, indent=2)
+ json.dump(dict_selected_gmms, f, indent=2)
class OccurrenceModel_Wangetal2023: # noqa: D101
@@ -958,7 +957,7 @@ def get_selected_earthquake(self): # noqa: D102
)
if self.num_selected[self.selected_alpha_ind] == 0:
- sys.exit(
+ raise RuntimeError(
'ERROR: Zero scenarios/ground motions are selected in Wang et al. (2023).\n' # noqa: ISC003
+ f'The tunnling parameter used is {self.alphas[self.selected_alpha_ind]}.\n'
+ 'Try using a smaller tuning parameter.'
@@ -992,4 +991,4 @@ def export_sampled_gmms( # noqa: D102
if output_dir is not None:
with open(os.path.join(output_dir, 'InfoSampledGM.json'), 'w') as f: # noqa: PTH118, PTH123
- python_json.dump(dict_selected_gmms, f, indent=2)
+ json.dump(dict_selected_gmms, f, indent=2)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulation.py b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulation.py
index a400da3a3..e225e29d2 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulation.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulation.py
@@ -38,12 +38,9 @@
#
import argparse
-import importlib
import json
import os
import shutil
-import subprocess
-import sys
import time
import numpy as np
@@ -51,7 +48,7 @@
R2D = True
-OPENSHA_JAR = 'opensha-all.jar' # version 25.4.1 (invoked in 04/2026)
+OPENSHA_JAR = 'opensha-all.jar' # version 26.1.1 (invoked in 05/2026)
def site_job(hazard_info): # noqa: C901, D103
@@ -214,7 +211,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
)
if not filePath_ini:
# Error in ini file
- sys.exit(
+ raise RuntimeError(
'HazardSimulation: errors in preparing the OpenQuake configuration file.'
)
if scenario_info['EqRupture']['Type'] in [
@@ -240,7 +237,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
)
)
print(err_msg) # noqa: T201
- sys.exit(err_msg)
+ raise RuntimeError(err_msg)
else:
print('HazardSimulation: OpenQuake Classical PSHA completed.') # noqa: T201
if scenario_info['EqRupture'].get('UHS', False):
@@ -264,7 +261,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
print('HazardSimulation: OpenQuake Scenario calculation completed.') # noqa: T201
else:
- sys.exit(
+ raise NotImplementedError(
'HazardSimulation: OpenQuakeClassicalPSHA, OpenQuakeUserConfig and OpenQuakeScenario are supported.'
)
@@ -372,7 +369,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
id_selected_scens = [int(x / num_gm_per_site) for x in id_selected_gmms]
id_selected_simus = [x % num_gm_per_site for x in id_selected_gmms]
# export sampled earthquakes
- _ = export_sampled_gmms( # noqa: F405
+ _ = export_sampled_gmms( # noqa: F821 -- export_sampled_gmms is a method, not a free function
id_selected_gmms, id_selected_scens, P_gmm, output_dir
)
num_site = ln_im_mr[0].shape[0]
@@ -439,7 +436,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
if runtag:
print('HazardSimulation: the ground motion list saved.') # noqa: T201
else:
- sys.exit(
+ raise RuntimeError(
'HazardSimulation: warning - issues with saving the ground motion list.'
)
# Downloading records
@@ -447,14 +444,14 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
user_password = event_info.get('UserPassword', None)
if (user_name is not None) and (user_password is not None) and (not R2D):
print('HazardSimulation: downloading ground motion records.') # noqa: T201
- raw_dir = download_ground_motion( # noqa: F405
+ raw_dir = download_ground_motion( # noqa: F821 -- see import block
gm_id, user_name, user_password, output_dir
)
if raw_dir:
print('HazardSimulation: ground motion records downloaded.') # noqa: T201
# Parsing records
print('HazardSimulation: parsing records.') # noqa: T201
- record_dir = parse_record( # noqa: F405, F841
+ record_dir = parse_record( # noqa: F821, F841 -- see import block
gm_file,
raw_dir,
output_dir,
@@ -524,33 +521,25 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
except: # noqa: E722
oq_flag = False
- # dependencies
- if R2D:
- packages = ['tqdm', 'psutil', 'PuLP', 'requests']
- else:
- packages = ['selenium', 'tqdm', 'psutil', 'PuLP', 'requests']
- for p in packages:
- if importlib.util.find_spec(p) is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', p]) # noqa: S603
-
# set up environment
import socket
if 'stampede2' not in socket.gethostname():
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'JPype1']) # noqa: S603
+ import GlobalVariable
+
import jpype
- from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
- #jpype.addClassPath('./lib/OpenSHA-1.5.2.jar') # not supported by opensha starting from 02/2026
jpype.addClassPath('./lib/{}'.format(OPENSHA_JAR))
try:
- jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
+ jpype.startJVM(
+ f'-Xmx{memory_request}G',
+ convertStrings=False,
+ jvmpath=GlobalVariable.find_compatible_jvm_path(),
+ )
except: # noqa: E722
print( # noqa: T201
- #f'StartJVM of ./lib/OpenSHA-1.5.2.jar with {memory_request} GB Memory fails. Try again after releasing some memory'
f'StartJVM of ./lib/{OPENSHA_JAR} with {memory_request} GB Memory fails. Try again after releasing some memory'
)
if oq_flag:
@@ -574,18 +563,37 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
shutil.rmtree(os.environ.get('OQ_DATADIR'))
os.makedirs(f"{os.environ.get('OQ_DATADIR')}") # noqa: PTH103
- # import modules
- # from ComputeIntensityMeasure import *
- # from CreateScenario import *
- # from CreateStation import *
-
- # # KZ-08/23/22: adding hazard occurrence model
- # from HazardOccurrence import *
- # from SelectGroundMotion import *
+ from ComputeIntensityMeasure import compute_im, export_im
+ from CreateScenario import (
+ create_earthquake_scenarios,
+ create_wind_scenarios,
+ load_earthquake_scenarios,
+ )
+ from CreateStation import create_stations
+ from GMSimulators import simulate_ground_motion
+ from HazardOccurrence import (
+ configure_hazard_occurrence,
+ export_sampled_earthquakes,
+ get_im_exceedance_probability_gm,
+ get_im_exceedance_probility,
+ sample_earthquake_occurrence,
+ )
+ from SelectGroundMotion import (
+ output_all_ground_motion_info,
+ select_ground_motion,
+ )
+ # The NGAWest2 record-download branch below calls download_ground_motion
+ # and parse_record, which are not defined in SelectGroundMotion's public
+ # surface; the F821 noqa annotations at those sites mark that branch as
+ # known-dead.
if oq_flag:
- # import FetchOpenQuake
- from FetchOpenQuake import * # noqa: F403
+ from FetchOpenQuake import (
+ OpenQuakeHazardCalc,
+ openquake_config,
+ oq_read_uhs_classical_psha,
+ oq_run_classical_psha,
+ )
# Initial process list
import psutil
@@ -603,6 +611,3 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
site_job(hazard_info)
else:
print('HazardSimulation: --job_type = Hazard or Site (please check).') # noqa: T201
-
- # Closing the current process
- sys.exit(0)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ.py b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ.py
index eb3333e40..8e67f930c 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ.py
@@ -38,12 +38,9 @@
#
import argparse
-import importlib
import json
import os
import shutil
-import subprocess
-import sys
import time
import numpy as np
@@ -52,7 +49,7 @@
R2D = True
-OPENSHA_JAR = 'opensha-all.jar' # version 25.4.1 (invoked in 04/2026)
+OPENSHA_JAR = 'opensha-all.jar' # version 26.1.1 (invoked in 05/2026)
def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
from CreateScenario import load_ruptures_openquake
@@ -145,7 +142,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
)
if not filePath_ini:
# Error in ini file
- sys.exit(
+ raise RuntimeError(
'HazardSimulation: errors in preparing the OpenQuake configuration file.'
)
if scenario_info['EqRupture']['Type'] in [
@@ -155,11 +152,11 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
]:
# Calling openquake to run classical PSHA
# oq_version = scenario_info['EqRupture'].get('OQVersion',default_oq_version)
- oq_run_flag = oq_run_classical_psha( # noqa: F405
+ oq_run_flag = oq_run_classical_psha(
filePath_ini,
exports='csv',
oq_version=oq_ver_loaded,
- dir_info=dir_info, # noqa: F405
+ dir_info=dir_info, # noqa: F821 -- dir_info is never defined in this scope
)
if oq_run_flag:
err_msg = 'HazardSimulation: OpenQuake Classical PSHA failed.'
@@ -171,14 +168,14 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
)
)
print(err_msg) # noqa: T201
- sys.exit(err_msg)
+ raise RuntimeError(err_msg)
else:
print('HazardSimulation: OpenQuake Classical PSHA completed.') # noqa: T201
if scenario_info['EqRupture'].get('UHS', False):
- ln_im_mr, mag_maf, im_list = oq_read_uhs_classical_psha( # noqa: F405
+ ln_im_mr, mag_maf, im_list = oq_read_uhs_classical_psha(
scenario_info,
event_info,
- dir_info, # noqa: F405
+ dir_info, # noqa: F821 -- dir_info is never defined in this scope
)
else:
ln_im_mr = []
@@ -188,7 +185,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
elif scenario_info['EqRupture']['Type'] == 'oqSourceXML':
# Creating and conducting OpenQuake calculations
- oq_calc = OpenQuakeHazardCalc( # noqa: F405
+ oq_calc = OpenQuakeHazardCalc(
filePath_ini,
event_info,
oq_ver_loaded,
@@ -200,7 +197,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
print('HazardSimulation: OpenQuake Scenario calculation completed.') # noqa: T201
else:
- sys.exit(
+ raise NotImplementedError(
'HazardSimulation: OpenQuakeClassicalPSHA, OpenQuakeUserConfig and OpenQuakeScenario are supported.'
)
@@ -407,7 +404,7 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
if runtag:
print('HazardSimulation: the ground motion list saved.') # noqa: T201
else:
- sys.exit(
+ raise RuntimeError(
'HazardSimulation: warning - issues with saving the ground motion list.'
)
# Downloading records
@@ -415,14 +412,14 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
user_password = event_info.get('UserPassword', None)
if (user_name is not None) and (user_password is not None) and (not R2D):
print('HazardSimulation: downloading ground motion records.') # noqa: T201
- raw_dir = download_ground_motion( # noqa: F405
+ raw_dir = download_ground_motion( # noqa: F821 -- see import block
gm_id, user_name, user_password, output_dir
)
if raw_dir:
print('HazardSimulation: ground motion records downloaded.') # noqa: T201
# Parsing records
print('HazardSimulation: parsing records.') # noqa: T201
- record_dir = parse_record( # noqa: F405, F841
+ record_dir = parse_record( # noqa: F821, F841 -- see import block
gm_file,
raw_dir,
output_dir,
@@ -550,15 +547,6 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
except: # noqa: E722
oq_flag = False
- # dependencies
- if R2D:
- packages = ['tqdm', 'psutil', 'PuLP', 'requests']
- else:
- packages = ['selenium', 'tqdm', 'psutil', 'PuLP', 'requests']
- for p in packages:
- if importlib.util.find_spec(p) is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', p]) # noqa: S603
-
# set up environment
import socket
@@ -567,21 +555,17 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
if GlobalVariable.JVM_started is False:
GlobalVariable.JVM_started = True
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call( # noqa: S603
- [sys.executable, '-m', 'pip', 'install', 'JPype1']
- ) # noqa: RUF100, S603
import jpype
-
- # from jpype import imports
import jpype.imports
- from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
- #jpype.addClassPath('./lib/OpenSHA-1.5.2.jar') # not supported by opensha starting from 02/2026
jpype.addClassPath('./lib/{}'.format(OPENSHA_JAR))
- jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
+ jpype.startJVM(
+ f'-Xmx{memory_request}G',
+ convertStrings=False,
+ jvmpath=GlobalVariable.find_compatible_jvm_path(),
+ )
if oq_flag:
# clear up old db.sqlite3 if any
if os.path.isfile(os.path.expanduser('~/oqdata/db.sqlite3')): # noqa: PTH111, PTH113
@@ -603,18 +587,33 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
shutil.rmtree(os.environ.get('OQ_DATADIR'))
os.makedirs(f"{os.environ.get('OQ_DATADIR')}") # noqa: PTH103
- # import modules
- from ComputeIntensityMeasure import * # noqa: F403
- from CreateScenario import * # noqa: F403
- from CreateStation import * # noqa: F403
+ from ComputeIntensityMeasure import compute_im, export_im
+ from CreateScenario import load_earthquake_rupFile
# KZ-08/23/22: adding hazard occurrence model
- from HazardOccurrence import * # noqa: F403
- from SelectGroundMotion import * # noqa: F403
+ from HazardOccurrence import (
+ configure_hazard_occurrence,
+ export_sampled_earthquakes,
+ get_im_exceedance_probability_gm,
+ get_im_exceedance_probility,
+ sample_earthquake_occurrence,
+ )
+ from SelectGroundMotion import (
+ output_all_ground_motion_info,
+ select_ground_motion,
+ )
+ # The NGAWest2 record-download branch below calls download_ground_motion
+ # and parse_record, which are not defined in SelectGroundMotion's public
+ # surface; the F821 noqa annotations at those sites mark that branch as
+ # known-dead.
if oq_flag:
- # import FetchOpenQuake
- from FetchOpenQuake import * # noqa: F403
+ from FetchOpenQuake import (
+ OpenQuakeHazardCalc,
+ openquake_config,
+ oq_read_uhs_classical_psha,
+ oq_run_classical_psha,
+ )
# untar site databases
# site_database = ['global_vs30_4km.tar.gz','global_zTR_4km.tar.gz','thompson_vs30_4km.tar.gz']
@@ -633,6 +632,3 @@ def hazard_job(hazard_info): # noqa: C901, D103, PLR0915
]
hazard_job(hazard_info)
-
- # Closing the current process
- sys.exit(0)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ_Stampede3.py b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ_Stampede3.py
index f10d4d454..407ef46af 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ_Stampede3.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/HazardSimulationEQ_Stampede3.py
@@ -1,6 +1,5 @@
import ujson as json
import os
-import sys
import pickle
from tqdm import tqdm
import numpy as np
@@ -145,7 +144,7 @@ def is_valid_path_format(string):
if returncode != 0:
print(result)
- sys.exit(f'return code: {returncode}')
+ raise RuntimeError(f'return code: {returncode}')
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/LoadRupFile.py b/modules/performRegionalEventSimulation/regionalGroundMotion/LoadRupFile.py
index b73cab8c7..761ce09ea 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/LoadRupFile.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/LoadRupFile.py
@@ -2,7 +2,6 @@
from tqdm import tqdm
import geopandas as gpd
import shapely
-import sys
import ujson as json
def get_rups_to_run(scenario_info, num_scenarios): # noqa: C901, D103
@@ -30,7 +29,7 @@ def get_rups_to_run(scenario_info, num_scenarios): # noqa: C901, D103
np.where(np.isin(rups_requested, rups_available))[0]
]
else:
- sys.exit(
+ raise NotImplementedError(
f'The scenario selection method {scenario_info["Generator"].get("method", None)} is not available'
)
return rups_to_run
@@ -42,11 +41,13 @@ def load_earthquake_rupFile(scenario_info, rupFilePath): # noqa: N802, N803, D1
with open(rupFilePath) as f: # noqa: PTH123
user_scenarios = json.load(f)
except: # noqa: E722
- sys.exit(f'CreateScenario: source file {rupFilePath} not found.')
+ raise FileNotFoundError( # noqa: B904
+ f'CreateScenario: source file {rupFilePath} not found.'
+ )
# number of features (i.e., ruptures)
num_scenarios = len(user_scenarios.get('features', []))
if num_scenarios < 1:
- sys.exit('CreateScenario: source file is empty.')
+ raise ValueError('CreateScenario: source file is empty.')
rups_to_run = get_rups_to_run(scenario_info, num_scenarios)
# get rupture and source ids
scenario_data = {}
@@ -115,7 +116,7 @@ def load_earthquake_rup_scenario(scenario_info, user_scenarios): # noqa: N802,
# number of features (i.e., ruptures)
num_scenarios = len(user_scenarios)
if num_scenarios < 1:
- sys.exit('CreateScenario: source file is empty.')
+ raise ValueError('CreateScenario: source file is empty.')
rups_to_run = get_rups_to_run(scenario_info, num_scenarios)
# get rupture and source ids
scenario_data = {}
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/ScenarioForecast.py b/modules/performRegionalEventSimulation/regionalGroundMotion/ScenarioForecast.py
index 46207bfc0..a80670221 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/ScenarioForecast.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/ScenarioForecast.py
@@ -38,16 +38,13 @@
# Jinyan Zhao
import argparse
-import importlib
import json
import os
-import subprocess
-import sys
import tarfile
import psutil
-OPENSHA_JAR = 'opensha-all.jar' # version 25.4.1 (invoked in 04/2026)
+OPENSHA_JAR = 'opensha-all.jar' # version 26.1.1 (invoked in 05/2026)
if __name__ == '__main__':
# parse arguments
@@ -84,16 +81,6 @@
except: # noqa: E722
oq_flag = False
- # dependencies
- packages = ['tqdm', 'psutil', 'pulp', 'requests']
- for p in packages:
- if importlib.util.find_spec(p) is None:
- # print(f"""The Python package {p} is required but not found.
- # Please install it by running
- # "{sys.executable} -m pip install -q {p}"
- # in your terminal or command prompt""")
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', '-q', p]) # noqa: S603
-
# set up environment
import socket
@@ -102,21 +89,17 @@
if GlobalVariable.JVM_started is False:
GlobalVariable.JVM_started = True
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call( # noqa: S603
- [sys.executable, '-m', 'pip', 'install', 'JPype1']
- ) # noqa: RUF100, S603
import jpype
-
- # from jpype import imports
import jpype.imports
- from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
- #jpype.addClassPath('./lib/OpenSHA-1.5.2.jar') # not supported by opensha starting from 02/2026
jpype.addClassPath('./lib/{}'.format(OPENSHA_JAR))
- jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
+ jpype.startJVM(
+ f'-Xmx{memory_request}G',
+ convertStrings=False,
+ jvmpath=GlobalVariable.find_compatible_jvm_path(),
+ )
from CreateScenario import (
create_earthquake_scenarios,
create_wind_scenarios,
@@ -139,10 +122,6 @@
# shutil.rmtree(os.environ.get('OQ_DATADIR'))
# os.makedirs(f"{os.environ.get('OQ_DATADIR')}")
- if oq_flag:
- # import FetchOpenQuake
- from FetchOpenQuake import * # noqa: F403
-
# untar site databases
site_database = [
'global_vs30_4km.tar.gz',
@@ -227,6 +206,3 @@
print('HazardSimulation: currently only supports EQ and Wind simulations.') # noqa: T201
# print(scenarios)
print('HazardSimulation: scenarios created.') # noqa: T201
-
- # Closing the current process
- sys.exit(0)
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/SelectGroundMotion.py b/modules/performRegionalEventSimulation/regionalGroundMotion/SelectGroundMotion.py
index 884a2d1cc..97268295c 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/SelectGroundMotion.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/SelectGroundMotion.py
@@ -233,7 +233,9 @@ def select_ground_motion( # noqa: C901, D103
filename.extend(tmp_filename)
# print(tmp_min_err)
else:
- sys.exit('SelectGroundMotion: currently only supporting NGAWest2.')
+ raise NotImplementedError(
+ 'SelectGroundMotion: currently only supporting NGAWest2.'
+ )
# output data
station_name = ['site' + str(j) + '.csv' for j in range(len(stations))]
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/calculation_single_proc.py b/modules/performRegionalEventSimulation/regionalGroundMotion/calculation_single_proc.py
index fe104580d..3ffbb309d 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/calculation_single_proc.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/calculation_single_proc.py
@@ -1,10 +1,11 @@
+# WARNING: this script targets the legacy OpenSHA v1.5.2 jar and is not
+# invoked by R2D. To use it, switch to opensha-all.jar and verify
+# compatibility with FetchOpenSHA.py.
import ujson as json
import os
import sys
import psutil
-import importlib
-import subprocess
-# Add the folder containing the script to the system path
+
script_path = os.path.dirname(os.path.realpath(__file__))
opensha_path = os.path.join(script_path, 'lib', 'OpenSHA-1.5.2.jar')
sys.path.append(script_path)
@@ -17,15 +18,8 @@
import GlobalVariable
if GlobalVariable.JVM_started is False:
GlobalVariable.JVM_started = True
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call( # noqa: S603
- [sys.executable, '-m', 'pip', 'install', 'JPype1']
- ) # noqa: RUF100, S603
import jpype
-
- # from jpype import imports
import jpype.imports
- # from jpype.types import * # noqa: F403
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.25)
@@ -33,34 +27,25 @@
jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
print('JVM started')
from ComputeIntensityMeasure import compute_im # noqa: F403
-# print('debug 3')
from GMSimulators import simulate_ground_motion
-# print('debug 1')
from LoadRupFile import load_earthquake_rup_scenario # noqa: F403
-# print('debug 1')
from SelectGroundMotion import select_ground_motion
from SelectGroundMotion import output_all_ground_motion_info
-# print('debug 1')
from ComputeIntensityMeasure import export_im
-# print('debug 1')
-# # @profile
+
def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
-
hazard_info = copy.deepcopy(hazard_info_orig)
- # Select the range of scenarios to simulate
- hazard_info['Scenario']['Generator']['RuptureFilter'] = str(sce_idx+1) # The index in the r2d ground motion simulation starts from 1
-
- # Below are initializing HazardSimulationEQ
- # directory (back compatibility here)
+ # R2D rupture indices are 1-based.
+ hazard_info['Scenario']['Generator']['RuptureFilter'] = str(sce_idx+1)
+
work_dir = hazard_info['Directory']
output_dir = os.path.join(work_dir, f'Output_{sce_idx}')
try:
os.mkdir(f'{output_dir}') # noqa: PTH102
except: # noqa: E722
print('HazardSimulation: output folder already exists.') # noqa: T201
- # Read Site .csv
site_file = hazard_info['Site']['siteFile']
try:
stations = pd.read_csv(site_file).to_dict(orient='records')
@@ -73,16 +58,10 @@ def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
scenarios = load_earthquake_rup_scenario(scenario_info, all_rups) # noqa: F405
print(f'Proc {procID}: HazardSimulation: scenarios loaded.') # noqa: T201
- # current, peak = tracemalloc.get_traced_memory()
- # print(f"Current memory usage: {current / 10**6} MB")
- # print(f"Peak memory usage: {peak / 10**6} MB")
-
selected_scen_ids = sorted(list(scenarios.keys())) # noqa: C414
- # Computing intensity measures
print(f'Proc {procID}: HazardSimulation: computing intensity measures.') # noqa: T201
- # Computing uncorrelated Sa
event_info = hazard_info['Event']
- # When vector IM is used. The PGA/SA needs to be computed before PGV
+ # PGA/SA must be computed before PGV when a vector IM is used.
im_info = event_info['IntensityMeasure']
im_raw_path, im_list = compute_im( # noqa: F405
scenarios,
@@ -94,13 +73,8 @@ def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
output_dir,
mth_flag=False,
)
-
- # current, peak = tracemalloc.get_traced_memory()
- # print(f"Current memory usage: {current / 10**6} MB")
- # print(f"Peak memory usage: {peak / 10**6} MB")
num_gm_per_site = event_info['NumberPerSite']
- # Computing correlated IMs
ln_im_mr, mag_maf = simulate_ground_motion(
stations,
im_raw_path,
@@ -112,9 +86,7 @@ def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
selected_scen_ids,
)
- # Selecting ground motion records
if scenario_info['Type'] == 'Earthquake':
- # Selecting records
data_source = event_info.get('Database', 0)
if data_source:
print('HazardSimulation: selecting ground motion records.') # noqa: T201
@@ -135,7 +107,6 @@ def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
print( # noqa: T201
f'HazardSimulation: ground motion records selected ({time.time() - start_time} s).'
)
- # print(gm_id)
gm_id = [int(i) for i in np.unique(gm_id)]
gm_file = [i for i in np.unique(gm_file)] # noqa: C416
runtag = output_all_ground_motion_info( # noqa: F405
@@ -144,7 +115,7 @@ def run_scenario_i(sce_idx, hazard_info_orig, all_rups, procID = 0):
if runtag:
print('HazardSimulation: the ground motion list saved.') # noqa: T201
else:
- sys.exit(
+ raise RuntimeError(
'HazardSimulation: warning - issues with saving the ground motion list.'
)
else:
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py b/modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py
index 2d2a0a45a..8bf5e4d79 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/gmpe/openSHAGMPE.py
@@ -38,7 +38,6 @@
# Transferred from openSHA to achieve better performance in r2d
import os
-import sys
import time
import numpy as np
@@ -81,10 +80,9 @@ def setIMT(self, imt): # noqa: N802, D102
supported_imt = [
f'SA{x}s' if isinstance(x, float) else x for x in self.supportedImt
]
- sys.exit(
+ raise ValueError(
f'The IM type {imt} is not supported by Chiou and Young (2014). \n The supported IM types are {supported_imt}'
)
- return False
self.c1 = self.coeff['c1'][imt]
self.c1a = self.coeff['c1a'][imt]
self.c1b = self.coeff['c1b'][imt]
@@ -350,10 +348,9 @@ def setIMT(self, imt): # noqa: N802, D102
supported_imt = [
f'SA{x}s' if isinstance(x, float) else x for x in self.supportedImt
]
- sys.exit(
+ raise ValueError(
f'The IM type {imt} is not supported by Abrahamson, Silva, and Kamai (2014). \n The supported IM types are {supported_imt}'
)
- return
self.imt = imt
self.a1 = self.coeff['a1'][imt]
self.a2 = self.coeff['a2'][imt]
@@ -663,10 +660,9 @@ def setIMT(self, imt): # noqa: N802, D102
supported_imt = [
f'SA{x}s' if isinstance(x, float) else x for x in self.supportedImt
]
- sys.exit(
+ raise ValueError(
f'The IM type {imt} is not supported by Boore, Stewart, Seyhan & Atkinson (2014). \n The supported IM types are {supported_imt}'
)
- return
self.imt = imt
self.e0 = self.coeff['e0'][imt]
self.e1 = self.coeff['e1'][imt]
@@ -893,10 +889,9 @@ def setIMT(self, imt): # noqa: N802, D102
supported_imt = [
f'SA{x}s' if isinstance(x, float) else x for x in self.supportedImt
]
- sys.exit(
+ raise ValueError(
f'The IM type {imt} is not supported by Campbell & Bozorgnia (2014). \n The supported IM types are {supported_imt}'
)
- return
self.imt = imt
self.c0 = self.coeff['c0'][imt]
self.c1 = self.coeff['c1'][imt]
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/landslide.py b/modules/performRegionalEventSimulation/regionalGroundMotion/landslide.py
index 288ee7e29..5fea551b4 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/landslide.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/landslide.py
@@ -27,7 +27,9 @@ def sampleRaster( # noqa: N802
f'More than one band in the file {raster_file_path}, the first band is used.'
)
except: # noqa: E722
- sys.exit(f'Can not read data from {raster_file_path}')
+ raise OSError( # noqa: B904
+ f'Can not read data from {raster_file_path}'
+ )
if xy_crs != raster_crs:
# make transformer for reprojection
transformer_xy_to_data = Transformer.from_crs(
@@ -84,7 +86,7 @@ def sampleVector(vector_file_path, vector_crs, x, y, dtype=None): # noqa: ARG00
xy_crs = CRS.from_user_input(4326)
vector_gdf = gpd.read_file(vector_file_path)
if vector_gdf.crs != vector_crs:
- sys.exit(
+ raise ValueError(
f"The CRS of vector file {vector_file_path} is {vector_gdf.crs}, and doesn't match the input CRS ({xy_crs}) defined for liquefaction triggering models"
)
if xy_crs != vector_crs:
@@ -522,11 +524,9 @@ def run( # noqa: D102
ln_im_data[scenario_id] = im_data_scen
im_list = im_list + output_keys # noqa: PLR6104, RUF100
else:
- sys.exit(
+ raise ValueError(
f"'PGA' is missing in the selected intensity measures and the landslide model 'BrayMacedo2019' can not be computed." # noqa: F541
)
- # print(f"At least one of 'PGA' and 'PGV' is missing in the selected intensity measures and the liquefaction trigging model 'ZhuEtal2017' can not be computed."\
- # , file=sys.stderr)
# sys.stderr.write("test")
# sys.exit(-1)
return (
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/lib/OpenSHA-1.5.2.jar b/modules/performRegionalEventSimulation/regionalGroundMotion/lib/OpenSHA-1.5.2.jar
deleted file mode 100644
index c18f07151..000000000
Binary files a/modules/performRegionalEventSimulation/regionalGroundMotion/lib/OpenSHA-1.5.2.jar and /dev/null differ
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/lib/opensha-all.jar b/modules/performRegionalEventSimulation/regionalGroundMotion/lib/opensha-all.jar
index 74ac517e0..607967390 100644
Binary files a/modules/performRegionalEventSimulation/regionalGroundMotion/lib/opensha-all.jar and b/modules/performRegionalEventSimulation/regionalGroundMotion/lib/opensha-all.jar differ
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/liquefaction.py b/modules/performRegionalEventSimulation/regionalGroundMotion/liquefaction.py
index 5856190f4..aae81183c 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/liquefaction.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/liquefaction.py
@@ -36,7 +36,9 @@ def sampleRaster( # noqa: N802
f'More than one band in the file {raster_file_path}, the first band is used.'
)
except: # noqa: E722
- sys.exit(f'Can not read data from {raster_file_path}')
+ raise OSError( # noqa: B904
+ f'Can not read data from {raster_file_path}'
+ )
if xy_crs != raster_crs:
# make transformer for reprojection
transformer_xy_to_data = Transformer.from_crs(
@@ -96,7 +98,7 @@ def sampleVector(vector_file_path, vector_crs, x, y, dtype=None): # noqa: ARG00
try:
user_crs_input = CRS.from_user_input(vector_crs).to_epsg()
if vector_gdf.crs.to_epsg() != user_crs_input:
- sys.exit(
+ raise ValueError(
f"The CRS of vector file {vector_file_path} is {vector_gdf.crs}, and doesn't match the input CRS ({xy_crs}) defined for liquefaction triggering models"
)
except: # noqa: E722
@@ -362,11 +364,9 @@ def run(self, ln_im_data, eq_data, im_list, output_keys, additional_output_keys)
else:
additional_output.update({key: item})
else:
- sys.exit(
+ raise ValueError(
"At least one of 'PGA' and 'PGV' is missing in the selected intensity measures and the liquefaction trigging model 'ZhuEtal2017' can not be computed."
)
- # print(f"At least one of 'PGA' and 'PGV' is missing in the selected intensity measures and the liquefaction trigging model 'ZhuEtal2017' can not be computed."\
- # , file=sys.stderr)
# sys.stderr.write("test")
# sys.exit(-1)
return ln_im_data, eq_data, im_list, additional_output
@@ -579,7 +579,7 @@ def run(self, ln_im_data, eq_data, im_list, output_keys, additional_output_keys)
else:
additional_output.update({key: item})
else:
- sys.exit(
+ raise ValueError(
"'PGA'is missing in the selected intensity measures and the liquefaction trigging model 'Hazus2020' can not be computed."
)
return ln_im_data, eq_data, im_list, additional_output
@@ -921,7 +921,7 @@ def run(self, ln_im_data, eq_data, im_list): # noqa: D102
ln_im_data[scenario_id] = im_data_scen
im_list = im_list + output_keys
else:
- sys.exit(
+ raise ValueError(
"At least one of 'PGA' and 'PGV' is missing in the selected intensity measures and the liquefaction trigging model 'ZhuEtal2017' can not be computed."
)
return ln_im_data, eq_data, im_list
@@ -1090,7 +1090,7 @@ def run(self, ln_im_data, eq_data, im_list): # noqa: D102
ln_im_data[scenario_id] = im_data_scen
im_list = im_list + output_keys
else:
- sys.exit(
+ raise ValueError(
"At least one of 'liq_susc' and 'liq_prob' is missing in the selected intensity measures and the liquefaction trigging model 'ZhuEtal2017' can not be computed."
)
return ln_im_data, eq_data, im_list
diff --git a/modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA_25.4.py b/modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA.py
similarity index 93%
rename from modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA_25.4.py
rename to modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA.py
index 6fbbacae5..663693ac7 100644
--- a/modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA_25.4.py
+++ b/modules/performRegionalEventSimulation/regionalGroundMotion/testFetchOpenSHA.py
@@ -1,47 +1,42 @@
##############################################################
-# OpenSHA v25.4 Compatibility Test Suite
+# OpenSHA v26.1 Compatibility Test Suite
# Run with: python test_opensha25.py
##############################################################
import numpy as np
import socket
-import subprocess
-import importlib
-import sys
import psutil
import GlobalVariable
if 'stampede2' not in socket.gethostname():
- import GlobalVariable
-
if GlobalVariable.JVM_started is False:
GlobalVariable.JVM_started = True
- if importlib.util.find_spec('jpype') is None:
- subprocess.check_call([sys.executable, '-m', 'pip', 'install', 'JPype1'])
import jpype
import jpype.imports
- from jpype.types import *
memory_total = psutil.virtual_memory().total / (1024.0**3)
memory_request = int(memory_total * 0.75)
jpype.addClassPath('./lib/opensha-all.jar')
- jpype.startJVM(f'-Xmx{memory_request}G', convertStrings=False)
-
-from java.io import *
-from java.lang import *
-from java.lang.reflect import *
-from java.util import *
-from org.opensha.commons.data import *
-from org.opensha.commons.data.function import *
-from org.opensha.commons.data.siteData import *
-from org.opensha.commons.geo import *
-from org.opensha.commons.param import *
-from org.opensha.commons.param.constraint import *
-from org.opensha.commons.param.event import *
-from org.opensha.sha.calc import *
-from org.opensha.sha.earthquake import *
-from org.opensha.sha.earthquake.param import *
+ jpype.startJVM(
+ f'-Xmx{memory_request}G',
+ convertStrings=False,
+ jvmpath=GlobalVariable.find_compatible_jvm_path(),
+ )
+
+from java.util import ArrayList
+from org.opensha.commons.data import Site
+from org.opensha.commons.data.siteData import OrderedSiteDataProviderList
+from org.opensha.commons.geo import Location
+from org.opensha.commons.param.event import ParameterChangeWarningListener
+from org.opensha.sha.earthquake import EqkRupture
+from org.opensha.sha.earthquake.param import (
+ BPTAveragingTypeOptions,
+ BackgroundRupType,
+ IncludeBackgroundOption,
+ MagDependentAperiodicityOptions,
+ ProbabilityModelOptions,
+)
from org.opensha.sha.earthquake.rupForecastImpl.Frankel02 import (
Frankel02_AdjustableEqkRupForecast,
)
@@ -52,43 +47,49 @@
from org.opensha.sha.earthquake.rupForecastImpl.WGCEP_UCERF_2_Final.MeanUCERF2 import (
MeanUCERF2,
)
-from org.opensha.sha.faultSurface import *
-from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrection
+from org.opensha.sha.faultSurface import PointSurface
from org.opensha.sha.faultSurface.utils import PointSourceDistanceCorrections
-from org.opensha.sha.imr import *
-from org.opensha.sha.imr.attenRelImpl import *
-from org.opensha.sha.imr.attenRelImpl.ngaw2 import *
-from org.opensha.sha.imr.attenRelImpl.ngaw2.NGAW2_Wrappers import *
-from org.opensha.sha.imr.param.IntensityMeasureParams import *
-from org.opensha.sha.imr.param.OtherParams import *
-from org.opensha.sha.util import *
+from org.opensha.sha.gcim.calc import GcimCalculator
+from org.opensha.sha.gcim.imr.attenRelImpl import (
+ BommerEtAl_2009_AttenRel,
+ KS_2006_AttenRel,
+)
+from org.opensha.sha.gcim.imr.param.IntensityMeasureParams import (
+ Ds575_Param,
+ Ds595_Param,
+)
+from org.opensha.sha.imr.attenRelImpl import AfshariStewart_2016_AttenRel
+from org.opensha.sha.imr.attenRelImpl.ngaw2 import (
+ ASK_2014,
+ BSSA_2014,
+ CB_2014,
+ CY_2014,
+)
+from org.opensha.sha.imr.param.IntensityMeasureParams import PeriodParam, SA_Param
+from org.opensha.sha.imr.param.OtherParams import StdDevTypeParam
+from org.opensha.sha.util import SiteTranslator
+# UCERF3 lives under the upstream `scratch.*` namespace.
try:
from scratch.UCERF3.erf.mean import MeanUCERF3
except ModuleNotFoundError:
MeanUCERF3 = jpype.JClass('scratch.UCERF3.erf.mean.MeanUCERF3')
-from org.opensha.sha.gcim.calc import *
-from org.opensha.sha.gcim.imr.attenRelImpl import *
-from org.opensha.sha.gcim.imr.param.EqkRuptureParams import *
-from org.opensha.sha.gcim.imr.param.IntensityMeasureParams import *
-
-# Import the 10 functions from the main module
-from FetchOpenSHA_25_4 import (
- getERF,
+from FetchOpenSHA import (
CreateIMRInstance,
- get_site_prop,
get_IM,
get_PointSource_info_CY2014,
- horzDistanceFast,
- getPtSrcDistCorr,
+ get_site_prop,
get_site_vs30_from_opensha,
get_site_z1pt0_from_opensha,
get_site_z2pt5_from_opensha,
+ getERF,
+ getPtSrcDistCorr,
+ horzDistanceFast,
)
def run_tests():
- """Run all tests to verify OpenSHA v25.4 compatibility."""
+ """Run all tests to verify OpenSHA v26.1 compatibility."""
from jpype import JInt
results = {'passed': 0, 'failed': 0, 'skipped': 0}
@@ -312,7 +313,6 @@ def report(test_name, passed, msg=''):
report(f'{test_id} {name_str} via Wrapper', False, str(e))
try:
- from org.opensha.commons.param.event import ParameterChangeWarningListener
@jpype.JImplements(ParameterChangeWarningListener)
class DummyListener:
@@ -330,7 +330,6 @@ def parameterChangeWarning(self, event):
report('3.5 KS_2006 (gcim)', False, str(e))
try:
- from org.opensha.commons.param.event import ParameterChangeWarningListener
@jpype.JImplements(ParameterChangeWarningListener)
class DummyListener2:
@@ -730,7 +729,7 @@ def parameterChangeWarning(self, event):
# ==========================================================
print('\n' + '=' * 70)
- print('TEST GROUP 10: v25.4 PointSource Distance Correction API')
+ print('TEST GROUP 10: PointSource Distance Correction API')
print('=' * 70)
# ==========================================================
@@ -786,7 +785,7 @@ def parameterChangeWarning(self, event):
print('=' * 70)
if results['failed'] == 0:
- print('\nđ All tests passed! OpenSHA v25.4 integration is working.\n')
+ print('\nđ All tests passed! OpenSHA v26.1 integration is working.\n')
else:
print(f'\nâ ď¸ {results["failed"]} test(s) failed. See above.\n')
@@ -795,6 +794,6 @@ def parameterChangeWarning(self, event):
if __name__ == '__main__':
print('\n' + '=' * 70)
- print('OpenSHA v25.4 Compatibility Test Suite')
+ print('OpenSHA v26.1 Compatibility Test Suite')
print('=' * 70)
run_tests()
diff --git a/modules/systemPerformance/ResidualDemand/run_residual_demand.py b/modules/systemPerformance/ResidualDemand/run_residual_demand.py
index 765afaa82..fec90a8a5 100644
--- a/modules/systemPerformance/ResidualDemand/run_residual_demand.py
+++ b/modules/systemPerformance/ResidualDemand/run_residual_demand.py
@@ -100,7 +100,7 @@ def select_realizations_to_run(damage_input, run_dir):
else:
rlzs_requested.append(int(rlzs))
rlzs_requested = np.array(rlzs_requested)
- rlzs_in_available = np.in1d(rlzs_requested, rlzs_available) # noqa: NPY201
+ rlzs_in_available = np.isin(rlzs_requested, rlzs_available)
if rlzs_in_available.sum() != 0:
rlzs_to_run = rlzs_requested[np.where(rlzs_in_available)[0]]
else:
diff --git a/modules/systemPerformance/ResidualDemand/transportation.py b/modules/systemPerformance/ResidualDemand/transportation.py
index 809354e61..6cceef009 100644
--- a/modules/systemPerformance/ResidualDemand/transportation.py
+++ b/modules/systemPerformance/ResidualDemand/transportation.py
@@ -44,6 +44,18 @@
from __future__ import annotations # noqa: I001
+import ctypes
+import sys
+
+# Programmatic fix for Apple Silicon macOS flat namespace errors
+if sys.platform == "darwin":
+ try:
+ # Force-load libomp globally so its symbols are visible to dlopen()
+ ctypes.CDLL('/opt/homebrew/opt/libomp/lib/libomp.dylib', mode=ctypes.RTLD_GLOBAL)
+ except OSError:
+ # Fallback for systems where Homebrew is not installed or configured differently
+ pass
+
import gc
import os
import json
@@ -63,7 +75,8 @@
import geopandas as gpd
import numpy as np
-import pandana.network as pdna
+#import pandana.network as pdna
+import pandarm as pdna
import pandas as pd
from sklearn.feature_extraction.text import TfidfVectorizer
from scipy.spatial.distance import cdist