diff --git a/src/seismic_hazard_forecasting.py b/src/seismic_hazard_forecasting.py
index 12e9abe..257c0f2 100644
--- a/src/seismic_hazard_forecasting.py
+++ b/src/seismic_hazard_forecasting.py
@@ -519,13 +519,17 @@ verbose: {verbose}")
                 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 if there are at least 10 grid values with range greater than 0.1
+            # 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
                 iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
                                     mode='reflect', anti_aliasing=False)
+                interp_method = 'bilinear'
+                trim_thresh = vmin
+                iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
             else:
                 iml_grid_hd = iml_grid_tmp
+                interp_method = 'none'
 
             iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
 
@@ -533,8 +537,8 @@ verbose: {verbose}")
             vmax_hd = max(x for x in iml_grid_hd.flatten() if not isnan(x))
 
             # trim edges so the grid is not so blocky
-            trim_thresh = vmin
-            iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
+            #trim_thresh = vmin
+            #iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
 
             # generate image overlay
             north, south = lat.max(), lat.min()  # Latitude range
@@ -547,7 +551,7 @@ 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)
+            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax, interpolation=interp_method)
             ax.axis('off')
 
             # Save the figure