Merge pull request #10410 from gem/mmi_values

Storing `mmi_values` in pre_execute

openquake/calculators/base.py

@@ -728,6 +728,11 @@ def pre_execute(self):
             with self.monitor('importing inputs', measuremem=True):
                 self.read_inputs()
                 self.save_crmodel()
+            if oq.impact and 'mmi' in oq.inputs:
+                logging.info('Computing MMI-aggregated values')
+                if mmi_values := self.assetcol.get_mmi_values(
+                        oq.aggregate_by, oq.inputs['mmi']):
+                    self.datastore['mmi_values'] = mmi_values
 
     def pre_execute_from_parent(self):
         """

openquake/commonlib/oqvalidation.py

@@ -1886,7 +1886,7 @@ def ruptures_hdf5(self):
     @property
     def impact(self):
         """
-        Return True if we are in Aristotle mode, i.e. there is an HDF5
+        Return True if we are in OQImpact mode, i.e. there is an HDF5
         exposure with a known structure
         """
         exposures = self.inputs.get('exposure', [])
@@ -1895,7 +1895,8 @@ def impact(self):
             if not self.quantiles:
                 self.quantiles = [0.05, 0.95]
             if not self.aggregate_by:
-                self.aggregate_by = [['ID_1'], ['OCCUPANCY']]
+                # self.aggregate_by = [['ID_1'], ['OCCUPANCY']]
+                self.aggregate_by = [['ID_1']]
         return yes
 
     @property

openquake/engine/tests/impact_test.py

@@ -20,8 +20,6 @@
 import pathlib
 import unittest
 import pytest
-import fiona
-from shapely.geometry import shape
 from openquake.calculators.checkers import check
 from openquake.calculators.export import export
 
@@ -56,17 +54,9 @@ def test_impact(n):
 
 
 def test_impact5():
-    # NB: expecting exposure in oq-engine and not in mosaic_dir!
+    # this is a case where there are no assets inside the MMI multipolygons
     if not os.path.exists(expo := cd.parent.parent.parent / 'exposure.hdf5'):
         raise unittest.SkipTest(f'Missing {expo}')
 
     # importing the exposure around Nepal and aggregating it
-    calc = check(cd / 'impact5/job.ini')
-    agg_values = calc.assetcol.get_agg_values
-
-    # this is a case where there are no assets inside the MMI multipolygons
-    shapes = calc.oqparam.inputs['mmi']
-    with fiona.open(f'zip://{shapes}!mi.shp') as f:
-        for feat in f:
-            values = agg_values([['ID_1']], shape(feat.geometry))
-            assert values['number'].sum() == 0
+    check(cd / 'impact5/job.ini')

openquake/risklib/asset.py

@@ -24,8 +24,8 @@
 
 import numpy
 import pandas
-from shapely import contains_xy
-
+import fiona
+from shapely import geometry, contains_xy
 
 from openquake.baselib import hdf5, general, config
 from openquake.baselib.node import Node, context
@@ -48,6 +48,19 @@
               'business_interruption'}
 
 
+def to_mmi(value, MMIs=('I', 'II', 'III', 'IV', 'V', 'VI', 'VII',
+                        'VIII', 'IX', 'X')):
+    """
+    :param value: float in the range 1..10
+    :returns: string "I" .. "X" representing a MMI
+    """
+    if value >= 10.5:
+        raise ValueError(f'{value} is too large to be an MMI')
+    elif value < 0.5:
+        raise ValueError(f'{value} is too small to be an MMI')
+    return MMIs[round(value) - 1]
+
+
 def add_dupl_fields(df, oqfields):
     """
     Add duplicated fields to the DataFrame, if any.
@@ -476,6 +489,25 @@ def get_agg_values(self, aggregate_by, geometry=None):
             agg_values[K] = tuple(dataf[vfields].sum())
         return agg_values
 
+    def get_mmi_values(self, aggregate_by, mmi_file):
+        """
+        :param aggregate_by:
+            a list of lists of tag names (i.e. [['NAME_1']])
+        :param mmi_file:
+            shapefile containing MMI geometries and values
+        :returns:
+            a dictionary MMI -> array with the value fields
+        """
+        out = {}
+        with fiona.open(f'zip://{mmi_file}!mi.shp') as f:
+            for feat in f:
+                geom = geometry.shape(feat.geometry)
+                mmi = to_mmi(feat.properties['PARAMVALUE'])
+                values = self.get_agg_values(aggregate_by, geom)
+                if values['number'].any():
+                    out[mmi] = values
+        return out
+
     def build_aggids(self, aggregate_by):
         """
         :param aggregate_by: list of Ag lists of strings