ISEPOS-2450 Added AOI_extent argument to shf_wrapper

src/seismic_hazard_forecasting.py

@@ -1,260 +1,4 @@
 # -*- 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
-from matplotlib.ticker import MultipleLocator
-import matplotlib.pyplot as plt
-import numpy as np
-
-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([None, None])
-#AOI_lon = np.array([None, None])
-
-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.savefig("activity_rate.png", dpi=600)
-    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):
-    """
-    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
-
-    # 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):
-    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)
-
-    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,
@@ -478,7 +222,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]
@@ -506,7 +250,8 @@ 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
@@ -521,16 +266,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)
 
-        logger.debug(f"Grid X range: {x_range}, Y range: {y_range}")
-        
         t_windowed = time
         r_windowed = [[x, y]]
 
@@ -605,14 +344,11 @@ verbose: {verbose}")
 
     elif rate_select:
         logger.info(f"Activity rate modeling selected")
-                
-        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:
+        time, mag_dummy, lat_dummy, lon_dummy, depth_dummy = read_mat_cat(catalog_file, output_datenum=True)
-            indices = np.argwhere(mag_data < mc)
+
-            mag_data = np.delete(mag_data, indices)
+        datenum_data = time  # REMEMBER THE DECIMAL DENOTES DAYS
-            datenum_data = np.delete(datenum_data, indices)
+
-        
         if time_unit == 'hours':
             multiplicator = 24
         elif time_unit == 'days':
@@ -628,46 +364,32 @@ verbose: {verbose}")
             logger.error(msg)
             raise Exception(msg)
 
-        #-----------data are sorted in case they were not-----------------
+        # Selects dates in datenum format and procceeds to forecast value
-        sorted_pairs = sorted(zip(datenum_data, mag_data), key=lambda x: x[0])
+        start_date = datenum_data[-1] - (2 * time_win_duration / multiplicator)
-        datenum_data, mag_data = map(list, zip(*sorted_pairs))
+        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)
 
-        #-------split the data into bins and apply BEAST for changepoint detection--------------------
+        # FINAL VALUES OF RATE AND ITS UNCERTAINTY IN THE 5-95 PERCENTILE
-        act_rate, bin_counts, bin_edges, bin_edges_dt, out, cps, rt, boundaries, bin_dur, time_unit = bins_and_beast(
+        unc_bca05 = [ci.low for ci in bca_conf95];
-            np.array(datenum_data), time_unit, time_win_duration, multiplicator)
+        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]
+
+        # 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))
+
+        # 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')
         
-        #------Forecasted rate is taken from BEAST or is equal to last value if no changepoints detected-----
+        # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1])
-        if len(cps) > 0:
+        logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}")
-            rate_forecast = rt[-1]
+        np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage",
-            last_cp_bin = int(cps[-1])
-        else:
-            rate_forecast = act_rate[-1]
-            last_cp_bin = len(act_rate) - 1
-
-        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)
-
-        #------------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)
-
-        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], dtype='d')
-        
-        np.savetxt('activity_rate.csv', lambdas, header=f"Activity Rate (Events per {time_unit[:-1]})",
                    delimiter=',', fmt='%1.4f')
 
     if forecast_select:
@@ -694,7 +416,7 @@ verbose: {verbose}")
             logger.error(msg)
             raise Exception(msg)
         
-        if lambdas == None:
+        if lambdas[0] == None:
             msg = "Activity rate modeling was not selected and custom activity rate was not provided; cannot continue..."
             logger.error(msg)
             raise Exception(msg)
@@ -712,22 +434,18 @@ 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)  # grid points
+        xx, yy = np.meshgrid(x_range, y_range, indexing='ij')  # 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(num_points):
+        for i in range(nx * ny):
             # 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)
@@ -737,13 +455,13 @@ 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:
             indices = list(np.where(fxy.flatten() > thresh_fxy)[0])
         else:
-            indices = np.arange(num_points)
+            indices = range(0, len(distances))
 
         if use_AOI:
             # Filter out receivers outside the AOI; Find indices where values are OUTSIDE the AOI
@@ -751,17 +469,6 @@ verbose: {verbose}")
             indices_outside_y = np.where((y_rx < y_AOI[0]) | (y_rx > y_AOI[1]))[0]      
             indices_outside_AOI = np.unique(np.concatenate((indices_outside_x, indices_outside_y)))
             indices_filtered = np.setdiff1d(indices, indices_outside_AOI)
-
-            # AOI grid extent
-            AOI_rx_x = x_rx[indices_filtered]
-            AOI_rx_y = y_rx[indices_filtered]
-            
-            AOI_rx_lat, AOI_rx_lon = utm.to_latlon(AOI_rx_x, AOI_rx_y, utm_zone_number, utm_zone_letter)
-    
-            logger.debug(f"Receiver UTM X range: {AOI_rx_x.min()} to {AOI_rx_x.max()}")
-            logger.debug(f"Receiver UTM Y range: {AOI_rx_y.min()} to {AOI_rx_y.max()}")
-            logger.debug(f"Receiver lat range: {AOI_rx_lat.min()} to {AOI_rx_lat.max()}")
-            logger.debug(f"Receiver lon range: {AOI_rx_lon.min()} to {AOI_rx_lon.max()}")
         else:
             indices_filtered = indices
         
@@ -770,6 +477,8 @@ 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
@@ -815,8 +524,6 @@ 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)
@@ -826,111 +533,65 @@ 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:
+            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
 
         for j in range(0, len(products)):
-            # Reassemble the grid cleanly using the original shape
+            vmin = min(x for x in iml_grid_prep[j] if x is not np.nan)
-            # Initialize a flat array filled entirely with NaNs
+            vmax = max(x for x in iml_grid_prep[j] if x is not np.nan)
-            iml_grid_flat = np.full(num_points, np.nan, dtype=np.float64)
+            iml_grid[j] = np.reshape(iml_grid_prep[j], (nx, ny)).astype(
-            
+                dtype=np.float64)  # this reduces values to 8 decimal places
-            # Assign the computed values to their exact original 1D index positions
+            iml_grid_tmp = np.nan_to_num(iml_grid[j])  # change nans to zeroes
-            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)):
-
-            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])
-
-            #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:
+            if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1:
-            #    up_factor = 1
+                up_factor = 4
-            #    iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
+                iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
-            #                        mode='reflect', anti_aliasing=False)
+                                    mode='reflect', anti_aliasing=False)
-            #    trim_thresh = vmin
+                trim_thresh = vmin
-            #    iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
+                iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
-            #else:
+            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 = np.nanmax(iml_grid_hd)
+            vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x))
 
             # generate image overlay
-            north, south = grid_lat_max, grid_lat_min  # Latitude range
+            north, south = lat.max(), lat.min()  # Latitude range
-            east, west = grid_lon_max, grid_lon_min  # Longitude range
+            east, west = lon.max(), 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(figsize=(6, 6))
-            fig, ax = plt.subplots(layout=None)
+            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation='bilinear')
-            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"
-            plt.savefig(overlay_filename, pad_inches=0, transparent=True)
+            plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True)
             plt.close(fig)
 
             # Embed geographic bounding box into the SVG