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

src/seismic_hazard_forecasting.py

@@ -88,9 +88,12 @@ 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
+    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([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(
@@ -215,6 +218,25 @@ verbose: {verbose}")
         utm_zone_letter = u[3]
         logger.debug(f"Latitude / Longitude coordinates correspond to UTM zone {utm_zone_number}{utm_zone_letter}")
 
+        if (None not in AOI_lat) and (None not in AOI_lon):
+            use_AOI = True
+            
+            #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
+            
+        else:
+            use_AOI = False
+        
             # define corners of grid based on global dataset
             x_min = x.min()
             y_min = y.min()
@@ -439,10 +461,19 @@ 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))
 
@@ -452,7 +483,7 @@ verbose: {verbose}")
 
         if use_pp:         # use dask parallel computing
             mp.set_start_method("fork", force=True)
-            iter = indices
+            iter = indices_filtered
             iml_grid_raw = []  # raw ground motion grids
             for imt in products:
                 logger.info(f"Estimating {imt}")
@@ -471,7 +502,7 @@ verbose: {verbose}")
 
         else:
             iml_grid_raw = []
-            iter = indices
+            iter = indices_filtered
             for imt in products:
 
                 if imt == "PGV":
@@ -500,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)):