Merge pull request 'fixes' (#13) from fixes into master

Reviewed-on: #13
Reviewed-by: mieszkomakuch <mieszkomakuch@noreply.example.org>

src/seismic_hazard_forecasting.py

@@ -48,7 +48,7 @@ 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
+    from math import ceil, floor, isnan
     import numpy as np
     import scipy
     import obspy
@@ -70,6 +70,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     import matplotlib.pyplot as plt
     from matplotlib.ticker import MultipleLocator
     from matplotlib.contour import ContourSet
+    import xml.etree.ElementTree as ET
+    import json
 
     logger = getDefaultLogger('igfash')
 
@@ -239,6 +241,9 @@ 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
         ny = int((grid_y_max - grid_y_min) / grid_dim) + 1
@@ -434,22 +439,39 @@ verbose: {verbose}")
 
         PGA = np.zeros(shape=(nx * ny))
 
-        # use dask parallel computing
         start = timer()
-        pbar = ProgressBar()
-        pbar.register()
-        # iter = range(0,len(distances))
-        iter = indices
-        iml_grid_raw = []  # raw ground motion grids
-        for imt in products:
-            logger.info(f"Estimating {imt}")
 
-            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=2.0,
-                                                         rx_label=i) for i in iter]
-            iml = dask.compute(*imls)
-            iml_grid_raw.append(list(iml))
+        use_pp = False
+
+        if use_pp:         # use dask parallel computing
+            pbar = ProgressBar()
+            pbar.register()
+            # iter = range(0,len(distances))
+            iter = indices
+            iml_grid_raw = []  # raw ground motion grids
+            for imt in products:
+                logger.info(f"Estimating {imt}")
+
+                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=2.0, rx_label=i,
+                                                            rtol=0.1, use_cython=False) for i in iter]
+                iml = dask.compute(*imls)
+                iml_grid_raw.append(list(iml))
+
+        else:
+            iml_grid_raw = []
+            iter = indices
+            for imt in products:
+                iml = []
+                for i in iter:
+                    iml_i = 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, imt=imt, IMT_min = 0.0,
+                                                    IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=False)
+                    iml.append(iml_i)
+                    logger.info(f"Estimated {imt} at rx {i} is {iml_i}")
+                iml_grid_raw.append(iml)
+
         end = timer()
         logger.info(f"Ground motion exceedance computation time: {round(end - start, 1)} seconds")
 
@@ -481,6 +503,12 @@ verbose: {verbose}")
                                  mode='reflect', anti_aliasing=False)
             iml_grid_hd[iml_grid_hd == 0.0] = np.nan  # change zeroes back to nan
 
+            # trim edges so the grid is not so blocky
+            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))
+            trim_thresh = vmin
+            iml_grid_hd[iml_grid_hd < trim_thresh] = np.nan
+
             # generate image overlay
             north, south = lat.max(), lat.min()  # Latitude range
             east, west = lon.max(), lon.min()  # Longitude range
@@ -489,17 +517,27 @@ verbose: {verbose}")
             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))
-            ax.imshow(iml_grid_hd, origin='lower', cmap='viridis')
+            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax)
             ax.axis('off')
 
             # Save the figure
             fig.canvas.draw()
-            plt.savefig("overlay_" + str(j) + ".svg", bbox_inches="tight", pad_inches=0, transparent=True)
+            overlay_filename = f"overlay_{j}.svg"
+            plt.savefig(overlay_filename, bbox_inches="tight", pad_inches=0, transparent=True)
             plt.close(fig)
 
+            # 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))))
+            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)
+            logger.info(f"Saved geographic bounds to SVG metadata (data-map-bounds): {overlay_filename} → {map_bounds}")
+
             # Make the color bar
-            cmap_name = 'viridis'
             width = 50
             height = 500
 
@@ -507,16 +545,20 @@ verbose: {verbose}")
             gradient = np.vstack((gradient, gradient)).T
             gradient = np.tile(gradient, (1, width))
 
+            colorbar_title = products[j]
+            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=plt.get_cmap(cmap_name),
-                      extent=[0, 1, vmin, vmax])  # Note: extent order is different for vertical
+            ax.imshow(gradient, aspect='auto', cmap=cmap.reversed(),
+                      extent=[0, 1, vmin, vmax_hd])  # 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, num_ticks)
+            tick_positions = np.linspace(vmin, vmax_hd, num_ticks)
             ax.set_yticks(tick_positions)
             ax.set_yticklabels([f"{tick:.2f}" for tick in tick_positions])  # format tick labels
-            ax.set_title(imt, pad=15)
+            ax.set_title(colorbar_title, loc='right', pad=15) 
             fig.subplots_adjust(left=0.25, right=0.75, bottom=0.05, top=0.95)  # Adjust Layout
             fig.savefig("colorbar_" + str(j) + ".svg", bbox_inches='tight')
             plt.close(fig)

src/shf_wrapper.py

@@ -23,6 +23,15 @@ from seismic_hazard_forecasting import main as shf
 
 
 def main(argv):
+
+    def str2bool(v):
+        if v.lower() in ("True", "TRUE", "yes", "true", "t", "y", "1"):
+            return True
+        elif v.lower() in ("False", "FALSE", "no", "false", "f", "n", "0"):
+            return False
+        else:
+            raise argparse.ArgumentTypeError("Boolean value expected.")
+
     parser = argparse.ArgumentParser()
 
     parser.add_argument("catalog_file", help="Path to input file of type 'catalog'")
@@ -30,23 +39,23 @@ def main(argv):
     parser.add_argument("--pdf_file", help="Path to input file of type 'PDF'", type=str, default=None, required=False)
     parser.add_argument("--m_file", help="Path to input file of type 'm'", type=str, default=None, required=False)
 
-    parser.add_argument("--m_select", type=bool)
+    parser.add_argument("--m_select", type=str2bool)
     parser.add_argument("--mag_label", type=str, default=None, required=False)
     parser.add_argument("--mc", type=float, default=None, required=False)
     parser.add_argument("--m_max", type=float, default=None, required=False)
     parser.add_argument("--m_kde_method", type=str)
-    parser.add_argument("--xy_select", type=bool)
+    parser.add_argument("--xy_select", type=str2bool)
     parser.add_argument("--grid_dim", type=int)
-    parser.add_argument("--xy_win_method", type=bool)
-    parser.add_argument("--rate_select", type=bool)
+    parser.add_argument("--xy_win_method", type=str2bool)
+    parser.add_argument("--rate_select", type=str2bool)
     parser.add_argument("--time_win_duration", type=float)
-    parser.add_argument("--forecast_select", type=bool)
+    parser.add_argument("--forecast_select", type=str2bool)
     parser.add_argument("--custom_rate", type=float, default=None, required=False)
     parser.add_argument("--forecast_len", type=float)
     parser.add_argument("--time_unit", type=str)
     parser.add_argument("--model", type=str)
     parser.add_argument("--products_string", type=str)
-    parser.add_argument("--verbose", type=bool)
+    parser.add_argument("--verbose", type=str2bool)
 
     args = parser.parse_args()
     shf(**vars(args))