diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 6b1968d..95af87a 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -709,74 +709,43 @@ verbose: {verbose}") m_cdf = get_cdf(m_pdf) - fxy = xy_kde[0] - logger.debug(f"Normalization check; sum of all f(x,y) values = {np.sum(fxy)}") + centroids_utm = grid_gdf_utm.geometry.centroid.values #extract the centroid of each cell + num_points = len(grid_gdf_utm) - xx, yy = np.meshgrid(x_range, y_range) # grid points - - # 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) + 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 + # Select only cells of the grid that are inside the AOI if use_AOI: - # Filter out receivers outside the AOI; Find indices where values are OUTSIDE the AOI - 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] + centroids_latlon = grid_gdf_latlon.geometry.centroid - AOI_rx_lat, AOI_rx_lon = utm.to_latlon(AOI_rx_x, AOI_rx_y, utm_zone_number, utm_zone_letter) - - logger.debug(f"Receiver UTM X range: {AOI_rx_x.min()} to {AOI_rx_x.max()}") - logger.debug(f"Receiver UTM Y range: {AOI_rx_y.min()} to {AOI_rx_y.max()}") - 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()}") + # Mark grid cells that are within the AOI using vectorized boundary checks + grid_gdf_latlon['AOI'] = ( + (centroids_latlon.x >= AOI_lon[0]) & (centroids_latlon.x <= AOI_lon[1]) & + (centroids_latlon.y >= AOI_lat[0]) & (centroids_latlon.y <= AOI_lat[1]) + ) + + 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: - indices_filtered = indices - - fr = fxy.flatten() + grid_gdf_latlon['AOI']=True #set entire grid to be the area of interest + distances_sel = grid_gdf_latlon['distance_matrix'].to_numpy() + 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 - logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...") - - start = timer() + logger.info(f"Estimating ground motion intensity at {len(distances_sel)} grid points...") use_pp = True - + + start = timer() if use_pp: # use dask parallel computing mp.set_start_method("fork", force=True) - iter = indices_filtered + iter = range(0,len(distances_sel)) iml_grid_raw = [] # raw ground motion grids for imt in products: logger.info(f"Estimating {imt}") @@ -786,7 +755,7 @@ verbose: {verbose}") else: 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, 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] @@ -795,9 +764,9 @@ verbose: {verbose}") else: iml_grid_raw = [] - iter = indices_filtered + iter = range(0,len(distances_sel)) + for imt in products: - if imt == "PGV": IMT_max = 200 # search interval max for velocity (cm/s) else: @@ -805,7 +774,7 @@ verbose: {verbose}") iml = [] 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, IMT_max = IMT_max, rx_label = i, rtol = 0.1, use_cython=True) iml.append(iml_i) @@ -815,127 +784,57 @@ verbose: {verbose}") end = timer() 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(): msg = "No valid ground motion intensity measures were forecasted. Try a different ground motion model." logger.error(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: - # trim grid to remove all nan values - - # Create a boolean mask of non-NaN values - # ~np.isnan() returns True for values and False for NaNs - nan_mask = ~np.isnan(iml_grid_prep[j]) - - # Identify valid rows and columns - # .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)] + grid_gdf_latlon_clean = grid_gdf_latlon.dropna(subset=[imt]) # remove null values from grid + + x_plot = grid_gdf_latlon_clean.geometry.centroid.x.values + y_plot = grid_gdf_latlon_clean.geometry.centroid.y.values + z_plot = grid_gdf_latlon_clean[imt].values + + vmin = np.nanmin(z_plot) + vmax = np.nanmax(z_plot) - vmin = np.nanmin(iml_grid_prep[j]) - vmax = np.nanmax(iml_grid_prep[j]) - - #iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(dtype=np.float64) # this reduces values to 8 decimal places - #iml_grid_tmp = np.nan_to_num(iml_grid[j]) # change nans to zeroes - - # upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1 - #if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1: - # up_factor = 1 - # iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)), - # mode='reflect', anti_aliasing=False) - # trim_thresh = vmin - # 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() + # Generate Image Overlay + fig, ax = plt.subplots() + contour = ax.tricontourf( + x_plot, + y_plot, + z_plot, + levels=200, #linear scale + cmap="YlOrRd", + ) + ax.set_aspect('equal') # keep geographic coordinates from stretching + ax.set_axis_off() + fig.patch.set_visible(False); ax.patch.set_visible(False) 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) + # 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 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.getroot().set("data-map-bounds", json.dumps(map_bounds)) tree.write(overlay_filename, encoding="utf-8", xml_declaration=True) @@ -949,14 +848,14 @@ verbose: {verbose}") gradient = np.vstack((gradient, gradient)).T gradient = np.tile(gradient, (1, width)) - colorbar_title = products[j] + colorbar_title = imt if "PGA" in colorbar_title or "SA" in colorbar_title: colorbar_title = colorbar_title + " (g)" fig, ax = plt.subplots(figsize=((width + 40) / 100.0, (height + 20) / 100.0), dpi=100) # Increase fig size for labels 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 num_ticks = 11 # Show more ticks tick_positions = np.linspace(vmin, vmax_hd, num_ticks)