Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

abs - surrogate aided Bayesian calibration algorithm implementation #331

Merged
merged 4 commits into from
Sep 30, 2024
Merged
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 9 additions & 12 deletions modules/createEVENT/shakeMapEvent/shakeMapEvent.py
Original file line number Diff line number Diff line change
@@ -45,10 +45,9 @@


def create_shakemap_event(eventDirectory, eventPath, IMTypes): # noqa: D103
IMTypesList = eval(IMTypes)

IMTypesList = eval(IMTypes)

print("Creating shakemap event")
print('Creating shakemap event')

xml_file_path = Path(eventDirectory) / eventPath / 'grid.xml'

@@ -69,7 +68,7 @@ def create_shakemap_event(eventDirectory, eventPath, IMTypes): # noqa: D103
lon, lat = float(values[0]), float(values[1])
point = Point(lon, lat)
points.append(point)

# Store only the specified attributes
attr = {}
attribute_mapping = {
@@ -78,37 +77,35 @@ def create_shakemap_event(eventDirectory, eventPath, IMTypes): # noqa: D103
'MMI': 4,
'PSA03': 5,
'PSA10': 6,
'PSA30': 7
'PSA30': 7,
}

for im_type in IMTypesList:
if im_type in attribute_mapping:
attr[im_type] = float(values[attribute_mapping[im_type]])

attributes.append(attr)

# Create GeoDataFrame
gdf = gpd.GeoDataFrame(attributes, geometry=points, crs="EPSG:4326")
gdf = gpd.GeoDataFrame(attributes, geometry=points, crs='EPSG:4326')

# Display the first few rows
print("Saving shakemap to gpkg")
print('Saving shakemap to gpkg')

# Save as a GeoPackage file
gdf_path = Path(eventDirectory) / 'EventGrid.gpkg'
gdf.to_file(gdf_path, driver="GPKG")
gdf.to_file(gdf_path, driver='GPKG')

return


if __name__ == '__main__':

parser = argparse.ArgumentParser()
parser.add_argument('--input', help='Input file')
parser.add_argument('--Directory', help='Directory path')
parser.add_argument('--EventPath', help='Event path')
parser.add_argument('--IntensityMeasureType', help='types of intensity measures')


args = parser.parse_args()

create_shakemap_event(args.Directory, args.EventPath, args.IntensityMeasureType)
Original file line number Diff line number Diff line change
@@ -43,39 +43,42 @@

from RasterEvent import create_event as create_raster_event


def is_raster_file(filename):
# Define a set of common raster file extensions
raster_extensions = {'.jpg', '.jpeg', '.png', '.bmp', '.gif', '.tiff', '.tif'}

# Create a Path object from the filename
file_path = Path(filename)

# Extract the file extension and check if it is in the set of raster extensions
return file_path.suffix.lower() in raster_extensions


def is_xml_file(filename):
# Check if the file has an .xml extension
if not filename.lower().endswith('.xml'):
return False

# Try to parse the file as XML
try:
ET.parse(filename)
return True
except ET.ParseError:
return False

def create_event(asset_file : str, event_grid_file: str): # noqa: C901, N803, D103

def create_event(asset_file: str, event_grid_file: str): # noqa: C901, N803, D103
if is_raster_file(event_grid_file):
return create_raster_event(asset_file, event_grid_file)
elif is_xml_file(event_grid_file):
# Here you would call a function to handle XML files
# For now, we'll just raise a NotImplementedError
raise NotImplementedError("XML file handling is not yet implemented.")
raise NotImplementedError('XML file handling is not yet implemented.')
else:
raise ValueError(f"{event_grid_file} is not a raster. Only rasters are currently supported.")

raise ValueError(
f'{event_grid_file} is not a raster. Only rasters are currently supported.'
)


if __name__ == '__main__':
@@ -84,6 +87,4 @@ def create_event(asset_file : str, event_grid_file: str): # noqa: C901, N803, D
parser.add_argument('--filenameEVENTgrid')
args = parser.parse_args()

create_event(
args.assetFile, args.filenameEVENTgrid
)
create_event(args.assetFile, args.filenameEVENTgrid)
56 changes: 24 additions & 32 deletions modules/performRegionalMapping/GISSpecifiedEvents/RasterEvent.py
Original file line number Diff line number Diff line change
@@ -51,16 +51,15 @@ def sample_raster_at_latlon(src, lat, lon):

# Ensure the indices are within the bounds of the raster
if row < 0 or row >= src.height or col < 0 or col >= src.width:
raise IndexError("Transformed coordinates are out of raster bounds")
raise IndexError('Transformed coordinates are out of raster bounds')

# Read the raster value at the given row and column
raster_value = src.read(1)[row, col]

return raster_value

def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103


def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103
# read the event grid data file
event_grid_path = Path(event_grid_file).resolve()
event_dir = event_grid_path.parent
@@ -82,22 +81,21 @@ def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103
asset_dict = json.load(f)

data_final = [
['GP_file','Latitude','Longitude'],
]
['GP_file', 'Latitude', 'Longitude'],
]

# Iterate through each asset
for asset in asset_dict:
asset_id = asset['id']
asset_file_path = asset['file']

# Load the corresponding file for each asset
with open(asset_file_path, encoding='utf-8') as asset_file :

# Load the corresponding file for each asset
with open(asset_file_path, encoding='utf-8') as asset_file:
# Load the asset data
asset_data = json.load(asset_file)

im_tag = asset_data['RegionalEvent']['intensityMeasures'][0]

# Extract the latitude and longitude
lat = float(asset_data['GeneralInformation']['location']['latitude'])
lon = float(asset_data['GeneralInformation']['location']['longitude'])
@@ -107,34 +105,33 @@ def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103

# Check if the transformed coordinates are within the raster bounds
bounds = src.bounds
if (bounds.left <= lon_transformed <= bounds.right and
bounds.bottom <= lat_transformed <= bounds.top):
if (
bounds.left <= lon_transformed <= bounds.right
and bounds.bottom <= lat_transformed <= bounds.top
):
try:
val = sample_raster_at_latlon(src=src,
lat=lat_transformed,
lon=lon_transformed)

data = [
[im_tag],
[val]
]

val = sample_raster_at_latlon(
src=src, lat=lat_transformed, lon=lon_transformed
)

data = [[im_tag], [val]]

# Save the simcenter file name
file_name = f'Site_{asset_id}.csvx{0}x{int(asset_id):05d}'

data_final.append([file_name,lat,lon])
data_final.append([file_name, lat, lon])

csv_save_path = event_dir / f'Site_{asset_id}.csv'
with open(csv_save_path, 'w', newline='') as file:
# Create a CSV writer object
writer = csv.writer(file)

# Write the data to the CSV file
writer.writerows(data)

# prepare a dictionary of events
event_list_json = [[file_name, 1.0]]

asset_data['Events'] = [{}]
asset_data['Events'][0] = {
'EventFolderPath': str(event_dir),
@@ -145,12 +142,10 @@ def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103
with open(asset_file_path, 'w', encoding='utf-8') as f: # noqa: PTH123
json.dump(asset_data, f, indent=2)


except IndexError as e:
print(f"Error for asset ID {asset_id}: {e}")
print(f'Error for asset ID {asset_id}: {e}')
else:
print(f"Asset ID: {asset_id} is outside the raster bounds")

print(f'Asset ID: {asset_id} is outside the raster bounds')

# # save the event dictionary to the BIM
# asset_data['Events'] = [{}]
@@ -165,13 +160,12 @@ def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103
# with open(asset_file, 'w', encoding='utf-8') as f: # noqa: PTH123
# json.dump(asset_data, f, indent=2)


# Save the final event grid
csv_save_path = event_dir / 'EventGrid.csv'
with open(csv_save_path, 'w', newline='') as file:
# Create a CSV writer object
writer = csv.writer(file)

# Write the data to the CSV file
writer.writerows(data_final)

@@ -185,6 +179,4 @@ def create_event(asset_file, event_grid_file): # noqa: C901, N803, D103
parser.add_argument('--filenameEVENTgrid')
args = parser.parse_args()

create_event(
args.assetFile, args.filenameEVENTgrid
)
create_event(args.assetFile, args.filenameEVENTgrid)
43 changes: 29 additions & 14 deletions modules/performRegionalMapping/NearestNeighborEvents/NNE.py
Original file line number Diff line number Diff line change
@@ -48,6 +48,7 @@
from sklearn.neighbors import NearestNeighbors
import geopandas as gpd


def find_neighbors( # noqa: C901, D103
asset_file,
event_grid_file,
@@ -89,7 +90,7 @@ def find_neighbors( # noqa: C901, D103
if file_extension == '.csv':
# Existing code for CSV files
grid_df = pd.read_csv(event_dir / event_grid_file, header=0)

# store the locations of the grid points in X
lat_E = grid_df['Latitude'] # noqa: N806
lon_E = grid_df['Longitude'] # noqa: N806
@@ -103,15 +104,15 @@ def find_neighbors( # noqa: C901, D103
else:
# Else assume GIS files - works will all gis files that geopandas supports
gdf = gpd.read_file(event_dir / event_grid_file)

# Ensure the GIS file is in a geographic coordinate system
if not gdf.crs.is_geographic:
gdf = gdf.to_crs(epsg=4326) # Convert to WGS84

# Extract coordinates from the geometry
gdf['Longitude'] = gdf.geometry.x
gdf['Latitude'] = gdf.geometry.y

# store the locations of the grid points in X
lat_E = gdf['Latitude'] # noqa: N806
lon_E = gdf['Longitude'] # noqa: N806
@@ -121,7 +122,7 @@ def find_neighbors( # noqa: C901, D103
grid_extra_keys = list(
gdf.drop(['geometry', 'Longitude', 'Latitude'], axis=1).columns
)

# Convert GeoDataFrame to regular DataFrame for consistency with the rest of the code
grid_df = pd.DataFrame(gdf.drop(columns='geometry'))

@@ -131,7 +132,9 @@ def find_neighbors( # noqa: C901, D103
else:
neighbors_to_get = neighbors

nbrs = NearestNeighbors(n_neighbors=neighbors_to_get, algorithm='ball_tree').fit(X)
nbrs = NearestNeighbors(n_neighbors=neighbors_to_get, algorithm='ball_tree').fit(
X
)

# load the building data file
with open(asset_file, encoding='utf-8') as f: # noqa: PTH123
@@ -307,34 +310,46 @@ def find_neighbors( # noqa: C901, D103
] * e

scale_list = np.ones(len(event_list))
else :
else:
event_list = []
scale_list = []
event_type = 'intensityMeasure'

# Determine event_count (number of IMs per grid point)
im_columns = [col for col in grid_df.columns if col not in ['geometry', 'Longitude', 'Latitude']]
im_columns = [
col
for col in grid_df.columns
if col not in ['geometry', 'Longitude', 'Latitude']
]
event_count = len(im_columns)

# for each neighbor
for sample_j, nbr in enumerate(nbr_samples):

# make sure we resample events if samples > event_count
event_j = sample_j % event_count

# get the index of the nth neighbor
nbr_index = ind_list[nbr]

# For GIS files, create a new CSV file
csv_filename = f'Site_{sample_j}.csv'

csv_path = event_dir / csv_filename

# Create a CSV file with data from the GIS file
# Use actual data from the GIS file if available, otherwise use dummy data
im_columns = [col for col in grid_df.columns if col not in ['geometry', 'Longitude', 'Latitude']]

im_data = pd.DataFrame({col: [grid_df.iloc[nbr_index][col]] * event_count for col in im_columns})
im_columns = [
col
for col in grid_df.columns
if col not in ['geometry', 'Longitude', 'Latitude']
]

im_data = pd.DataFrame(
{
col: [grid_df.iloc[nbr_index][col]] * event_count
for col in im_columns
}
)

im_data.to_csv(csv_path, index=False)
# save the collection file name and the IM row id
1 change: 1 addition & 0 deletions modules/performUQ/common/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1 @@
# noqa: D104
Loading
Loading