Replaced usage of removed variable

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

@@ -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,
@@ -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
@@ -310,26 +306,23 @@ 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
     import xml.etree.ElementTree as ET
     import json
     import multiprocessing as mp
+    import geopandas as gpd
+    import shapely
 
     logger = getDefaultLogger('igfash')
 
@@ -349,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:])
     
@@ -362,7 +352,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
 xy_select: {xy_select}\n grid_dim: {grid_dim}\n xy_win_method: {xy_win_method}\n rate_select: {rate_select}\n time_win_duration: {time_win_duration}\n \
 forecast_select: {forecast_select}\n custom_rate: {custom_rate}\n forecast_len: {forecast_len}\n time_unit: {time_unit}\n model: {model}\n products: {products_string}\n \
 verbose: {verbose}")
-
+    logger.debug(f"Area of interest selected by user - Latitude: {AOI_lat}, Longitude: {AOI_lon}")
     # print key package version numbers
     logger.debug(f"Python version {sys.version}")
     logger.debug(f"Numpy version {version('numpy')}")
@@ -471,68 +461,65 @@ 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
 
-            # define corners of grid based on global dataset
-            x_min = x.min()
-            y_min = y.min()
-            x_max = x.max()
-            y_max = y.max()
+        # 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
         
-        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)
+        # 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
         
+        # Shift the shorter axis outward symmetrically
+        grid_x_min, grid_x_max = (grid_x_min - delta, grid_x_max + delta) if ext_w < ext_h else (grid_x_min, grid_x_max)
+        grid_y_min, grid_y_max = (grid_y_min - delta, grid_y_max + delta) if ext_h < ext_w else (grid_y_min, grid_y_max)
         
+        # make grid points
+        x_range = np.arange(grid_x_min, grid_x_max + grid_dim, grid_dim)
+        nx = len(x_range)
+        y_range = np.arange(grid_y_min, grid_y_max + grid_dim, grid_dim)
+        ny = len(y_range)
 
-        # rectangular grid
-        nx = int((grid_x_max - grid_x_min) / grid_dim) + 1
-        ny = int((grid_y_max - grid_y_min) / grid_dim) + 1
+        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")
 
-        # 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)
-        
-        # 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}")
+        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,6 +543,7 @@ verbose: {verbose}")
 
         xy_kale = output_kale[0]
         xy_kde = output_kde[0]
+        grid_gdf_latlon['location_PDF'] = xy_kde[0].flatten() # insert location PDF as a column of the latlon GDF
         
         # plot location PDF
         xy_kale_km = type(xy_kale)(xy_kale.dataset / 1000)
@@ -712,74 +700,41 @@ 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
+        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
 
-        # set every grid point to be a receiver
-        grid_shape = xx.shape
-        x_rx = xx.flatten()
-        y_rx = yy.flatten()
+        # Select only cells of the grid that are inside the AOI
+        if use_AOI:
+            centroids_latlon = grid_gdf_latlon.geometry.centroid
             
-        num_points = x_rx.size
-        distances = np.zeros(shape=(num_points, grid_shape[0], grid_shape[1]))
+            # 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])
+            )                               
+        else:
+            grid_gdf_latlon['AOI']=True #set entire grid to be the area of interest
 
-        # 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)
+        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
         
-        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
+        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 = 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)
-
-        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]
-            
-            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
-        
-        fr = fxy.flatten()
+        #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}")
@@ -789,7 +744,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]
@@ -798,9 +753,9 @@ verbose: {verbose}")
 
         else:
             iml_grid_raw = []
-            iter = indices_filtered
-            for imt in products:
+            iter = range(0,len(distances_sel))
             
+            for imt in products:
                 if imt == "PGV":
                     IMT_max = 200 # search interval max for velocity (cm/s)
                 else:
@@ -808,7 +763,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)
@@ -818,127 +773,51 @@ 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
+        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 or exclude_low_fxy:
+            grid_gdf_latlon_clean = grid_gdf_latlon.dropna(subset=[imt]) # remove null values from grid
 
-        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)
+            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
 
-            # Assign the computed values to their exact original 1D index positions
-            iml_grid_flat[indices_filtered] = iml_grid_raw[j]
+            vmin = np.nanmin(z_plot)
+            vmax = np.nanmax(z_plot)
             
-            # 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:
-            #    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
+            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"),
-                            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)
@@ -952,17 +831,17 @@ 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
+            ax.imshow(gradient, aspect='auto', cmap="YlOrRd_r",
+                      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)