Update src/seismic_hazard_forecasting.py

try distances/1000

src/seismic_hazard_forecasting.py

@@ -193,12 +193,6 @@ def apply_beast(act_rate):
     valid_mask = (cps > mirror_len) & (cps <= mirror_len + len(act_rate))
     cps = cps[valid_mask] - mirror_len
 
-    # Discard changepoints too close to the end (artifacts of end-mirroring). 
-    # bins_after_cp sets the minimum number of bins required after a changepoint.
-    bins_after_cp = 1
-    if len(cps) > 0:
-        cps = cps[cps <= len(act_rate) - bins_after_cp]
-
     return beast_result, np.sort(cps)
 
 def bins_and_beast(dates, unit, bin_dur, multiplicator):
@@ -259,6 +253,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):
     """
     Python application that reads an earthquake catalog and performs seismic hazard forecasting.
@@ -751,17 +746,6 @@ 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
         
@@ -858,35 +842,7 @@ 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])
 
@@ -919,13 +875,10 @@ 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, ax = plt.subplots(layout=None)
-            ax.margins(0) # clear any data margins
+            fig, ax = plt.subplots(figsize=(6, 6))
+            fig.add_axes([0, 0, 1, 1])
             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()