From ab19294be5bc48448cbb7af04b14f39ab9acac83 Mon Sep 17 00:00:00 2001 From: ftong <95+ftong@noreply.example.org> Date: Tue, 7 Jul 2026 16:40:00 +0200 Subject: [PATCH] Update src/seismic_hazard_forecasting.py use geopandas for coordinate conversion before Location PDF estimation --- src/seismic_hazard_forecasting.py | 119 +++++++++++++++++++------------------- 1 file changed, 58 insertions(+), 61 deletions(-) diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py index 5e2b674..6b1968d 100644 --- a/src/seismic_hazard_forecasting.py +++ b/src/seismic_hazard_forecasting.py @@ -12,10 +12,6 @@ ncp_choice = 'default' tcp_max = 5 torder_min = 0 torder_max = 1 -#AOI_lat = np.array([51.48, 51.54]) -#AOI_lon = np.array([16.15, 16.24]) -#AOI_lat = np.array([None, None]) -#AOI_lon = np.array([None, None]) def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries, bin_dur, unit, multiplicator, @@ -164,7 +160,7 @@ def apply_beast(act_rate): tseg_minlength=2, mcmc_chains=10, mcmc_thin=mcmc_th, mcmc_seed=10) -# User-driven ncp selection + # User-driven ncp selection if ncp_choice == 'median': ncp = beast_result.trend.ncp_median if np.isnan(ncp) or ncp == 0: @@ -188,7 +184,7 @@ def apply_beast(act_rate): ncp = int(ncp) # Filter NaNs first — BEAST fills unused cp slots with nan - valid_cps = beast_result.trend.cp[~np.isnan(beast_result.trend.cp)] + # valid_cps = beast_result.trend.cp[~np.isnan(beast_result.trend.cp)] cps = beast_result.trend.cp[:ncp] # Discard mirrored zone changepoints and correct indices @@ -330,6 +326,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max import xml.etree.ElementTree as ET import json import multiprocessing as mp + import geopandas as gpd + import shapely logger = getDefaultLogger('igfash') @@ -362,7 +360,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max xy_select: {xy_select}\n grid_dim: {grid_dim}\n xy_win_method: {xy_win_method}\n rate_select: {rate_select}\n time_win_duration: {time_win_duration}\n \ forecast_select: {forecast_select}\n custom_rate: {custom_rate}\n forecast_len: {forecast_len}\n time_unit: {time_unit}\n model: {model}\n products: {products_string}\n \ verbose: {verbose}") - + logger.debug(f"Area of interest selected by user - Latitude: {AOI_lat}, Longitude: {AOI_lon}") # print key package version numbers logger.debug(f"Python version {sys.version}") logger.debug(f"Numpy version {version('numpy')}") @@ -471,69 +469,67 @@ verbose: {verbose}") time, mag, lat, lon, depth = read_mat_cat(catalog_file) - # convert to UTM - u = utm.from_latlon(lat, lon) - x = u[0] - y = u[1] - utm_zone_number = u[2] - utm_zone_letter = u[3] - logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}") + # Create a GeoDataFrame for the catalog data (initially WGS84 EPSG:4326) + catalog_gdf = gpd.GeoDataFrame( + {'depth': depth, 'time': time}, + geometry=gpd.points_from_xy(lon, lat), + crs="EPSG:4326" + ) + utm_crs = catalog_gdf.estimate_utm_crs() # Find the UTM EPSG code for this location + catalog_gdf_utm = catalog_gdf.to_crs(utm_crs) # Convert the catalog entirely to UTM meters + logger.debug(f"Latitude / Longitude event coordinates converted to UTM zone {utm_crs}") + + # Extract event coordinates directly from the vector geometry + x = catalog_gdf_utm.geometry.x.values + y = catalog_gdf_utm.geometry.y.values + + # Handle Area of Interest (AOI) if (None not in AOI_lat) and (None not in AOI_lon): use_AOI = True - logger.info(f"Area of Interest (AOI) selected with latitutde range {AOI_lat} and longitude range {AOI_lon}") - #convert AOI to UTM - u_AOI = utm.from_latlon(AOI_lat, AOI_lon) - x_AOI = u_AOI[0] - y_AOI = u_AOI[1] - - # make sure grid contains the user's AOI - x_min = np.concatenate((x, x_AOI)).min() - y_min = np.concatenate((y, y_AOI)).min() - x_max = np.concatenate((x, x_AOI)).max() - y_max = np.concatenate((y, y_AOI)).max() - exclude_low_fxy = False # don't exclude any points because we need to analyze all grid points in the AOI + # Create an AOI GeoDataFrame and project it to the same UTM CRS + aoi_gdf_utm = gpd.GeoDataFrame( + geometry=gpd.points_from_xy(AOI_lon, AOI_lat), + crs="EPSG:4326" + ).to_crs(utm_crs) + + # Combine dataframes to extract the collective bounding box limits + combined_gdf = gpd.GeoDataFrame(geometry=gpd.pd.concat([catalog_gdf_utm.geometry, aoi_gdf_utm.geometry])) + x_min, y_min, x_max, y_max = combined_gdf.total_bounds else: use_AOI = False + x_min, y_min, x_max, y_max = catalog_gdf_utm.total_bounds + + # round up grid dimensions + grid_x_min = (x_min // grid_dim) * grid_dim + grid_x_max = ((x_max + grid_dim - 1) // grid_dim) * grid_dim + grid_y_min = (y_min // grid_dim) * grid_dim + grid_y_max = ((y_max + grid_dim - 1) // grid_dim) * grid_dim - # define corners of grid based on global dataset - x_min = x.min() - y_min = y.min() - x_max = x.max() - y_max = y.max() - - grid_x_max = int(ceil(x_max / grid_dim) * grid_dim) - grid_x_min = int(floor(x_min / grid_dim) * grid_dim) - grid_y_max = int(ceil(y_max / grid_dim) * grid_dim) - grid_y_min = int(floor(y_min / grid_dim) * grid_dim) - - - - # rectangular grid - nx = int((grid_x_max - grid_x_min) / grid_dim) + 1 - ny = int((grid_y_max - grid_y_min) / grid_dim) + 1 - - # ensure a square grid is used - if nx > ny: # enlarge y dimension to match x - ny = nx - grid_y_max = int(grid_y_min + (ny - 1) * grid_dim) - - else: # enlarge x dimension to match y - nx = ny - grid_x_max = int(grid_x_min + (nx - 1) * grid_dim) - - # update grid extent in lat/lon - grid_lat_max, grid_lon_max = utm.to_latlon(grid_x_max, grid_y_max, utm_zone_number, utm_zone_letter) - grid_lat_min, grid_lon_min = utm.to_latlon(grid_x_min, grid_y_min, utm_zone_number, utm_zone_letter) + # expand extent until it is square + ext_w, ext_h = grid_x_max - grid_x_min, grid_y_max - grid_y_min + delta = abs(ext_w - ext_h) / 2 - # new x and y range - x_range = np.linspace(grid_x_min, grid_x_max, nx) - y_range = np.linspace(grid_y_min, grid_y_max, ny) - - logger.debug(f"Grid X range: {x_range}, Y range: {y_range}") + # Shift the shorter axis outward symmetrically + grid_x_min, grid_x_max = (grid_x_min - delta, grid_x_max + delta) if ext_w < ext_h else (grid_x_min, grid_x_max) + grid_y_min, grid_y_max = (grid_y_min - delta, grid_y_max + delta) if ext_h < ext_w else (grid_y_min, grid_y_max) + # make grid points + x_range = np.arange(grid_x_min, grid_x_max + grid_dim, grid_dim) + nx = len(x_range) + y_range = np.arange(grid_y_min, grid_y_max + grid_dim, grid_dim) + ny = len(y_range) + + X, Y = np.meshgrid(x_range, y_range) + cells = shapely.box(X, Y, X + grid_dim, Y + grid_dim) + grid_gdf_utm = gpd.GeoDataFrame(geometry=cells.flatten(), crs=utm_crs) + grid_gdf_latlon = grid_gdf_utm.to_crs("EPSG:4326") + + logger.debug(f"Grid extent in UTM XY {grid_gdf_utm.total_bounds}") + logger.debug(f"Grid extent in lat lon {grid_gdf_latlon.total_bounds}") + t_windowed = time r_windowed = [[x, y]] @@ -556,7 +552,8 @@ verbose: {verbose}") xy_kale = output_kale[0] xy_kde = output_kde[0] - + grid_gdf_latlon['location_PDF'] = xy_kde[0].flatten() # insert location PDF as a column of the latlon GDF + # plot location PDF xy_kale_km = type(xy_kale)(xy_kale.dataset / 1000) corner = kale.corner(xy_kale_km, quantiles=[0.025, 0.16, 0.50, 0.84, 0.975], cmap='hot')