Update src/seismic_hazard_forecasting.py

use geopandas for coordinate conversion before Location PDF estimation

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,68 +469,66 @@ verbose: {verbose}")
 
         time, mag, lat, lon, depth = read_mat_cat(catalog_file)
 
-        # convert to UTM
+        # Create a GeoDataFrame for the catalog data (initially WGS84 EPSG:4326)
-        u = utm.from_latlon(lat, lon)
+        catalog_gdf = gpd.GeoDataFrame(
-        x = u[0]
+            {'depth': depth, 'time': time},
-        y = u[1]
+            geometry=gpd.points_from_xy(lon, lat),
-        utm_zone_number = u[2]
+            crs="EPSG:4326"
-        utm_zone_letter = u[3]
+        )
-        logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}")
 
+        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
 
-            # define corners of grid based on global dataset
+        # round up grid dimensions 
-            x_min = x.min()
+        grid_x_min = (x_min // grid_dim) * grid_dim
-            y_min = y.min()
+        grid_x_max = ((x_max + grid_dim - 1) // grid_dim) * grid_dim
-            x_max = x.max()
+        grid_y_min = (y_min // grid_dim) * grid_dim
-            y_max = y.max()
+        grid_y_max = ((y_max + grid_dim - 1) // grid_dim) * grid_dim
         
-        grid_x_max = int(ceil(x_max / grid_dim) * grid_dim)
+        # expand extent until it is square
-        grid_x_min = int(floor(x_min / grid_dim) * grid_dim)
+        ext_w, ext_h = grid_x_max - grid_x_min, grid_y_max - grid_y_min
-        grid_y_max = int(ceil(y_max / grid_dim) * grid_dim)
+        delta = abs(ext_w - ext_h) / 2
-        grid_y_min = int(floor(y_min / grid_dim) * grid_dim)
         
+        # 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)
 
-        # rectangular grid
+        X, Y = np.meshgrid(x_range, y_range)
-        nx = int((grid_x_max - grid_x_min) / grid_dim) + 1
+        cells = shapely.box(X, Y, X + grid_dim, Y + grid_dim)
-        ny = int((grid_y_max - grid_y_min) / grid_dim) + 1
+        grid_gdf_utm = gpd.GeoDataFrame(geometry=cells.flatten(), crs=utm_crs)
+        grid_gdf_latlon = grid_gdf_utm.to_crs("EPSG:4326")
 
-        # ensure a square grid is used
+        logger.debug(f"Grid extent in UTM XY {grid_gdf_utm.total_bounds}")
-        if nx > ny:  # enlarge y dimension to match x
+        logger.debug(f"Grid extent in lat lon {grid_gdf_latlon.total_bounds}")
-            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)
-        
-        # 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}")
                 
         t_windowed = time
         r_windowed = [[x, y]]
@@ -556,6 +552,7 @@ 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)