Using the circulant embedding regular grid in OQImpact

openquake/calculators/base.py

@@ -1209,6 +1209,15 @@ def _read_risk3(self):
                 assoc_dist = (oq.region_grid_spacing * 1.414
                               if oq.region_grid_spacing else ASSOC_DIST)
                 sm = readinput.get_site_model(oq, self.datastore.hdf5)
+                if oq.impact:
+                    # the grid sites have no parameters and req_site_params
+                    # is not known yet: add all the site model ones, both
+                    # to the sitecol and to the complete sitecol
+                    for name in sm.dtype.names:
+                        if name in site.site_param_dt:
+                            dt = site.site_param_dt[name]
+                            self.sitecol.add_col(name, dt)
+                            self.sitecol.complete.add_col(name, dt)
                 if oq.prefer_global_site_params and not numpy.isnan(
                         oq.reference_vs30_value):
                     self.sitecol.set_global_params(oq, oq.req_site_params)
@@ -1221,7 +1230,7 @@ def _read_risk3(self):
                             or oq.siteid else 'strict')
                     self.sitecol.assoc(sm, assoc_dist, mode,
                                        station_sids=station_sids)
-                    if 'station_data' in oq.inputs:
+                    if 'station_data' in oq.inputs or oq.impact:
                         # the complete sitecol is required
                         self.sitecol.complete.assoc(sm, assoc_dist, mode,
                                                     station_sids=station_sids)

openquake/commonlib/oqvalidation.py

@@ -743,6 +743,8 @@
 
 region_grid_spacing:
   Used together with the *region* option to generate the hazard sites.
+  In impact mode it is mandatory and it defines the spacing in km of the
+  regular grid of sites around the rupture (circulant embedding).
   Example: *region_grid_spacing = 10*.
   Default: None
 
@@ -1781,6 +1783,7 @@ def validate(self):
         """
         super().validate()
         self.check_source_model()
+        self.check_impact()
         if isinstance(self.hazard_calculation_id, str):
             path = self.hazard_calculation_id
             if not os.path.isabs(path):
@@ -2710,6 +2713,12 @@ def check_uniform_hazard_spectra(self):
                 'There is a single IMT, the uniform_hazard_spectra plot will '
                 'contain a single point')
 
+    def check_impact(self):
+        # in impact mode the sites are a regular grid with this spacing
+        if self.impact and not self.region_grid_spacing:
+            self.raise_invalid('region_grid_spacing must be set in impact '
+                               'mode')
+
     def check_source_model(self):
         if ('hazard_curves' in self.inputs or 'gmfs' in self.inputs
                 or 'rupture_model' in self.inputs

openquake/commonlib/readinput.py

@@ -69,6 +69,7 @@
     build_planar_rupture_from_dict, get_ruptures, get_ebrupture)
 from openquake.hazardlib.map_array import MapArray
 from openquake.hazardlib.geo.utils import hex6
+from openquake.hazardlib.geo.regular_grid import regular_grid
 from openquake.hazardlib.shakemap.parsers import convert_to_oq_xml
 from openquake.hazardlib.countries import country2code, MODELS, ALIASES
 from openquake.pfd.gsim import get_pfd_gsim_lt
@@ -568,9 +569,12 @@ def get_mesh_exp(oqparam, h5=None):
     """
     exposure = get_exposure(oqparam, h5)
     if oqparam.impact:
+        # the sites are a regular UTM grid enclosing the site model around
+        # the rupture, as required by circulant embedding; the site model
+        # parameters are associated to the grid in base.py
         sm = get_site_model(oqparam, h5)
         mesh = geo.Mesh(sm['lon'], sm['lat'])
-        return mesh, exposure
+        return regular_grid(mesh, oqparam.region_grid_spacing), exposure
     if oqparam.sites:
         mesh = geo.Mesh.from_coords(oqparam.sites)
         return mesh, exposure
@@ -859,15 +863,12 @@ def get_site_collection(oqparam, h5=None):
             req_site_params = set()   # no parameters are required
         else:
             req_site_params = oqparam.req_site_params
-        if h5 and 'site_model' in h5:
+        if oqparam.impact:
+            # the grid is not the site model: the parameters are associated
+            # to the sites in base.py
+            sm = None
+        elif h5 and 'site_model' in h5:
             sm = h5['site_model'][:]
-        elif oqparam.impact and (
-                    not oqparam.infrastructure_connectivity_analysis):
-            # filter the far away sites
-            rup = get_rupture(oqparam)
-            dist = oqparam.maximum_distance('*')(rup.mag)
-            [expo_hdf5] = oqparam.inputs['exposure']
-            sm = get_site_model_around(expo_hdf5, rup, dist)
         elif (not h5 and 'site_model' in oqparam.inputs and
               'exposure' not in oqparam.inputs):
             # tested in test_with_site_model