From 6fac004caccd8ba7d4d9d57317875d648428bc15 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 10 Jun 2026 18:13:04 +0200
Subject: [PATCH 01/50] Update src/seismic_hazard_forecasting.py

when the 4 values specifying the lat and lon range of the area of interest (AOI) are provided, only do forecasting for grid points within the AOI
---
 src/seismic_hazard_forecasting.py | 59 ++++++++++++++++++++++++++++++++-------
 1 file changed, 49 insertions(+), 10 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 9e43458..ba5f1aa 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -88,9 +88,12 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     else:
         logger.setLevel(logging.INFO)
 
-    exclude_low_fxy = True # skip low probability areas of the map
+    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([51.48, 51.54]) # temporary hard-coding to area of Zelazny Most. To be replaced with user-defined lat and lon range 
+    AOI_lon = np.array([16.15, 16.24])
+    
     # 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}")
     logger.debug(
@@ -215,11 +218,30 @@ verbose: {verbose}")
         utm_zone_letter = u[3]
         logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}")
 
-        # define corners of grid based on global dataset
-        x_min = x.min()
-        y_min = y.min()
-        x_max = x.max()
-        y_max = y.max()
+        if (None not in AOI_lat) and (None not in AOI_lon):
+            use_AOI = True
+            
+            #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
+            
+        else:
+            use_AOI = False
+        
+            # 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)
@@ -439,10 +461,19 @@ verbose: {verbose}")
         else:
             indices = range(0, len(distances))
 
+        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)
+        else:
+            indices_filtered = indices
+        
         fr = fxy.flatten()
 
         # For each receiver compute estimated ground motion values
-        logger.info(f"Estimating ground motion intensity at {len(indices)} grid points...")
+        logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...")
 
         PGA = np.zeros(shape=(nx * ny))
 
@@ -452,7 +483,7 @@ verbose: {verbose}")
 
         if use_pp:         # use dask parallel computing
             mp.set_start_method("fork", force=True)
-            iter = indices
+            iter = indices_filtered
             iml_grid_raw = []  # raw ground motion grids
             for imt in products:
                 logger.info(f"Estimating {imt}")
@@ -471,7 +502,7 @@ verbose: {verbose}")
 
         else:
             iml_grid_raw = []
-            iter = indices
+            iter = indices_filtered
             for imt in products:
 
                 if imt == "PGV":
@@ -500,7 +531,15 @@ verbose: {verbose}")
         iml_grid = [[] for _ in range(len(products))]  # final ground motion grids
         iml_grid_prep = iml_grid.copy()  # temp ground motion grids
 
-        if exclude_low_fxy:
+        if use_AOI:
+            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)):

From 7e89f75c844e1f7994c172668d52b995a0ccde08 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 15 Jun 2026 11:05:04 +0200
Subject: [PATCH 02/50] Update src/seismic_hazard_forecasting.py

put AOI_extent as a parameter of the main function
---
 src/seismic_hazard_forecasting.py | 8 +++++---
 1 file changed, 5 insertions(+), 3 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index ba5f1aa..3602e52 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -3,7 +3,7 @@
 
 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.
     Arguments:
@@ -33,6 +33,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
                     forecasting.
         forecast_len: Length of the forecast for seismic hazard assessment.
         time_unit: Times units for the inputs Time Window Duration, Custom Activity Rate, and Forecast Length.
+        AOI_extent: The forecast geographical area of interest specified as a latitude and longitude range in decimal degrees
+                    in the form [lat_min, lat_max, lon_min, lon_max].
         model: Select from the following ground motion models available. Other models in the Openquake library are
                     available but have not yet been tested.
         products_string: The ground motion intensity types to output. Use a space between names to select more than
@@ -91,8 +93,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([51.48, 51.54]) # temporary hard-coding to area of Zelazny Most. To be replaced with user-defined lat and lon range 
-    AOI_lon = np.array([16.15, 16.24])
+    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}")

From 46ff6b8a6c805a950eb8ff8e6254fe1d650bcd04 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 14:29:44 +0200
Subject: [PATCH 03/50] Update src/seismic_hazard_forecasting.py

use new activity rate forecast method
---
 src/seismic_hazard_forecasting.py | 314 +++++++++++++++++++++++++++++++++++---
 1 file changed, 289 insertions(+), 25 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3602e52..369ed95 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -1,9 +1,246 @@
 # -*- coding: utf-8 -*-
+from eqdist.rate import datenum_to_datetime
+import Rbeast as rb;
+from scipy.stats import bootstrap
+from matplotlib.dates import DateFormatter, AutoDateLocator
+
+def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
+                  bin_dur, unit, multiplicator,
+                  rate_forecast, rate_unc_high, rate_unc_low,
+                  datenum_data, mag_data):
+
+    end_date = bin_edges[-1]
+
+    fig, ax = plt.subplots(figsize=(14, 6))
+    ax.plot(bin_edges_dt[1:], act_rate, '-o', linewidth=2.5, markersize=6.5, label='Activity rate')
+
+    if rate_forecast is not None:
+        next_date = end_date + (bin_dur / multiplicator)
+        ax.plot(datenum_to_datetime(next_date), rate_forecast, 
+                'ro', label='Forecasted Rate', markersize=6.5)
+        ax.plot([bin_edges_dt[-1], datenum_to_datetime(next_date)], [act_rate[-1], rate_forecast], 'r-', linewidth=2.5)
+        ax.vlines(datenum_to_datetime(next_date), rate_unc_low, rate_unc_high, colors='r',
+                  linewidth=2, label='Bootstrap uncertainty')
+
+    ax.xaxis.set_major_locator(AutoDateLocator())
+    ax.xaxis.set_major_formatter(DateFormatter('%d-%b-%Y'))
+    plt.xticks(rotation=45)
+    plt.title(f'Activity rate (Time Unit: {unit}, Bin Duration: {bin_dur} {unit})',fontsize=18)
+    # plt.title(f'Activity rate (Bin Duration: {bin_dur} {unit})',fontsize=18)
+    plt.xlabel('Time (Bin Center Date)', fontsize=16)
+    ax.set_ylabel('Activity rate per selected time period',fontsize=16)
+    plt.grid(True)
+
+    if len(rt) > 0:
+        for i in range(len(rt)):
+            ax.plot(bin_edges_dt[1:][boundaries[i]:boundaries[i+1]],
+                    [rt[i]] * (boundaries[i+1] - boundaries[i]),
+                    linewidth=2, label=f'Rate period {i+1}')
+
+    # ---- Magnitude scatter on right y-axis ----
+    ax2 = ax.twinx()
+    event_dates = [datenum_to_datetime(d) for d in datenum_data]
+
+    #-------------extract magnitude bins from data---------------------
+    mags = np.array(mag_data)
+    min_mag = mags.min()
+    max_mag = mags.max()
+
+    low_thresh = int(np.floor(min_mag))
+    high_thresh = int(np.floor(max_mag))
+
+    thresholds = list(range(low_thresh, high_thresh + 1))
+
+    base_size = 15
+    size_step = 35
+
+    bins_def = []
+    for idx, t in enumerate(thresholds):
+        low = t
+        if idx < len(thresholds) - 1:
+            high = thresholds[idx + 1]
+            label = f'{low:.1f} \u2264 M < {high:.1f}'
+        else:
+            high = np.inf
+            label = f'M \u2265 {low:.1f}'
+        size = base_size + idx * size_step
+        bins_def.append((low, high, size, label))
+
+    for low, high, size, label in bins_def:
+        mask = (mags >= low) & (mags < high)
+        if np.any(mask):
+            sel_dates = [d for d, m in zip(event_dates, mask) if m]
+            sel_mags = mags[mask]
+            ax2.scatter(sel_dates, sel_mags, s=size,
+                        facecolor='purple', edgecolor='black',
+                        alpha=0.15, linewidth=1, label=label)
+
+    ax2.set_ylabel('Magnitude', color='purple',fontsize=16)
+    ax2.yaxis.set_major_locator(MultipleLocator(0.5))
+    ax2.yaxis.set_minor_locator(MultipleLocator(0.1))
+    ax2.spines['right'].set_color('purple')
+    ax2.tick_params(axis='y', colors='purple')
+
+    h1, l1 = ax.get_legend_handles_labels()
+    h2, l2 = ax2.get_legend_handles_labels()
+    handles = h1 + h2
+    labels = l1 + l2
+    n_legend = len(handles)
+
+    ncols = max(1, int(np.ceil(n_legend / 5)))  # ~5 entries per column
+    
+    #-------add 20% headroom above the data to make space for legend------
+    ymin, ymax = ax.get_ylim()
+    ax.set_ylim(ymin, ymax * 1.20) 
+    
+    ax.legend(handles, labels, loc='best',
+          ncol=ncols, borderaxespad=0,framealpha=0.7)
+
+    ax.set_zorder(ax2.get_zorder() + 1) # put scatter plot behind the line plot
+    ax.patch.set_visible(False)
+    
+    fig.tight_layout()
+    plt.show()
+
+def bootstrap_forecast(data):
+
+    window_data=data
+
+    if len(window_data) >= 5:
+        res95 = bootstrap((window_data,), np.mean, confidence_level=0.95,
+                          method='BCa', n_resamples=1000)
+    else:
+        res95 = bootstrap((window_data,), np.mean, confidence_level=0.95,
+                          method='BCa', n_resamples=int(len(window_data) ** len(window_data)))
+
+    forecast = np.mean(res95.bootstrap_distribution)
+    bca_conf95 = res95.confidence_interval
+    return forecast, bca_conf95
+
+def calc_rates(act_rate, cps):
+    """
+    Calculates mean activity rates between changepoints.
+    cps : sorted array of changepoint indices into act_rate
+    Returns rt (list of rates) and segment boundaries
+    """
+    boundaries = [0] + list(cps.astype(int)) + [len(act_rate)]
+    rt = [np.mean(act_rate[boundaries[i]:boundaries[i+1]]) 
+          for i in range(len(boundaries)-1)]
+    return rt, boundaries
+
+def apply_beast(act_rate, **kwargs):
+    """
+    Applies BEAST to the smmothed rate data using different smoothing windows.
+    Input
+    act_rate : The activity rate data array to smooth and apply BEAST.
+    Output
+    out : A list of BEAST results for each smoothed rate array.
+    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])
+
+    mcmc_th = int(np.clip(np.ceil(len(act_rate) / 100), 2, 15))
+    beast_result = rb.beast(act_rate_mirrored, 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
+    if ncp_choice == 'median':
+        ncp = beast_result.trend.ncp_median
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])
+    elif ncp_choice == 'mode':
+        ncp = beast_result.trend.ncp_mode
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])   
+    elif ncp_choice == 'pct90':
+        ncp = beast_result.trend.ncp_pct90
+        if np.isnan(ncp) or ncp == 0:
+            return beast_result, np.array([])
+    else:  # default: median with mode and pct90 fallback
+        ncp = beast_result.trend.ncp_median
+        if np.isnan(ncp) or ncp == 0:
+            ncp = beast_result.trend.mode
+        if np.isnan(ncp) or ncp == 0:
+            ncp = beast_result.trend.ncp_pct90
+        if np.isnan(ncp) or ncp == 0:    
+            return beast_result, np.array([])
+
+    ncp = int(ncp)
+    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
+
+    return beast_result, np.sort(cps)
+
+def bins_and_beast(dates, unit, bin_dur, multiplicator, **kwargs):
+    start_date = dates.min()
+    end_date = dates.max()
+
+    valid_units = ['hours', 'days']
+    if unit not in valid_units:
+        unit = 'days'
+        bin_dur = 15
+    if (end_date - start_date) < 15 and unit == 'days':
+        unit = 'hours'
+        bin_dur = 12
+
+    bin_edges = [end_date]
+    while bin_edges[-1] > start_date:
+        bin_edges.append(bin_edges[-1] - (bin_dur / multiplicator))
+    bin_edges = bin_edges[::-1]
+
+    #-------Drop first bin or keep it if >80% of set duration------
+    first_width_days = bin_edges[1] - start_date
+    first_width_units = first_width_days * multiplicator
+
+    if first_width_units >= 0.8 * bin_dur:
+        bin_edges[0] = start_date  # edge of first bin is at data start
+    else:
+        bin_edges = bin_edges[1:]  # drop bin 0 (and its events)
+    
+    #------------Error if remaining bins are fewer than 2------------
+    if len(bin_edges) < 2:
+        raise ValueError(
+            f"Not enough data to form at least one full bin of duration "
+            f"{bin_dur} {unit}(s) after dropping the partial first bin "
+            f"({first_width_units:.2f} {unit}(s), below the 80% threshold). "
+            f"Try a shorter bin_dur or check your input data range."
+        )
+    
+    bin_edges_dt = [datenum_to_datetime(d) for d in bin_edges]
+    bin_counts, _ = np.histogram(dates, bins=bin_edges)
+
+    act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator / bin_dur)
+                for i, count in enumerate(bin_counts)]
+
+    out, cps = apply_beast(act_rate, **kwargs)
+
+    if len(cps) > 0:
+        rt, boundaries = calc_rates(act_rate, cps)
+        print(f'Changepoints detected at bins: {cps}')
+    else:
+        rt = []
+        boundaries = []
+        print('-----------------------------------------------------')
+        print('No changepoints detected by BEAST (Zhao et al., 2019)')
+        print('-----------------------------------------------------')
+
+    return act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, unit
 
 
 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, AOI_extent, model, products_string, verbose):
+         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.
     Arguments:
@@ -93,8 +330,11 @@ 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:])
+    AOI_lat = np.array([51.48, 51.54])
+    AOI_lon = np.array([16.15, 16.24])
+
     
     # 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}")
@@ -344,7 +584,12 @@ verbose: {verbose}")
 
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
-
+        
+        ncp_choice = 'default'
+        tcp_max = 5
+        torder_min = 0
+        torder_max = 1
+        
         time, mag_dummy, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
 
         datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
@@ -364,31 +609,50 @@ verbose: {verbose}")
             logger.error(msg)
             raise Exception(msg)
 
-        # Selects dates in datenum format and procceeds to forecast value
-        start_date = datenum_data[-1] - (2 * time_win_duration / multiplicator)
-        dates_calc = [date for date in datenum_data if start_date <= date <= datenum_data[-1]]
-        forecasts, bca_conf95, rate_mean = bootstrap_forecast_rolling(dates_calc, multiplicator)
+        #-----------data are sorted in case they were not-----------------
+        sorted_pairs = sorted(zip(datenum_data, mag_data), key=lambda x: x[0])
+        datenum_data, mag_data = map(list, zip(*sorted_pairs))
 
-        # FINAL VALUES OF RATE AND ITS UNCERTAINTY IN THE 5-95 PERCENTILE
-        unc_bca05 = [ci.low for ci in bca_conf95];
-        unc_bca95 = [ci.high for ci in bca_conf95]
-        rate_unc_high = multiplicator / np.array(unc_bca05);
-        rate_unc_low = multiplicator / np.array(unc_bca95);
-        rate_forecast = multiplicator / np.median(forecasts)  # [per time unit]
+        #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
+        if len(cps) > 0:
+            rate_forecast = rt[-1]
+            last_cp_bin = int(cps[-1])
+        else:
+            rate_forecast = act_rate[-1]
+            last_cp_bin = len(act_rate) - 1
 
-        # Plot of forecasted activity rate with previous binned activity rate
-        act_rate, bin_counts, bin_edges, out, pprs, rt, idx, u_e = calc_bins(np.array(datenum_data), time_unit,
-                                                                             time_win_duration, dates_calc,
-                                                                             rate_forecast, rate_unc_high, rate_unc_low,
-                                                                             multiplicator, quiet=True, figsize=(14,9))
+        last_cp_datenum = bin_edges[last_cp_bin]
+        dates_calc = [date for date in datenum_data if last_cp_datenum <= date <= datenum_data[-1]]
+        interevent_times = np.diff(dates_calc)
 
-        # Assign probabilities 
-        lambdas, lambdas_perc = lambda_probs(act_rate, dates_calc, bin_edges)
-        lambdas = np.array(lambdas, dtype='d')
-        lambdas_perc = np.array(lambdas_perc, dtype='d')
+        #------------Use BCa for uncertainty intervals-----------------
+        forecast, bca_conf95 = bootstrap_forecast(interevent_times)
+        rate_unc_high = bin_dur / (bca_conf95.low * multiplicator)
+        rate_unc_low = bin_dur / (bca_conf95.high * multiplicator)
+
+        #----------------------Plot------------------------------------
+        plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
+                      bin_dur, time_unit, multiplicator,
+                      rate_forecast, rate_unc_high, rate_unc_low,
+                      datenum_data, mag_data)
+
+        print("\n----------------- Forecast Summary -----------------")
+        print(f"Forecasted activity rate (next {bin_dur} {time_unit}(s)): {rate_forecast:.4f}")
+        print(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
+        print("------------------------------------------------------")
         
-        # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1])
-        logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}")
+        lambdas = np.array(rate_forecast/bin_dur, dtype='d')
+        lambdas_perc = np.array(1.0, dtype='d')
+
+
+
+
+
+
+
+        
+        
+        logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit}")
         np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage",
                    delimiter=',', fmt='%1.4f')
 

From bb3184a126cb2c64243e4b3ad4e51b03dcb101bd Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 14:36:49 +0200
Subject: [PATCH 04/50] Update src/seismic_hazard_forecasting.py

correct activity rate to be per time unit before fed to forecasting
---
 src/seismic_hazard_forecasting.py | 19 +++++--------------
 1 file changed, 5 insertions(+), 14 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 369ed95..e1fa212 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -636,24 +636,15 @@ verbose: {verbose}")
                       rate_forecast, rate_unc_high, rate_unc_low,
                       datenum_data, mag_data)
 
-        print("\n----------------- Forecast Summary -----------------")
-        print(f"Forecasted activity rate (next {bin_dur} {time_unit}(s)): {rate_forecast:.4f}")
-        print(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
-        print("------------------------------------------------------")
+        logger.info("\n----------------- Forecast Summary -----------------")
+        logger.info(f"Forecasted activity rate (next {bin_dur} {time_unit}(s)): {rate_forecast:.4f}")
+        logger.info(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
+        logger.info("------------------------------------------------------")
         
         lambdas = np.array(rate_forecast/bin_dur, dtype='d')
         lambdas_perc = np.array(1.0, dtype='d')
 
-
-
-
-
-
-
-        
-        
-        logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit}")
-        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage",
+        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header=f"Activity Rate (Events per {time_unit})", "percentage",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From b51e5b9f434806b360f1eb6f72f38bbd0e5d9df2 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 14:38:00 +0200
Subject: [PATCH 05/50] Update src/seismic_hazard_forecasting.py

correct mag_data variable name
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index e1fa212..3651871 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -590,7 +590,7 @@ verbose: {verbose}")
         torder_min = 0
         torder_max = 1
         
-        time, mag_dummy, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
+        time, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
 
         datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
 

From 3b797e41cbe3c0ad7ab376138c5fb09f59e205b7 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 14:40:44 +0200
Subject: [PATCH 06/50] Update src/seismic_hazard_forecasting.py

fix syntax error
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3651871..a5e560b 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -644,7 +644,7 @@ verbose: {verbose}")
         lambdas = np.array(rate_forecast/bin_dur, dtype='d')
         lambdas_perc = np.array(1.0, dtype='d')
 
-        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header=f"Activity Rate (Events per {time_unit})", "percentage",
+        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header=f"Activity Rate (Events per {time_unit}), percentage",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From 09cf4b0e6abf1f5b2b1f5dc2bfff86b7e181a2e3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 14:54:20 +0200
Subject: [PATCH 07/50] Update src/seismic_hazard_forecasting.py

fix mag_label
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index a5e560b..429a7c2 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -590,7 +590,7 @@ verbose: {verbose}")
         torder_min = 0
         torder_max = 1
         
-        time, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
+        time, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True)
 
         datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
 

From a00d4d6f521d68240507f826a05f1254cb7df9f3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:02:14 +0200
Subject: [PATCH 08/50] Update src/seismic_hazard_forecasting.py

fix missing call to bin_and_beast
---
 src/seismic_hazard_forecasting.py | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 429a7c2..e770159 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -100,6 +100,7 @@ def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
     ax.patch.set_visible(False)
     
     fig.tight_layout()
+    plt.savefig("activity_rate.png", dpi=600)
     plt.show()
 
 def bootstrap_forecast(data):
@@ -613,6 +614,10 @@ verbose: {verbose}")
         sorted_pairs = sorted(zip(datenum_data, mag_data), key=lambda x: x[0])
         datenum_data, mag_data = map(list, zip(*sorted_pairs))
 
+        #-------split the data into bins and apply BEAST for changepoint detection--------------------
+        act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
+            np.array(datenum_data), time_unit, time_range, multiplicator)
+        
         #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
         if len(cps) > 0:
             rate_forecast = rt[-1]

From 4212037a785474c72ae316f80fad85f7f7f08be5 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:08:17 +0200
Subject: [PATCH 09/50] Update src/seismic_hazard_forecasting.py

fix time_win_duration name
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index e770159..a578a0e 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -616,7 +616,7 @@ verbose: {verbose}")
 
         #-------split the data into bins and apply BEAST for changepoint detection--------------------
         act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
-            np.array(datenum_data), time_unit, time_range, multiplicator)
+            np.array(datenum_data), time_unit, time_win_duration, multiplicator)
         
         #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
         if len(cps) > 0:

From 41a900c3660c19c09d1cc675c4bb4405f7bd98ae Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:09:47 +0200
Subject: [PATCH 10/50] Update src/seismic_hazard_forecasting.py

import numpy
---
 src/seismic_hazard_forecasting.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index a578a0e..75b1b45 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -3,6 +3,7 @@ from eqdist.rate import datenum_to_datetime
 import Rbeast as rb;
 from scipy.stats import bootstrap
 from matplotlib.dates import DateFormatter, AutoDateLocator
+import numpy as np
 
 def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
                   bin_dur, unit, multiplicator,

From fb0013199752679ca4892da86bc694b14cd02c79 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:33:43 +0200
Subject: [PATCH 11/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 75b1b45..9f2b064 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -130,7 +130,7 @@ def calc_rates(act_rate, cps):
           for i in range(len(boundaries)-1)]
     return rt, boundaries
 
-def apply_beast(act_rate, **kwargs):
+def apply_beast(act_rate, *args):
     """
     Applies BEAST to the smmothed rate data using different smoothing windows.
     Input
@@ -183,7 +183,7 @@ def apply_beast(act_rate, **kwargs):
 
     return beast_result, np.sort(cps)
 
-def bins_and_beast(dates, unit, bin_dur, multiplicator, **kwargs):
+def bins_and_beast(dates, unit, bin_dur, multiplicator, *args):
     start_date = dates.min()
     end_date = dates.max()
 
@@ -224,7 +224,7 @@ def bins_and_beast(dates, unit, bin_dur, multiplicator, **kwargs):
     act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator / bin_dur)
                 for i, count in enumerate(bin_counts)]
 
-    out, cps = apply_beast(act_rate, **kwargs)
+    out, cps = apply_beast(act_rate, *args)
 
     if len(cps) > 0:
         rt, boundaries = calc_rates(act_rate, cps)
@@ -617,7 +617,7 @@ verbose: {verbose}")
 
         #-------split the data into bins and apply BEAST for changepoint detection--------------------
         act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
-            np.array(datenum_data), time_unit, time_win_duration, multiplicator)
+            np.array(datenum_data), time_unit, time_win_duration, multiplicator, tcp_max, torder_min, torder_max)
         
         #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
         if len(cps) > 0:

From 037863e975c183536c791e76bdcb49aeaae50b4c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:38:35 +0200
Subject: [PATCH 12/50] Update src/seismic_hazard_forecasting.py

temporary global variables since GUI not ready
---
 src/seismic_hazard_forecasting.py | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 9f2b064..67db6da 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -5,6 +5,14 @@ from scipy.stats import bootstrap
 from matplotlib.dates import DateFormatter, AutoDateLocator
 import numpy as np
 
+global ncp_choice, tcp_max, torder_min, torder_max, AOI_lat, AOI_lon
+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])
+
 def plot_results(act_rate, bin_edges, bin_edges_dt, rt, boundaries,
                   bin_dur, unit, multiplicator,
                   rate_forecast, rate_unc_high, rate_unc_low,
@@ -130,7 +138,7 @@ def calc_rates(act_rate, cps):
           for i in range(len(boundaries)-1)]
     return rt, boundaries
 
-def apply_beast(act_rate, *args):
+def apply_beast(act_rate):
     """
     Applies BEAST to the smmothed rate data using different smoothing windows.
     Input
@@ -183,7 +191,7 @@ def apply_beast(act_rate, *args):
 
     return beast_result, np.sort(cps)
 
-def bins_and_beast(dates, unit, bin_dur, multiplicator, *args):
+def bins_and_beast(dates, unit, bin_dur, multiplicator):
     start_date = dates.min()
     end_date = dates.max()
 
@@ -224,7 +232,7 @@ def bins_and_beast(dates, unit, bin_dur, multiplicator, *args):
     act_rate = [count / ((bin_edges[i + 1] - bin_edges[i]) * multiplicator / bin_dur)
                 for i, count in enumerate(bin_counts)]
 
-    out, cps = apply_beast(act_rate, *args)
+    out, cps = apply_beast(act_rate)
 
     if len(cps) > 0:
         rt, boundaries = calc_rates(act_rate, cps)
@@ -334,9 +342,6 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
 
 #    AOI_lat = np.array(AOI_extent[:2]) 
 #    AOI_lon = np.array(AOI_extent[2:])
-    AOI_lat = np.array([51.48, 51.54])
-    AOI_lon = np.array([16.15, 16.24])
-
     
     # 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}")
@@ -586,12 +591,7 @@ verbose: {verbose}")
 
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
-        
-        ncp_choice = 'default'
-        tcp_max = 5
-        torder_min = 0
-        torder_max = 1
-        
+                
         time, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True)
 
         datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
@@ -617,7 +617,7 @@ verbose: {verbose}")
 
         #-------split the data into bins and apply BEAST for changepoint detection--------------------
         act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
-            np.array(datenum_data), time_unit, time_win_duration, multiplicator, tcp_max, torder_min, torder_max)
+            np.array(datenum_data), time_unit, time_win_duration, multiplicator)
         
         #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
         if len(cps) > 0:

From f50a315ed7dda4b45229182a5a972393a326e73c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:39:41 +0200
Subject: [PATCH 13/50] Update src/seismic_hazard_forecasting.py

import matplotlib
---
 src/seismic_hazard_forecasting.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 67db6da..3b9fa7a 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -3,6 +3,7 @@ from eqdist.rate import datenum_to_datetime
 import Rbeast as rb;
 from scipy.stats import bootstrap
 from matplotlib.dates import DateFormatter, AutoDateLocator
+import matplotlib.pyplot as plt
 import numpy as np
 
 global ncp_choice, tcp_max, torder_min, torder_max, AOI_lat, AOI_lon

From 01ec215124c78a0099077a4812467736ea118c1d Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:41:33 +0200
Subject: [PATCH 14/50] Update src/seismic_hazard_forecasting.py

missing from matplotlib.ticker import MultipleLocator
---
 src/seismic_hazard_forecasting.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3b9fa7a..7ba0d77 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -3,6 +3,7 @@ from eqdist.rate import datenum_to_datetime
 import Rbeast as rb;
 from scipy.stats import bootstrap
 from matplotlib.dates import DateFormatter, AutoDateLocator
+from matplotlib.ticker import MultipleLocator
 import matplotlib.pyplot as plt
 import numpy as np
 

From e7a344bca32a0a97439cb8364d1a30e246fbddbf Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:45:23 +0200
Subject: [PATCH 15/50] Update src/seismic_hazard_forecasting.py

remove "percentage" in activity rate csv
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 7ba0d77..5b1bd86 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -652,7 +652,7 @@ verbose: {verbose}")
         lambdas = np.array(rate_forecast/bin_dur, dtype='d')
         lambdas_perc = np.array(1.0, dtype='d')
 
-        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header=f"Activity Rate (Events per {time_unit}), percentage",
+        np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit})",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From 0a09a2dc00d1082f51980921b0358066853946f3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:50:50 +0200
Subject: [PATCH 16/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 5b1bd86..3de67d9 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -652,7 +652,7 @@ verbose: {verbose}")
         lambdas = np.array(rate_forecast/bin_dur, dtype='d')
         lambdas_perc = np.array(1.0, dtype='d')
 
-        np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit})",
+        np.savetxt('activity_rate.csv', np.atleast_1d(lambdas), header=f"Activity Rate (Events per {time_unit})",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From 59955cf08596a514b0ceb4b6bfecb4fce5f5329c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:52:18 +0200
Subject: [PATCH 17/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 3de67d9..b980171 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -652,7 +652,7 @@ verbose: {verbose}")
         lambdas = np.array(rate_forecast/bin_dur, dtype='d')
         lambdas_perc = np.array(1.0, dtype='d')
 
-        np.savetxt('activity_rate.csv', np.atleast_1d(lambdas), header=f"Activity Rate (Events per {time_unit})",
+        np.savetxt('activity_rate.csv', np.atleast_1d(lambdas), header=f"Activity Rate (Events per {time_unit[:-1]})",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From 18e3fd0ca31a891af61ff4533c7c025eeaae8aa6 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:54:49 +0200
Subject: [PATCH 18/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index b980171..429fb37 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -679,7 +679,7 @@ verbose: {verbose}")
             logger.error(msg)
             raise Exception(msg)
         
-        if lambdas[0] == None:
+        if lambdas == None:
             msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..."
             logger.error(msg)
             raise Exception(msg)

From 504b553b9aa4540ad2625474231301aa0de9f4f2 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 15:59:50 +0200
Subject: [PATCH 19/50] Update src/seismic_hazard_forecasting.py

make lambdas at least 1d
---
 src/seismic_hazard_forecasting.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 429fb37..4c7821d 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -649,10 +649,10 @@ verbose: {verbose}")
         logger.info(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
         logger.info("------------------------------------------------------")
         
-        lambdas = np.array(rate_forecast/bin_dur, dtype='d')
+        lambdas = np.atleast_1d(np.array(rate_forecast/bin_dur, dtype='d'))
         lambdas_perc = np.array(1.0, dtype='d')
 
-        np.savetxt('activity_rate.csv', np.atleast_1d(lambdas), header=f"Activity Rate (Events per {time_unit[:-1]})",
+        np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit[:-1]})",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:

From 3ea88e0eb4e18e01cf8c36192bd0a3f159e64137 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 16:05:37 +0200
Subject: [PATCH 20/50] Update src/seismic_hazard_forecasting.py

fix array size of lambdas and lambdas_perc
---
 src/seismic_hazard_forecasting.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 4c7821d..a144433 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -649,9 +649,9 @@ verbose: {verbose}")
         logger.info(f"95% BCa confidence interval: [{rate_unc_low:.4f}, {rate_unc_high:.4f}]")
         logger.info("------------------------------------------------------")
         
-        lambdas = np.atleast_1d(np.array(rate_forecast/bin_dur, dtype='d'))
-        lambdas_perc = np.array(1.0, dtype='d')
-
+        lambdas = np.array([rate_forecast/bin_dur], dtype='d')
+        lambdas_perc = np.array([1], dtype='d')
+        
         np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit[:-1]})",
                    delimiter=',', fmt='%1.4f')
 

From 9a0b9444a378e22cfc7e458505605a7c7355da56 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 17:01:55 +0200
Subject: [PATCH 21/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index a144433..99b710c 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -14,6 +14,8 @@ 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,

From 9f78e5681c6b79971f0e86b23004cf83b9a89301 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Fri, 19 Jun 2026 18:07:55 +0200
Subject: [PATCH 22/50] Update src/seismic_hazard_forecasting.py

use upscale=1 while testing
---
 src/seismic_hazard_forecasting.py | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 99b710c..0799d15 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -14,8 +14,8 @@ 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])
+#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,
@@ -826,7 +826,7 @@ verbose: {verbose}")
 
             # 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 = 4
+                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

From 89f3f70c62678b9576e1c5840db2eda37437103b Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 12:28:23 +0200
Subject: [PATCH 23/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 0799d15..d8930f3 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -12,10 +12,10 @@ 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])
+#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,

From b7ee8d52ab02d18c6d78efd0334bdcba2a3b88a6 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 14:08:10 +0200
Subject: [PATCH 24/50] Update src/seismic_hazard_forecasting.py

log message about AOI activation and fix orientation if AOI selected
---
 src/seismic_hazard_forecasting.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index d8930f3..7cc3d0a 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -473,7 +473,7 @@ verbose: {verbose}")
 
         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]
@@ -836,6 +836,10 @@ verbose: {verbose}")
 
             iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
 
+            # correct orientation is AOI is used
+            if use_AOI:
+                iml_grid_hd = np.rot90(iml_grid_hd, k=-1)
+                
             #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
             vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x))
 

From 0ef41d72f67c5b8711c0f7e34f5406e23a993b3b Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 19:06:19 +0200
Subject: [PATCH 25/50] Update src/seismic_hazard_forecasting.py

fix orientation issue attempt 1
---
 src/seismic_hazard_forecasting.py | 15 ++++++++-------
 1 file changed, 8 insertions(+), 7 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 7cc3d0a..c3b0426 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -699,18 +699,22 @@ verbose: {verbose}")
         fxy = xy_kde[0]
         logger.debug(f"Normalization check; sum of all f(x,y) values = {np.sum(fxy)}")
 
-        xx, yy = np.meshgrid(x_range, y_range, indexing='ij')  # grid points
+        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))
+        #distances = np.zeros(shape=(nx * ny, nx, ny))
         rx_lat = np.zeros(nx * ny)
         rx_lon = np.zeros(nx * ny)
 
-        for i in range(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)
@@ -726,7 +730,7 @@ verbose: {verbose}")
         if exclude_low_fxy:
             indices = list(np.where(fxy.flatten() > thresh_fxy)[0])
         else:
-            indices = range(0, len(distances))
+            indices = np.arange(num_points)
 
         if use_AOI:
             # Filter out receivers outside the AOI; Find indices where values are OUTSIDE the AOI
@@ -836,9 +840,6 @@ verbose: {verbose}")
 
             iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
 
-            # correct orientation is AOI is used
-            if use_AOI:
-                iml_grid_hd = np.rot90(iml_grid_hd, k=-1)
                 
             #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
             vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x))

From 3d4f9138d73a34054dfc33327a43310c35c4a7b3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 19:23:13 +0200
Subject: [PATCH 26/50] Update src/seismic_hazard_forecasting.py

fix attempt 2
---
 src/seismic_hazard_forecasting.py | 13 +++++++++++++
 1 file changed, 13 insertions(+)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index c3b0426..1b7880e 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -803,6 +803,19 @@ verbose: {verbose}")
         iml_grid_prep = iml_grid.copy()  # temp ground motion grids
 
         if use_AOI:
+
+            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)):

From 337457d49cf0d5ef22981e8eda7ebe46ad0b5f99 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 22:58:10 +0200
Subject: [PATCH 27/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 21 ++++++++++-----------
 1 file changed, 10 insertions(+), 11 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 1b7880e..ee252b0 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -767,7 +767,7 @@ verbose: {verbose}")
                 imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, 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]
+                                                            rtol=0.1, use_cython=False) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
@@ -837,19 +837,18 @@ verbose: {verbose}")
         for j in range(0, len(products)):
             vmin = min(x for x in iml_grid_prep[j] if x is not np.nan)
             vmax = max(x for x in iml_grid_prep[j] if x is not np.nan)
-            iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(
-                dtype=np.float64)  # this reduces values to 8 decimal places
+            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
+            #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
 

From 144f7e1fcdb623c103062dbfcf89369c03933152 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 23:04:07 +0200
Subject: [PATCH 28/50] Update src/seismic_hazard_forecasting.py

remove events below mc for activity rate
---
 src/seismic_hazard_forecasting.py | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index ee252b0..c9aa165 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -596,10 +596,13 @@ verbose: {verbose}")
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
                 
-        time, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True)
-
-        datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
+        datenum_data, mag_data, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, mag_label=mag_label, output_datenum=True)
 
+        if trim_to_mc:
+            indices = np.argwhere(mag_data < mc)
+            mag_data = np.delete(mag_data, indices)
+            datenum_data = np.delete(datenum_data, indices)
+        
         if time_unit == 'hours':
             multiplicator = 24
         elif time_unit == 'days':

From d5e5435d8348b5d9c9db2312b6e799b8df8541e4 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Mon, 22 Jun 2026 23:36:33 +0200
Subject: [PATCH 29/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index c9aa165..cef5655 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -749,8 +749,6 @@ verbose: {verbose}")
         # For each receiver compute estimated ground motion values
         logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...")
 
-        PGA = np.zeros(shape=(nx * ny))
-
         start = timer()
 
         use_pp = True
@@ -796,6 +794,8 @@ verbose: {verbose}")
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
+        logger.info(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)

From 105a36893aa7d296f0081e92360b71c2f673e7a3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 00:04:34 +0200
Subject: [PATCH 30/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 34 +++++++++++++++++-----------------
 1 file changed, 17 insertions(+), 17 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index cef5655..12ee04d 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -341,7 +341,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     else:
         logger.setLevel(logging.INFO)
 
-    exclude_low_fxy = False # skip low probability areas of the map
+    exclude_low_fxy = True # 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]) 
@@ -805,7 +805,7 @@ verbose: {verbose}")
         iml_grid = [[] for _ in range(len(products))]  # final ground motion grids
         iml_grid_prep = iml_grid.copy()  # temp ground motion grids
 
-        if use_AOI:
+        if use_AOI or exclude_low_fxy:
 
             for j in range(0, len(products)):
                 # Reassemble the grid cleanly using the original shape
@@ -819,21 +819,21 @@ verbose: {verbose}")
                 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
+            #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
 

From c7f1dfa54848573e55b658ece08e27cc4e0cb5e3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 00:40:16 +0200
Subject: [PATCH 31/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 19 ++++++++++---------
 1 file changed, 10 insertions(+), 9 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 12ee04d..a544319 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -768,7 +768,7 @@ verbose: {verbose}")
                 imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, 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=False) for i in iter]
+                                                            rtol=0.1, use_cython=True) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
@@ -838,10 +838,11 @@ verbose: {verbose}")
             iml_grid_prep = iml_grid_raw
 
         for j in range(0, len(products)):
-            vmin = min(x for x in iml_grid_prep[j] if x is not np.nan)
-            vmax = max(x for x in iml_grid_prep[j] if x is not np.nan)
-            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
+            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:
@@ -851,13 +852,13 @@ verbose: {verbose}")
             #    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_tmp
 
-            iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
+            #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 = max(x for x in iml_grid_hd.flatten() if not isnan(x))
+            vmax_hd = np.nanmax(iml_grid_hd.flatten())
 
             # generate image overlay
             north, south = lat.max(), lat.min()  # Latitude range

From e3d7fca55c7a9e16e83b837bedd6145ae6a3470c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 02:05:27 +0200
Subject: [PATCH 32/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index a544319..8a83a18 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -341,7 +341,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     else:
         logger.setLevel(logging.INFO)
 
-    exclude_low_fxy = True # skip low probability areas of the map
+    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]) 

From fa192fd7a0e233e3f369003f778071e0fbb0d37b Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 02:20:53 +0200
Subject: [PATCH 33/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 8a83a18..4606470 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -768,7 +768,7 @@ verbose: {verbose}")
                 imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, 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]
+                                                            rtol=0.01, use_cython=True) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 

From 62e551ff6a7f7f64103ea77f3a3cf9d99385c64f Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 02:52:42 +0200
Subject: [PATCH 34/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 4606470..fda2af5 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -858,7 +858,7 @@ verbose: {verbose}")
 
             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.flatten())
+            vmax_hd = np.nanmax(iml_grid_hd)
 
             # generate image overlay
             north, south = lat.max(), lat.min()  # Latitude range

From 13b463a15a8ef6b1f18ba737abe60ab7f9cec72c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 02:53:44 +0200
Subject: [PATCH 35/50] Update src/seismic_hazard_forecasting.py

try distance in metres
---
 src/seismic_hazard_forecasting.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index fda2af5..b9e3ed7 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -727,7 +727,7 @@ verbose: {verbose}")
                                                  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
+        #distances = distances/1000.0
 
         # compute ground motion only at grid points that have minimum probability density of thresh_fxy
         if exclude_low_fxy:

From 92b6eb4880244272f9f4b0639a7e6c62e4d1f491 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 03:21:53 +0200
Subject: [PATCH 36/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 28 ++++++++++++++--------------
 1 file changed, 14 insertions(+), 14 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index b9e3ed7..8612321 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -768,7 +768,7 @@ verbose: {verbose}")
                 imls = [dask.delayed(compute_IMT_exceedance)(rx_lat[i], rx_lon[i], distances[i].flatten(), fr, 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.01, use_cython=True) for i in iter]
+                                                            rtol=0.1, use_cython=True) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
@@ -805,18 +805,18 @@ verbose: {verbose}")
         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:
+        #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 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:
@@ -834,8 +834,8 @@ verbose: {verbose}")
         #        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
+        #else:
+        #    iml_grid_prep = iml_grid_raw
 
         for j in range(0, len(products)):
             vmin = np.nanmin(iml_grid_prep[j])

From 3576b9eb963c3949909ab6bbe152c84dd307d5e1 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 03:58:33 +0200
Subject: [PATCH 37/50] Update src/seismic_hazard_forecasting.py

activate test zone AOI
---
 src/seismic_hazard_forecasting.py | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 8612321..ea718d4 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -12,10 +12,10 @@ 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])
+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,

From 72dc427a0ea79aed0b4bd6221f93d5669028e5b3 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 10:56:54 +0200
Subject: [PATCH 38/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index ea718d4..b6d2daf 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -501,8 +501,7 @@ verbose: {verbose}")
         grid_y_max = int(ceil(y_max / grid_dim) * grid_dim)
         grid_y_min = int(floor(y_min / grid_dim) * grid_dim)
 
-        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)
+       
 
         # rectangular grid
         nx = int((grid_x_max - grid_x_min) / grid_dim) + 1
@@ -517,6 +516,10 @@ verbose: {verbose}")
             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)
@@ -861,8 +864,8 @@ verbose: {verbose}")
             vmax_hd = np.nanmax(iml_grid_hd)
 
             # generate image overlay
-            north, south = lat.max(), lat.min()  # Latitude range
-            east, west = lon.max(), lon.min()  # Longitude range
+            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])]

From 73fd54cc769c147c45ad92ba19f5a92d07c48e1f Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 13:09:39 +0200
Subject: [PATCH 39/50] Update src/seismic_hazard_forecasting.py

add more logging
---
 src/seismic_hazard_forecasting.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index b6d2daf..1385f43 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -524,6 +524,8 @@ verbose: {verbose}")
         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]]
 
@@ -797,7 +799,7 @@ verbose: {verbose}")
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
-        logger.info(f"IMT values: {iml_grid_raw[0]}")
+        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."
@@ -868,6 +870,8 @@ verbose: {verbose}")
             east, west = grid_lon_max, grid_lon_min  # Longitude range
             bounds = [[south, west], [north, east]]
 
+            logger.debug(f"Grid extent North {north}, South {south}, East {east}, West {west}")
+            
             map_center = [np.mean([north, south]), np.mean([east, west])]
 
             # Create an image from the grid

From bb70ebbe24b33e52f1b03eeb320c512231753534 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 14:33:55 +0200
Subject: [PATCH 40/50] Update src/seismic_hazard_forecasting.py

stretch image to exact borders of svg file
---
 src/seismic_hazard_forecasting.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 1385f43..4b768fb 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -869,8 +869,6 @@ verbose: {verbose}")
             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]]
-
-            logger.debug(f"Grid extent North {north}, South {south}, East {east}, West {west}")
             
             map_center = [np.mean([north, south]), np.mean([east, west])]
 
@@ -878,13 +876,14 @@ verbose: {verbose}")
             cmap_name = 'YlOrRd'
             cmap = plt.get_cmap(cmap_name)
             fig, ax = plt.subplots(figsize=(6, 6))
-            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear')
+            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')
 
             # Save the figure
             fig.canvas.draw()
             overlay_filename = f"overlay_{j}.svg"
-            plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True)
+            plt.savefig(overlay_filename, pad_inches=0, transparent=True)
             plt.close(fig)
 
             # Embed geographic bounding box into the SVG

From 13509d37c186ddfc225d91816cc2433c81a27b76 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 17:16:07 +0200
Subject: [PATCH 41/50] Update src/shf_wrapper.py

update wrapper
---
 src/shf_wrapper.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/src/shf_wrapper.py b/src/shf_wrapper.py
index 7675fb3..3a033cf 100644
--- a/src/shf_wrapper.py
+++ b/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()

From e1924587df40780ebb3f55628c2160cc2f901fec Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 23 Jun 2026 17:18:52 +0200
Subject: [PATCH 42/50] Update src/seismic_hazard_forecasting.py

use AOI from user interface
---
 src/seismic_hazard_forecasting.py | 9 ++++-----
 1 file changed, 4 insertions(+), 5 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 4b768fb..15e2f63 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/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])
 
@@ -253,8 +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):
+         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:

From eae7602a6ae67c064647b0f3a549ce674cf5f73a Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 24 Jun 2026 11:02:53 +0200
Subject: [PATCH 43/50] Update src/seismic_hazard_forecasting.py

---
 src/seismic_hazard_forecasting.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 15e2f63..58a8968 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -343,8 +343,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}")

From 27396a7d824b2253c86c50b5cec298358c5f514c Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 24 Jun 2026 16:38:37 +0200
Subject: [PATCH 44/50] Update src/seismic_hazard_forecasting.py

update apply_best and plotting
---
 src/seismic_hazard_forecasting.py | 54 ++++++++++++++++++++++++++++++++++++---
 1 file changed, 51 insertions(+), 3 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 58a8968..b74cf34 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -193,6 +193,12 @@ 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):
@@ -745,6 +751,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
         
@@ -841,7 +858,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])
 
@@ -874,11 +919,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"

From 75a488a54ea061041f4871203ee6cf2c791bd0d6 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Tue, 30 Jun 2026 17:47:42 +0200
Subject: [PATCH 45/50] Update src/seismic_hazard_forecasting.py

update to activity rate calculation based on June 25 email
---
 src/seismic_hazard_forecasting.py | 29 ++++++++++++++++-------------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index b74cf34..5e2b674 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -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,17 +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 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]
+    # 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)
 

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 46/50] 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')

From a56f18dde0aa94fc80909f23f67a8d8021d7ccbe Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 8 Jul 2026 16:12:39 +0200
Subject: [PATCH 47/50] Update src/seismic_hazard_forecasting.py

change IMT plotting method
---
 src/seismic_hazard_forecasting.py | 241 +++++++++++---------------------------
 1 file changed, 70 insertions(+), 171 deletions(-)

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)

From 2ff3f3c8086b061746719a81fb36a44332b63bf0 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 8 Jul 2026 16:42:22 +0200
Subject: [PATCH 48/50] Update src/seismic_hazard_forecasting.py

use dataframe total_bounds as the bounding box of SVG file
---
 src/seismic_hazard_forecasting.py | 26 +++++++++-----------------
 1 file changed, 9 insertions(+), 17 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 95af87a..ef965e9 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -723,15 +723,13 @@ verbose: {verbose}")
             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:
             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
 
+        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
+        
         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
@@ -821,16 +819,10 @@ verbose: {verbose}")
             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]
+            north = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[3]
+            south = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[1]
+            east = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[2]
+            west = grid_gdf_latlon[grid_gdf_latlon['AOI']].total_bounds[0]
                        
             # Embed geographic bounding box into the SVG
             map_bounds = dict(zip(("south", "west", "north", "east"),
@@ -858,7 +850,7 @@ verbose: {verbose}")
                       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)
+            tick_positions = np.linspace(vmin, vmax, num_ticks)
             ax.set_yticks(tick_positions)
             ax.set_yticklabels([f"{tick:.2f}" for tick in tick_positions])  # format tick labels
             ax.set_title(colorbar_title, loc='right', pad=15) 

From b59cf0b5d1ae7cdc553c786194a0f62e2db826c1 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 8 Jul 2026 17:13:30 +0200
Subject: [PATCH 49/50] Update src/seismic_hazard_forecasting.py

cleanup package imports
---
 src/seismic_hazard_forecasting.py | 11 +++--------
 1 file changed, 3 insertions(+), 8 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index ef965e9..572321c 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -306,20 +306,15 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     import logging
     from base_logger import getDefaultLogger
     from timeit import default_timer as timer
-    from math import ceil, floor, isnan
     import numpy as np
     import dask
     import kalepy as kale
-    import utm
-    from skimage.transform import resize
     import igfash
-    from igfash.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf, read_csv
-    from igfash.window import win_CTL, win_CNE
+    from igfash.io import read_mat_cat, read_mat_m, read_mat_mc, read_mat_pdf
+    from igfash.window import win_CNE
     import igfash.kde as kde
     from igfash.gm import compute_IMT_exceedance
-    from igfash.compute import get_cdf, hellinger_dist, cols_to_rows
-    from igfash.rate import lambda_probs, calc_bins, bootstrap_forecast_rolling
-    from igfash.mc import estimate_mc
+    from igfash.compute import get_cdf
     import matplotlib.pyplot as plt
     from matplotlib.ticker import MultipleLocator
     from matplotlib.contour import ContourSet

From 5325c074041a5307ca7e8957d555dc306fdeecd8 Mon Sep 17 00:00:00 2001
From: ftong <95+ftong@noreply.example.org>
Date: Wed, 8 Jul 2026 17:16:21 +0200
Subject: [PATCH 50/50] Update src/seismic_hazard_forecasting.py

AOI feature makes exclude_low_fxy redundant as we now have way to save time by selecting small area to forecast
---
 src/seismic_hazard_forecasting.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 572321c..100dc14 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -342,9 +342,6 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     else:
         logger.setLevel(logging.INFO)
 
-    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:])
     
@@ -482,7 +479,6 @@ verbose: {verbose}")
         # Handle Area of Interest (AOI)
         if (None not in AOI_lat) and (None not in AOI_lon):
             use_AOI = True
-            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(