Update src/seismic_hazard_forecasting.py

update to activity rate calculation based on June 25 email

src/seismic_hazard_forecasting.py

@@ -7,13 +7,13 @@ from matplotlib.ticker import MultipleLocator
 import matplotlib.pyplot as plt
 import numpy as np
 
-global ncp_choice, tcp_max, torder_min, torder_max, AOI_lat, AOI_lon
+global ncp_choice, tcp_max, torder_min, torder_max
 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([51.48, 51.54])
+#AOI_lon = np.array([16.15, 16.24])
 #AOI_lat = np.array([None, None])
 #AOI_lon = np.array([None, None])
 
@@ -152,19 +152,19 @@ def apply_beast(act_rate):
     prob : A list of probabilities and change points extracted from BEAST results.
     """
 
-    mirror_len = int(np.ceil(0.20 * len(act_rate)))
-    left_mirror = act_rate[:mirror_len][::-1]
-    right_mirror = act_rate[-mirror_len:][::-1]
-    act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror])
+    #mirror_len = int(np.ceil(0.20 * len(act_rate)))
+    #left_mirror = act_rate[:mirror_len][::-1]
+    #right_mirror = act_rate[-mirror_len:][::-1]
+    #act_rate_mirrored = np.concatenate([left_mirror, act_rate, right_mirror])
 
     mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15))
-    beast_result = rb.beast(act_rate_mirrored, period=0,
+    beast_result = rb.beast(act_rate, period=0,
                             tcp_minmax=[0, tcp_max],
                             torder_minmax=[torder_min, torder_max],
                             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:
@@ -187,11 +187,20 @@ def apply_beast(act_rate):
             return beast_result, np.array([])
 
     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)]
     cps = beast_result.trend.cp[:ncp]
 
     # Discard mirrored zone changepoints and correct indices
-    valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate))
-    cps = cps[valid_mask] - mirror_len
+    #valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate))
+    #cps = cps[valid_mask] - mirror_len
+
+    # Discard changepoints too close to the start or end (artifacts of mirroring).
+    # bins_after_cp / bins_before_cp set the minimum buffer bins required at each end.
+    #bins_before_cp = 2
+    #bins_after_cp = 2
+    #if len(cps) > 0:
+    #    cps = cps[(cps >= bins_before_cp) & (cps <= len(act_rate) - bins_after_cp)]
 
     return beast_result, np.sort(cps)
 
@@ -253,8 +262,7 @@ def bins_and_beast(dates, unit, bin_dur, multiplicator):
 
 def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max,
          m_kde_method, xy_select, grid_dim, xy_win_method, rate_select, time_win_duration,
-         forecast_select, custom_rate, forecast_len, time_unit, model, products_string, verbose):
-#         forecast_select, custom_rate, forecast_len, time_unit, AOI_extent, model, products_string, verbose):
+         forecast_select, custom_rate, forecast_len, time_unit, AOI_extent, model, products_string, verbose):
     """
     Python application that reads an earthquake catalog and performs seismic hazard forecasting.
     Arguments:
@@ -344,8 +352,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     exclude_low_fxy = False # skip low probability areas of the map
     thresh_fxy = 1e-3  # minimum fxy value (location PDF) needed to do PGA estimation (to skip low probability areas); also should scale according to number of grid points
 
-#    AOI_lat = np.array(AOI_extent[:2]) 
-#    AOI_lon = np.array(AOI_extent[2:])
+    AOI_lat = np.array(AOI_extent[:2]) 
+    AOI_lon = np.array(AOI_extent[2:])
     
     # log user selections
     logger.debug(f"User input files\n Catalog: {catalog_file}\n Mc: {mc_file}\n Mag_PDF: {pdf_file}\n Mag: {m_file}")
@@ -746,6 +754,17 @@ verbose: {verbose}")
             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)
+    
+            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()}")
         else:
             indices_filtered = indices
         
@@ -842,7 +861,35 @@ verbose: {verbose}")
         #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)):
+
+            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)]
+            
             vmin = np.nanmin(iml_grid_prep[j])
             vmax = np.nanmax(iml_grid_prep[j])
 
@@ -875,11 +922,14 @@ verbose: {verbose}")
             # Create an image from the grid
             cmap_name = 'YlOrRd'
             cmap = plt.get_cmap(cmap_name)
-            fig, ax = plt.subplots(figsize=(6, 6))
-            fig.add_axes([0, 0, 1, 1])
+            #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"

src/shf_wrapper.py

@@ -31,6 +31,9 @@ def main(argv):
             return False
         else:
             raise argparse.ArgumentTypeError("Boolean value expected.")
+        
+    def float_or_none(v):
+        return None if v.lower() == "none" else float(v)
 
     parser = argparse.ArgumentParser()
 
@@ -55,6 +58,7 @@ def main(argv):
     parser.add_argument("--time_unit", type=str)
     parser.add_argument("--model", type=str)
     parser.add_argument("--products_string", type=str)
+    parser.add_argument("--AOI_extent", nargs=4, type=float_or_none, default=[None] * 4, required=False)
     parser.add_argument("--verbose", type=str2bool)
 
     args = parser.parse_args()