Update src/seismic_hazard_forecasting.py

src/seismic_hazard_forecasting.py

@@ -4,6 +4,27 @@
 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):
+    #catalog_file = 'LGCD_catalog_2024.mat'
+    mc_file = None
+    pdf_file = None
+    m_file = None
+    m_select = True
+    mag_label = 'Mw'
+    mc = 1.8
+    m_max = 4.5
+    m_kde_method = 'arviz-silverman'
+    xy_select = True
+    grid_dim = 5000
+    xy_win_method = False
+    rate_select = True
+    time_win_duration = 100
+    forecast_select = True
+    custom_rate = None
+    forecast_len = 100
+    time_unit = 'days'
+    model = 'Lasocki2013'
+    products_string = 'PGA'
+    verbose = True
     """
     Python application that reads an earthquake catalog and performs seismic hazard forecasting.
     Arguments:
@@ -52,6 +73,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     from math import ceil, floor, isnan
     import numpy as np
     import dask
+    from dask.diagnostics import ProgressBar  # use Dask progress bar
     import kalepy as kale
     import utm
     from skimage.transform import resize
@@ -88,7 +110,9 @@ 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
+    # temporary hard-coded configuration
+    # exclude_low_fxy = False
+    exclude_low_fxy = True
     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
 
     # log user selections
@@ -123,6 +147,10 @@ verbose: {verbose}")
             logger.info("No magnitude label of catalog specified, therefore try Mw by default")
             mag_label = 'Mw'
 
+        # if cat_label == None:
+        #   print("No magnitude label of catalog specified, therefore try 'Catalog' by default")
+        #   cat_label='Catalog'
+
         time, mag, lat, lon, depth = read_mat_cat(catalog_file, mag_label=mag_label, catalog_label='Catalog')
 
         # check for null magnitude values
@@ -252,7 +280,7 @@ verbose: {verbose}")
         # %% compute KDE and extract PDF
         start = timer()
 
-        if xy_win_method:
+        if xy_win_method == "TW":
             logger.info("Time weighting function selected")
 
             x_weights = np.linspace(0, 15, len(t_windowed))
@@ -313,7 +341,7 @@ verbose: {verbose}")
 
     # run activity rate modeling
     lambdas = [None]
-    if custom_rate != None and forecast_select and not rate_select:
+    if custom_rate != None and forecast_select:
         logger.info(f"Using activity rate specified by user: {custom_rate} per {time_unit}")
         lambdas = np.array([custom_rate], dtype='d')
         lambdas_perc = np.array([1], dtype='d')
@@ -452,19 +480,17 @@ verbose: {verbose}")
 
         if use_pp:         # use dask parallel computing
             mp.set_start_method("fork", force=True)
+            #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}")
 
-                if imt == "PGV":
-                    IMT_max = 200 # search interval max for velocity (cm/s)
-                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,
                                                             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,
+                                                            log_level=logging.DEBUG, imt=imt, IMT_min=0.0, IMT_max=2.0, rx_label=i,
                                                             rtol=0.1, use_cython=True) for i in iter]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
@@ -473,17 +499,11 @@ verbose: {verbose}")
             iml_grid_raw = []
             iter = indices
             for imt in products:
-
-                if imt == "PGV":
-                    IMT_max = 200 # search interval max for velocity (cm/s)
-                else:
-                    IMT_max = 2.0 # search interval max for acceleration (g)
-
                 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 = IMT_max, rx_label = i, rtol = 0.1, use_cython=True)
+                                                    IMT_max = 2.0, rx_label = i, rtol = 0.1, use_cython=True)
                     iml.append(iml_i)
                     logger.info(f"Estimated {imt} at rx {i} is {iml_i}")
                 iml_grid_raw.append(iml)
@@ -518,20 +538,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, 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:
+            # upscale the grid
             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)
-                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
 
-            #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
+            # 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
@@ -544,7 +561,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, interpolation='bilinear')
+            ax.imshow(iml_grid_hd, origin='lower', cmap=cmap, vmin=vmin, vmax=vmax)
             ax.axis('off')
 
             # Save the figure