Update src/seismic_hazard_forecasting.py

when the 4 values specifying the lat and lon range of the area of interest (AOI) are provided, only do forecasting for grid points within the AOI

gsim/lasocki_2013.py

@@ -0,0 +1,142 @@
+# -*- coding: utf-8 -*-
+# vim: tabstop=4 shiftwidth=4 softtabstop=4
+#
+# Copyright (C) 2013-2023 GEM Foundation
+#
+# OpenQuake is free software: you can redistribute it and/or modify it
+# under the terms of the GNU Affero General Public License as published
+# by the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# OpenQuake is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU Affero General Public License for more details.
+#
+# You should have received a copy of the GNU Affero General Public License
+# along with OpenQuake. If not, see <http://www.gnu.org/licenses/>.
+
+"""
+Module exports :class:`Lasocki2013`.
+"""
+import numpy as np
+from scipy.constants import g
+
+from openquake.hazardlib.gsim.base import GMPE, CoeffsTable
+from openquake.hazardlib import const
+from openquake.hazardlib.imt import PGA, PGV, SA
+
+def get_magnitude_energy(mag):
+    """
+    Converts magnitude to energy-based magnitude term.
+    """
+    return 1.15 + 1.96 * mag
+
+def get_distance_term(coeffs, repi):
+    """
+    Computes the distance term using the given GMPE coefficients.
+    """
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    return np.log10(R_h)
+
+def get_standard_deviation(coeffs, coeffs_cov, magE, repi):
+    """
+    Computes the standard deviation term.
+    """
+    Cb = np.array(list(coeffs_cov)).reshape(3,3)
+    R_h = np.sqrt(repi ** 2 + coeffs["c7"] ** 2)
+    X0 = np.array([1, magE[0], np.log10(R_h[0]**2)])
+    variance_term = np.sqrt(X0 @ Cb @ X0 + coeffs["sigma"]**2)
+    return variance_term
+
+class Lasocki2013(GMPE):
+    """
+    Implement equation developed by Lasocki in "REPORT ON THE ATTENUATION 
+    RELATIONS OF PEAK GROUND ACCELERATION AND SPECTRAL ORDINATES OF GROUND 
+    MOTION FOR MINING-INDUCED SEISMIC EVENTS IN THE REGION OF THE ŻELAZNY 
+    MOST REPOSITORY", 2009.
+    Equation coefficients provided for the random horizontal component
+    """
+    #: Supported tectonic region type is induced, given
+    #: that the equations have been derived for the LGCD mining area
+    DEFINED_FOR_TECTONIC_REGION_TYPE = const.TRT.INDUCED
+
+    #: Supported intensity measure types are spectral acceleration,
+    #: and peak ground acceleration
+    DEFINED_FOR_INTENSITY_MEASURE_TYPES = {PGA, SA}
+
+    #: Supported intensity measure component is random horizontal
+    #: :attr:`~openquake.hazardlib.const.IMC.RANDOM_HORIZONTAL`,
+    DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.RANDOM_HORIZONTAL
+    # DEFINED_FOR_INTENSITY_MEASURE_COMPONENT = const.IMC.VERTICAL
+
+    #: Supported standard deviation type is total
+    DEFINED_FOR_STANDARD_DEVIATION_TYPES = {const.StdDev.TOTAL}
+
+    #: site params are not required
+    REQUIRES_SITES_PARAMETERS = set()
+
+    #: Required rupture parameter is magnitude
+    REQUIRES_RUPTURE_PARAMETERS = {'mag'}
+
+    #: Required distance measure is epicentral distance
+    #: see paragraph 'Predictor Variables', page 6.
+    REQUIRES_DISTANCES = {'repi'}
+
+    def compute(self, ctx: np.recarray, imts, mean, sig, tau, phi):
+        """
+        Computes mean ground motion values and standard deviations for Lasocki (2013) GMPE.
+        
+        Parameters:
+            repi (float): Epicentral distance (m)
+            mag (float): Earthquake magnitude
+            imts (list of str): List of intensity measure types (e.g., 'PHA', 'PVA')
+            mean (np.array): Array to store computed mean values
+            sig (np.array): Array to store computed standard deviations
+        """
+        # Loop through each IMT and compute values
+        for i, imt in enumerate(imts):            
+            C = self.COEFFS[imt]
+            C_cov = self.COEFFS_COV[imt]
+            mag = ctx.mag
+            repi = ctx.repi*1000.0 # Convert distance from km to m
+            # Compute magnitude energy term
+            magE = get_magnitude_energy(mag)
+
+            # Compute GMPE terms
+            mag_term = C['c1'] + C['c2'] * magE
+            dist_term = C['c5'] * get_distance_term(C, repi)
+            
+            # Compute mean ground motion
+            imean = mag_term + dist_term
+            mean_value = np.log((10 ** imean) / g)  # Convert to natural log scale and divide by g
+            mean[i] = mean_value
+
+            # Compute standard deviation
+            sigma = get_standard_deviation(C, C_cov, magE, repi)
+            sig[i] = sigma
+
+    #: coefficient table provided by report
+    COEFFS = CoeffsTable(sa_damping=5, table="""\
+    IMT  c1     c2      c5      c7    sigma     N        
+    pga  1.25   0.31    -1.34   558   0.196     1196
+    0.6  -3.86  0.55    -0.65   183   0.242     1194
+    1.0  -3.94  0.62    -0.66   308   0.262     1197
+    2.0  -1.14  0.47    -1.02   741   0.235     1195
+    5.0  0.99   0.31    -1.15   690   0.234     1206
+    10.0 3.06   0.27    -1.65   906   0.203     1192
+    20.0 2.62   0.27    -1.60   435   0.196     1191
+    50.0 2.09   0.27    -1.48   375   0.204     1191
+    """)
+    
+    COEFFS_COV = CoeffsTable(sa_damping=5, table="""\
+    IMT  Cb00      Cb01      Cb02      Cb10      Cb11      Cb12      Cb20      Cb21      Cb22        
+    pga  0.005586  -0.000376 -0.000752 -0.000376 0.000103  -0.000111 -0.000752 -0.000111 0.000440
+    0.6  0.007509  -0.000662 -0.000688 -0.000662 0.000161  -0.000154 -0.000688 -0.000154 0.000516
+    1.0  0.009119  -0.00075  -0.000948 -0.00075  0.000189  -0.000187 -0.000948 -0.000187 0.000657                 
+    2.0  0.008563  -0.000514 -0.001282 -0.000514 0.000147  -0.000164 -0.001282 -0.000164 0.000696         
+    5.0  0.008283  -0.000516 -0.001202 -0.000516 0.000145  -0.000161 -0.001202 -0.000161 0.000668
+   10.0  0.006904  -0.00036  -0.001145 -0.00036  0.000108  -0.000126 -0.001145 -0.000126 0.000578
+   20.0  0.005389  -0.000396 -0.000658 -0.000396 0.000104  -0.000107 -0.000658 -0.000107 0.000408
+   50.0  0.005874  -0.000449 -0.000678 -0.000449 0.000114  -0.000114 -0.000678 -0.000114 0.000428
+    """)

src/seismic_hazard_forecasting.py

@@ -52,7 +52,6 @@ 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
@@ -69,6 +68,7 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     from matplotlib.contour import ContourSet
     import xml.etree.ElementTree as ET
     import json
+    import multiprocessing as mp
 
     logger = getDefaultLogger('igfash')
 
@@ -88,11 +88,12 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
     else:
         logger.setLevel(logging.INFO)
 
-    # temporary hard-coded configuration
+    exclude_low_fxy = False # skip low probability areas of the map
-    # 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
 
+    AOI_lat = np.array([51.48, 51.54]) # temporary hard-coding to area of Zelazny Most. To be replaced with user-defined lat and lon range 
+    AOI_lon = np.array([16.15, 16.24])
+    
     # log user selections
     logger.debug(f"User input files\n Catalog: {catalog_file}\n Mc: {mc_file}\n Mag_PDF: {pdf_file}\n Mag: {m_file}")
     logger.debug(
@@ -125,10 +126,6 @@ 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
@@ -221,11 +218,30 @@ verbose: {verbose}")
         utm_zone_letter = u[3]
         logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}")
 
-        # define corners of grid based on global dataset
+        if (None not in AOI_lat) and (None not in AOI_lon):
-        x_min = x.min()
+            use_AOI = True
-        y_min = y.min()
+            
-        x_max = x.max()
+            #convert AOI to UTM
-        y_max = y.max()
+            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
+            
+        else:
+            use_AOI = False
+        
+            # define corners of grid based on global dataset
+            x_min = x.min()
+            y_min = y.min()
+            x_max = x.max()
+            y_max = y.max()
 
         grid_x_max = int(ceil(x_max / grid_dim) * grid_dim)
         grid_x_min = int(floor(x_min / grid_dim) * grid_dim)
@@ -258,7 +274,7 @@ verbose: {verbose}")
         # %% compute KDE and extract PDF
         start = timer()
 
-        if xy_win_method == "TW":
+        if xy_win_method:
             logger.info("Time weighting function selected")
 
             x_weights = np.linspace(0, 15, len(t_windowed))
@@ -319,7 +335,7 @@ verbose: {verbose}")
 
     # run activity rate modeling
     lambdas = [None]
-    if custom_rate != None and forecast_select:
+    if custom_rate != None and forecast_select and not rate_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')
@@ -366,7 +382,9 @@ verbose: {verbose}")
 
         # 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')
+        
         # print("Forecasted activity rates: ", lambdas, "events per", time_unit[:-1])
         logger.info(f"Forecasted activity rates: {lambdas} events per {time_unit} with percentages {lambdas_perc}")
         np.savetxt('activity_rate.csv', np.vstack((lambdas, lambdas_perc)).T, header="lambda, percentage",
@@ -443,42 +461,60 @@ verbose: {verbose}")
         else:
             indices = range(0, len(distances))
 
+        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)
+        else:
+            indices_filtered = indices
+        
         fr = fxy.flatten()
 
         # For each receiver compute estimated ground motion values
-        logger.info(f"Estimating ground motion intensity at {len(indices)} grid points...")
+        logger.info(f"Estimating ground motion intensity at {len(indices_filtered)} grid points...")
 
         PGA = np.zeros(shape=(nx * ny))
 
         start = timer()
 
-        use_pp = False
+        use_pp = True
 
         if use_pp:         # use dask parallel computing
-            pbar = ProgressBar()
+            mp.set_start_method("fork", force=True)
-            pbar.register()
+            iter = indices_filtered
-            # 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=2.0, rx_label=i,
+                                                            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]
                 iml = dask.compute(*imls)
                 iml_grid_raw.append(list(iml))
 
         else:
             iml_grid_raw = []
-            iter = indices
+            iter = indices_filtered
             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 = 2.0, rx_label = i, rtol = 0.1, use_cython=True)
+                                                    IMT_max = IMT_max, 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)
@@ -495,7 +531,15 @@ verbose: {verbose}")
         iml_grid = [[] for _ in range(len(products))]  # final ground motion grids
         iml_grid_prep = iml_grid.copy()  # temp ground motion grids
 
-        if 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)):
@@ -513,17 +557,20 @@ 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
+            # upscale the grid, trim, and interpolate if there are at least 10 grid values with range greater than 0.1
-            up_factor = 4
+            if np.count_nonzero(iml_grid_tmp) >= 10 and vmax-vmin > 0.1:
-            iml_grid_hd = resize(iml_grid_tmp, (up_factor * len(iml_grid_tmp), up_factor * len(iml_grid_tmp)),
+                up_factor = 4
-                                 mode='reflect', anti_aliasing=False)
+                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
 
-            # trim edges so the grid is not so blocky
+            #vmin_hd = min(x for x in iml_grid_hd.flatten() if not isnan(x))
-            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
@@ -536,7 +583,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='bilinear')
             ax.axis('off')
 
             # Save the figure