official-apps/SeismicHazardForecasting
11
0
Fork 2
Code Issues Pull Requests 1 Activity

Update src/seismic_hazard_forecasting.py

Browse Source 
change IMT plotting method
This commit is contained in:
ftong
2026-07-08 16:12:39 +02:00
parent ab19294be5
commit a56f18dde0
1 changed files with 70 additions and 171 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/seismic_hazard_forecasting.py
+70 -171
View File
@@ -709,74 +709,43 @@ verbose: {verbose}")
m_cdf = get_cdf(m_pdf) m_cdf = get_cdf(m_pdf)
fxy = xy_kde[0] centroids_utm = grid_gdf_utm.geometry.centroid.values #extract the centroid of each cell
logger.debug(f"Normalization check; sum of all f(x,y) values = {np.sum(fxy)}") num_points = len(grid_gdf_utm)
xx, yy = np.meshgrid(x_range, y_range) # grid points distances = np.array([shapely.distance(centroids_utm[i], centroids_utm) for i in range(num_points)]) #compute distance between every grid point
grid_gdf_latlon['distance_matrix'] = [distances[i] for i in range(num_points)] #store the distance matrix in the GDF
# set every grid point to be a receiver
grid_shape = xx.shape
x_rx = xx.flatten()
y_rx = yy.flatten()
num_points = x_rx.size
distances = np.zeros(shape=(num_points, grid_shape[0], grid_shape[1]))
# compute distance matrix for each receiver
#distances = np.zeros(shape=(nx * ny, nx, ny))
rx_lat = np.zeros(nx * ny)
rx_lon = np.zeros(nx * ny)
for i in range(num_points):
# Compute the squared distances directly using NumPy's vectorized operations
squared_distances = (xx - x_rx[i]) ** 2 + (yy - y_rx[i]) ** 2
distances[i] = np.sqrt(squared_distances)
# create context object for receiver and append to list
rx_lat[i], rx_lon[i] = utm.to_latlon(x_rx[i], y_rx[i], utm_zone_number,
utm_zone_letter) # get receiver location as lat,lon
# convert distances from m to km because openquake ground motion models take input distances in kilometres
#distances = distances/1000.0
# compute ground motion only at grid points that have minimum probability density of thresh_fxy
if exclude_low_fxy:
indices = list(np.where(fxy.flatten() > thresh_fxy)[0])
else:
indices = np.arange(num_points)
# Select only cells of the grid that are inside the AOI
if use_AOI: if use_AOI:
# Filter out receivers outside the AOI; Find indices where values are OUTSIDE the AOI centroids_latlon = grid_gdf_latlon.geometry.centroid
indices_outside_x = np.where((x_rx < x_AOI[0]) | (x_rx > x_AOI[1]))[0]
indices_outside_y = np.where((y_rx < y_AOI[0]) | (y_rx > y_AOI[1]))[0]
indices_outside_AOI = np.unique(np.concatenate((indices_outside_x, indices_outside_y)))
indices_filtered = np.setdiff1d(indices, indices_outside_AOI)
# AOI grid extent
AOI_rx_x = x_rx[indices_filtered]
AOI_rx_y = y_rx[indices_filtered]
AOI_rx_lat, AOI_rx_lon = utm.to_latlon(AOI_rx_x, AOI_rx_y, utm_zone_number, utm_zone_letter) # Mark grid cells that are within the AOI using vectorized boundary checks
grid_gdf_latlon['AOI'] = (
logger.debug(f"Receiver UTM X range: {AOI_rx_x.min()} to {AOI_rx_x.max()}") (centroids_latlon.x >= AOI_lon[0]) & (centroids_latlon.x <= AOI_lon[1]) &
logger.debug(f"Receiver UTM Y range: {AOI_rx_y.min()} to {AOI_rx_y.max()}") (centroids_latlon.y >= AOI_lat[0]) & (centroids_latlon.y <= AOI_lat[1])
logger.debug(f"Receiver lat range: {AOI_rx_lat.min()} to {AOI_rx_lat.max()}") )
logger.debug(f"Receiver lon range: {AOI_rx_lon.min()} to {AOI_rx_lon.max()}")
distances_sel = grid_gdf_latlon.loc[grid_gdf_latlon['AOI']]['distance_matrix'].to_numpy()
centroids_sel = grid_gdf_latlon.loc[grid_gdf_latlon['AOI']].centroid
else: else:
indices_filtered = indices grid_gdf_latlon['AOI']=True #set entire grid to be the area of interest
distances_sel = grid_gdf_latlon['distance_matrix'].to_numpy()
fr = fxy.flatten() centroids_sel = grid_gdf_latlon.centroid
loc_pdf = grid_gdf_latlon['location_PDF'].to_numpy() # extract the previously created location PDF from the GDF
# convert distances from m to km because openquake ground motion models take input distances in kilometres
#distances_sel = distances_sel/1000.0
# For each receiver compute estimated ground motion values # For each receiver compute estimated ground motion values
logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...") logger.info(f"Estimating ground motion intensity at {len(distances_sel)} grid points...")
start = timer()
use_pp = True use_pp = True
start = timer()
if use_pp: # use dask parallel computing if use_pp: # use dask parallel computing
mp.set_start_method("fork", force=True) mp.set_start_method("fork", force=True)
iter = indices_filtered iter = range(0,len(distances_sel))
iml_grid_raw = [] # raw ground motion grids iml_grid_raw = [] # raw ground motion grids
for imt in products: for imt in products:
logger.info(f"Estimating {imt}") logger.info(f"Estimating {imt}")
@@ -786,7 +755,7 @@ verbose: {verbose}")
else: else:
IMT_max = 2.0 # search interval max for acceleration (g) IMT_max = 2.0 # search interval max for acceleration (g)
imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, imls = [dask.delayed(compute_IMT_exceedance)(centroids_sel.iloc[i].y, centroids_sel.iloc[i].x, distances_sel[i].flatten(), loc_pdf, p, lambdas,
forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model, forecast_len, lambdas_perc, m_range, m_pdf, m_cdf, model,
log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i, log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=IMT_max, rx_label=i,
rtol=0.1, use_cython=True) for i in iter] rtol=0.1, use_cython=True) for i in iter]
@@ -795,9 +764,9 @@ verbose: {verbose}")
else: else:
iml_grid_raw = [] iml_grid_raw = []
iter = indices_filtered iter = range(0,len(distances_sel))
for imt in products: for imt in products:
if imt == "PGV": if imt == "PGV":
IMT_max = 200 # search interval max for velocity (cm/s) IMT_max = 200 # search interval max for velocity (cm/s)
else: else:
@@ -805,7 +774,7 @@ verbose: {verbose}")
iml = [] iml = []
for i in iter: for i in iter:
iml_i = compute_IMT_exceedance(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, p, lambdas, forecast_len, iml_i = compute_IMT_exceedance(centroids_sel.iloc[i].y, centroids_sel.iloc[i].x, distances_sel[i].flatten(), loc_pdf, p, lambdas, forecast_len,
lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0, lambdas_perc, m_range, m_pdf, m_cdf, model, imt=imt, IMT_min = 0.0,
IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True) IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True)
iml.append(iml_i) iml.append(iml_i)
@@ -815,127 +784,57 @@ verbose: {verbose}")
end = timer() end = timer()
logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds") logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
logger.debug(f"IMT values: {iml_grid_raw[0]}")
if np.isnan(iml_grid_raw).all(): if np.isnan(iml_grid_raw).all():
msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model." msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model."
logger.error(msg) logger.error(msg)
raise Exception(msg) raise Exception(msg)
# create list of one empty list for each imt
iml_grid = [[] for _ in range(len(products))] # final ground motion grids
iml_grid_prep = iml_grid.copy() # temp ground motion grids
#if use_AOI or exclude_low_fxy:
for j in range(0, len(products)):
# Reassemble the grid cleanly using the original shape
# Initialize a flat array filled entirely with NaNs
iml_grid_flat = np.full(num_points, np.nan, dtype=np.float64)
# Assign the computed values to their exact original 1D index positions
iml_grid_flat[indices_filtered] = iml_grid_raw[j]
# Reshape back using the exact shape of your original xx/yy grids
iml_grid_prep[j] = iml_grid_flat.reshape(grid_shape)
#for i in indices:
# if i in indices_filtered:
# for j in range(0, len(products)):
# iml_grid_prep[j].append(iml_grid_raw[j].pop(0))
# else:
# list(map(lambda lst: lst.append(np.nan),
# iml_grid_prep)) # use np.nan to indicate grid point excluded
#elif exclude_low_fxy:
# for i in range(0, len(distances)):
# if i in indices:
# for j in range(0, len(products)):
# iml_grid_prep[j].append(iml_grid_raw[j].pop(0))
# else:
# list(map(lambda lst: lst.append(np.nan),
# iml_grid_prep)) # use np.nan to indicate grid point excluded
#else:
# iml_grid_prep = iml_grid_raw
if use_AOI:
# Update grid extents
grid_x_min = AOI_rx_x.min()
grid_x_max = AOI_rx_x.max()
grid_y_min = AOI_rx_y.min()
grid_y_max = AOI_rx_y.max()
grid_lat_min = AOI_rx_lat.min()
grid_lat_max = AOI_rx_lat.max()
grid_lon_min = AOI_rx_lon.min()
grid_lon_max = AOI_rx_lon.max()
for j in range(0, len(products)): for j, imt in enumerate(products): #generate image overlay for each IMT product
logger.debug(f"{products[j]} values: {iml_grid_raw[j]}")
grid_gdf_latlon.loc[grid_gdf_latlon['AOI'], imt] = iml_grid_raw[j] # insert computed imt grid into GDF
if use_AOI: grid_gdf_latlon_clean = grid_gdf_latlon.dropna(subset=[imt]) # remove null values from grid
# trim grid to remove all nan values
x_plot = grid_gdf_latlon_clean.geometry.centroid.x.values
# Create a boolean mask of non-NaN values y_plot = grid_gdf_latlon_clean.geometry.centroid.y.values
# ~np.isnan() returns True for values and False for NaNs z_plot = grid_gdf_latlon_clean[imt].values
nan_mask = ~np.isnan(iml_grid_prep[j])
vmin = np.nanmin(z_plot)
# Identify valid rows and columns vmax = np.nanmax(z_plot)
# .any(axis=1) checks each row; .any(axis=0) checks each column
row_mask = nan_mask.any(axis=1)
col_mask = nan_mask.any(axis=0)
# Extract the sub-array ---
# np.ix_ creates an open mesh from multiple boolean arrays so they can be broadcast together
iml_grid_prep[j] = iml_grid_prep[j][np.ix_(row_mask, col_mask)]
vmin = np.nanmin(iml_grid_prep[j]) # Generate Image Overlay
vmax = np.nanmax(iml_grid_prep[j]) fig, ax = plt.subplots()
contour = ax.tricontourf(
#iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(dtype=np.float64) # this reduces values to 8 decimal places x_plot,
#iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes y_plot,
z_plot,
# upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1 levels=200, #linear scale
#if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1: cmap="YlOrRd",
# up_factor = 1 )
# iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), ax.set_aspect('equal') # keep geographic coordinates from stretching
# mode='reflect', anti_aliasing=False) ax.set_axis_off()
# trim_thresh = vmin fig.patch.set_visible(False); ax.patch.set_visible(False)
# iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
#else:
#iml_grid_hd = iml_grid_tmp
#iml_grid_hd[iml_grid_hd == 0.0] = np.nan # change zeroes back to nan
iml_grid_hd = iml_grid_prep[j]
#vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
vmax_hd = np.nanmax(iml_grid_hd)
# generate image overlay
north, south = grid_lat_max, grid_lat_min # Latitude range
east, west = grid_lon_max, grid_lon_min # Longitude range
bounds = [[south, west], [north, east]]
map_center = [np.mean([north, south]), np.mean([east, west])]
# Create an image from the grid
cmap_name = 'YlOrRd'
cmap = plt.get_cmap(cmap_name)
#fig, ax = plt.subplots(figsize=(6, 6))
fig, ax = plt.subplots(layout=None)
ax.margins(0) # clear any data margins
ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear', aspect='auto')
ax.axis('off')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
# Save the figure
fig.canvas.draw()
overlay_filename = f"overlay_{j}.svg" overlay_filename = f"overlay_{j}.svg"
plt.savefig(overlay_filename, pad_inches=0, transparent=True) plt.savefig(overlay_filename, pad_inches=0, bbox_inches="tight", transparent=True)
plt.close(fig) plt.close(fig)
# set image map extent in geographic coordinates
if use_AOI:
north = AOI_lat[1]
south = AOI_lat[0]
east = AOI_lon[1]
west = AOI_lon[0]
else:
north = grid_gdf_latlon.total_bounds[3]
south = grid_gdf_latlon.total_bounds[1]
east = grid_gdf_latlon.total_bounds[2]
west = grid_gdf_latlon.total_bounds[0]
# Embed geographic bounding box into the SVG # Embed geographic bounding box into the SVG
map_bounds = dict(zip(("south", "west", "north", "east"), map_bounds = dict(zip(("south", "west", "north", "east"),
map(float, (grid_lat_min, grid_lon_min, grid_lat_max, grid_lon_max)))) map(float, (south, west, north, east))))
tree = ET.parse(overlay_filename) tree = ET.parse(overlay_filename)
tree.getroot().set("data-map-bounds", json.dumps(map_bounds)) tree.getroot().set("data-map-bounds", json.dumps(map_bounds))
tree.write(overlay_filename, encoding="utf-8", xml_declaration=True) tree.write(overlay_filename, encoding="utf-8", xml_declaration=True)
@@ -949,14 +848,14 @@ verbose: {verbose}")
gradient = np.vstack((gradient, gradient)).T gradient = np.vstack((gradient, gradient)).T
gradient = np.tile(gradient, (1, width)) gradient = np.tile(gradient, (1, width))
colorbar_title = products[j] colorbar_title = imt
if "PGA" in colorbar_title or "SA" in colorbar_title: if "PGA" in colorbar_title or "SA" in colorbar_title:
colorbar_title = colorbar_title + " (g)" colorbar_title = colorbar_title + " (g)"
fig, ax = plt.subplots(figsize=((width + 40) / 100.0, (height + 20) / 100.0), fig, ax = plt.subplots(figsize=((width + 40) / 100.0, (height + 20) / 100.0),
dpi=100) # Increase fig size for labels dpi=100) # Increase fig size for labels
ax.imshow(gradient, aspect='auto', cmap=cmap.reversed(), ax.imshow(gradient, aspect='auto', cmap=cmap.reversed(),
extent=[0, 1, vmin, vmax_hd]) # Note: extent order is different for vertical extent=[0, 1, vmin, vmax]) # Note: extent order is different for vertical
ax.set_xticks([]) # Remove x-ticks for vertical colorbar ax.set_xticks([]) # Remove x-ticks for vertical colorbar
num_ticks = 11 # Show more ticks num_ticks = 11 # Show more ticks
tick_positions = np.linspace(vmin, vmax_hd, num_ticks) tick_positions = np.linspace(vmin, vmax_hd, num_ticks)