Merge pull request #10426 from gem/mmi_values

Fixed `get_mmi_values`

openquake/calculators/base.py

@@ -730,13 +730,10 @@ def pre_execute(self):
                 self.read_inputs()
                 self.save_crmodel()
             if oq.impact and 'mmi' in oq.inputs:
-                if len(self.assetcol) < 10_000:
-                    logging.info('Computing MMI-aggregated values')
-                    if mmi_values := self.assetcol.get_mmi_values(
-                            oq.aggregate_by, oq.inputs['mmi']):
-                        self.datastore['mmi_tags'] = mmi_values
-                else:
-                    logging.info('Too many assets to aggregate by MMI!')
+                logging.info('Computing MMI-aggregated values')
+                if mmi_values := self.assetcol.get_mmi_values(
+                        oq.aggregate_by, oq.inputs['mmi']):
+                    self.datastore['mmi_tags'] = mmi_values
 
     def pre_execute_from_parent(self):
         """

openquake/risklib/asset.py

@@ -25,7 +25,7 @@
 import numpy
 import pandas
 import fiona
-from shapely import geometry, contains_xy
+from shapely import geometry, prepare, contains_xy
 
 from openquake.baselib import hdf5, general, config
 from openquake.baselib.node import Node, context
@@ -461,8 +461,7 @@ def get_agg_values(self, aggregate_by, geometry=None):
             self.tagcol.get_aggkey(aggregate_by))}
         K = len(aggkey)
         if geometry:
-            lonlats = numpy.column_stack([self['lon'], self['lat']])
-            array = self.array[contains_xy(geometry, lonlats)]
+            array = self.array[contains_xy(geometry, self['lon'], self['lat'])]
         else:
             array = self.array
 
@@ -502,12 +501,30 @@ def get_mmi_values(self, aggregate_by, mmi_file):
         with fiona.open(f'zip://{mmi_file}!mi.shp') as f:
             for feat in f:
                 geom = geometry.shape(feat.geometry)
+                prepare(geom)  # MANDATORY: gives an incredible speedup!
                 mmi = to_mmi(feat.properties['PARAMVALUE'])
                 values = self.get_agg_values(aggregate_by, geom)
                 if values['number'].any():
-                    out[mmi] = values
+                    if mmi not in out:
+                        out[mmi] = values
+                    else:
+                        for lt in values.dtype.names:
+                            out[mmi][lt] += values[lt]
         return out
 
+    # not used yet
+    def agg_by_site(self):
+        """
+        :returns: an array of aggregated values indexed by site ID
+        """
+        N = self['site_id'].max() + 1
+        vfields = self.fields + self.occfields
+        agg_values = numpy.zeros(N, [(f, F32) for f in vfields])
+        for vf in vfields:
+            arr = self['value-' + vf if vf in self.fields else vf]
+            agg_values[vf] = general.fast_agg(self['site_id'], arr)
+        return agg_values
+
     def build_aggids(self, aggregate_by):
         """
         :param aggregate_by: list of Ag lists of strings