ISEPOS-2450 Added AOI_extent argument to shf_wrapper

put AOI_extent as a parameter of the main function

src/seismic_hazard_forecasting.py

@@ -3,7 +3,7 @@
 
 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):
+         forecast_select, custom_rate, forecast_len, time_unit, AOI_extent, model, products_string, verbose):
     """
     Python application that reads an earthquake catalog and performs seismic hazard forecasting.
     Arguments:
@@ -33,6 +33,8 @@ def main(catalog_file, mc_file, pdf_file, m_file, m_select, mag_label, mc, m_max
                     forecasting.
         forecast_len: Length of the forecast for seismic hazard assessment.
         time_unit: Times units for the inputs Time Window Duration, Custom Activity Rate, and Forecast Length.
+        AOI_extent: The forecast geographical area of interest specified as a latitude and longitude range in decimal degrees
+                    in the form [lat_min, lat_max, lon_min, lon_max].
         model: Select from the following ground motion models available. Other models in the Openquake library are
                     available but have not yet been tested.
         products_string: The ground motion intensity types to output. Use a space between names to select more than
@@ -88,9 +90,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(AOI_extent[:2]) 
+    AOI_lon = np.array(AOI_extent[2:])
+    
     # 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,11 +220,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)
@@ -439,10 +463,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 +485,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 +504,7 @@ verbose: {verbose}")
 
         else:
             iml_grid_raw = []
-            iter = indices
+            iter = indices_filtered
             for imt in products:
 
                 if imt == "PGV":
@@ -500,7 +533,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)):

src/shf_wrapper.py

@@ -31,6 +31,9 @@ def main(argv):
             return False
         else:
             raise argparse.ArgumentTypeError("Boolean value expected.")
+        
+    def float_or_none(v):
+        return None if v.lower() == "none" else float(v)
 
     parser = argparse.ArgumentParser()
 
@@ -55,6 +58,7 @@ def main(argv):
     parser.add_argument("--time_unit", type=str)
     parser.add_argument("--model", type=str)
     parser.add_argument("--products_string", type=str)
+    parser.add_argument("--AOI_extent", nargs=4, type=float_or_none, default=[None] * 4, required=False)
     parser.add_argument("--verbose", type=str2bool)
 
     args = parser.parse_args()