Added test for consequence=losses for liquefaction and landslide

debian/changelog

@@ -1,4 +1,5 @@
   [Michele Simionato]
+  * Added support for consequence=losses for liquefaction and landslides
   * Added a check for missing secondary perils
   * Added loss types liquefaction and landslide
   * Removed support for XML consequences, after 3 years of deprecation

openquake/calculators/event_based_damage.py

@@ -18,6 +18,7 @@
 
 import os.path
 import logging
+from dataclasses import dataclass
 import numpy
 import pandas
 
@@ -36,6 +37,19 @@
 U32 = numpy.uint32
 F32 = numpy.float32
 
+@dataclass
+class Dparam:
+    """
+    Parameters for a damage calculation
+    """
+    eids: U32
+    aggids: U16
+    rlzs: U32
+    ci: dict
+    D: int
+    Dc: int
+    rng: scientific.MultiEventRNG
+
 
 def zero_dmgcsq(A, R, crmodel):
     """
@@ -66,71 +80,56 @@ def damage_from_gmfs(gmfslices, oqparam, dstore, monitor):
     return event_based_damage(df, oqparam, dstore, monitor)
 
 
-def _update(asset_df, gmf_df, aggids, allrlzs, sec_sims,
-            crmodel, ci, R, Dc, dmgcsq, dddict):
+def _gen_d4(asset_df, gmf_df, crmodel, dparam):
+    # yields (aids, d4) triples
     oq = crmodel.oqparam
-    loss_types = oq.loss_types
-    eids = gmf_df.eid.to_numpy()
-    if R > 1:
-        rlzs = allrlzs[eids]
-    if sec_sims or not oq.float_dmg_dist:
-        rng = scientific.MultiEventRNG(
-            oq.master_seed, numpy.unique(eids))
+    sec_sims = oq.secondary_simulations.items()
     for prob_field, num_sims in sec_sims:
         probs = gmf_df[prob_field].to_numpy()   # LiqProb
         if not oq.float_dmg_dist:
-            dprobs = rng.boolean_dist(probs, num_sims).mean(axis=1)
+            dprobs = dparam.rng.boolean_dist(probs, num_sims).mean(axis=1)
+    E = len(dparam.eids)
+    L = len(oq.loss_types)
     for taxo, adf in asset_df.groupby('taxonomy'):
         out = crmodel.get_output(adf, gmf_df)
         aids = adf.index.to_numpy()
+        A = len(aids)
         assets = adf.to_records()
         if oq.float_dmg_dist:
             number = assets['value-number']
         else:
             number = U32(assets['value-number'])
-        for lti, lt in enumerate(loss_types):
+        d4 = numpy.zeros((L, A, E, dparam.Dc), F32)
+        D = dparam.D
+        for lti, lt in enumerate(oq.loss_types):
             fractions = out[lt]
-            Asid, E, D = fractions.shape
-            assert len(eids) == E
-            d3 = numpy.zeros((Asid, E, Dc), F32)
             if oq.float_dmg_dist:
-                d3[:, :, :D] = fractions
-                for a in range(Asid):
-                    d3[a] *= number[a]
+                d4[lti, :, :, :D] = fractions
+                for a in range(A):
+                    d4[lti, a] *= number[a]
             else:
                 # this is a performance distaster; for instance
                 # the Messina test in oq-risk-tests becomes 12x
                 # slower even if it has only 25_736 assets
-                d3[:, :, :D] = rng.discrete_dmg_dist(
-                    eids, fractions, number)
+                d4[lti, :, :, :D] = dparam.rng.discrete_dmg_dist(
+                    dparam.eids, fractions, number)
 
             # secondary perils and consequences
             for a, asset in enumerate(assets):
                 if sec_sims:
                     for d in range(1, D):
                         # doing the mean on the secondary simulations
                         if oq.float_dmg_dist:
-                            d3[a, :, d] *= probs
+                            d4[lti, a, :, d] *= probs
                         else:
-                            d3[a, :, d] *= dprobs
+                            d4[lti, a, :, d] *= dprobs
 
                 csq = crmodel.compute_csq(
-                    asset, d3[a, :, :D] / number[a], lt,
+                    asset, d4[lti, a, :, :D] / number[a], lt,
                     oq.time_event)
                 for name, values in csq.items():
-                    d3[a, :, ci[name]] = values
-            if R == 1:
-                dmgcsq[aids, 0, lti] += d3.sum(axis=1)
-            else:
-                for e, rlz in enumerate(rlzs):
-                    dmgcsq[aids, rlz, lti] += d3[:, e]
-            tot = d3.sum(axis=0)  # sum on the assets
-            for e, eid in enumerate(eids):
-                dddict[eid, oq.K][lti] += tot[e]
-                if oq.K:
-                    for kids in aggids:
-                        for a, aid in enumerate(aids):
-                            dddict[eid, kids[aid]][lti] += d3[a, e]
+                    d4[lti, a, :, dparam.ci[name]] = values
+        yield aids, d4  # d4 has shape (L, A, E, Dc)
 
 
 def event_based_damage(df, oq, dstore, monitor):
@@ -144,27 +143,21 @@ def event_based_damage(df, oq, dstore, monitor):
     mon_risk = monitor('computing risk', measuremem=False)
     with monitor('reading gmf_data'):
         if oq.parentdir:
-            dstore = datastore.read(
-                oq.hdf5path, parentdir=oq.parentdir)
+            dstore = datastore.read(oq.hdf5path, parentdir=oq.parentdir)
         else:
             dstore.open('r')
         assetcol = dstore['assetcol']
-        if oq.K:
-            # TODO: move this in the controller!
-            aggids, _ = assetcol.build_aggids(
-                oq.aggregate_by, oq.max_aggregations)
-        else:
-            aggids = numpy.zeros(len(assetcol), U16)
         crmodel = monitor.read('crmodel')
-    sec_sims = oq.secondary_simulations.items()
+        aggids = monitor.read('aggids')
     dmg_csq = crmodel.get_dmg_csq()
     ci = {dc: i + 1 for i, dc in enumerate(dmg_csq)}
     dmgcsq = zero_dmgcsq(len(assetcol), oq.R, crmodel)
     _A, R, L, Dc = dmgcsq.shape
+    D = Dc - len(crmodel.get_consequences())
     if R > 1:
         allrlzs = dstore['events']['rlz_id']
     else:
-        allrlzs = [0]
+        allrlzs = U32([0])
     assert len(oq.loss_types) == L
     with mon_risk:
         dddict = general.AccumDict(accum=numpy.zeros((L, Dc), F32))  # eid, kid
@@ -173,8 +166,33 @@ def event_based_damage(df, oq, dstore, monitor):
             gmf_df = df[df.sid == sid]
             if len(gmf_df) == 0:
                 continue
-            _update(asset_df, gmf_df, aggids, allrlzs, sec_sims,
-                    crmodel, ci, R, Dc, dmgcsq, dddict)
+            oq = crmodel.oqparam
+            eids = gmf_df.eid.to_numpy()
+            if R > 1:
+                rlzs = allrlzs[eids]
+            else:
+                rlzs = allrlzs
+            if oq.secondary_simulations or not oq.float_dmg_dist:
+                rng = scientific.MultiEventRNG(
+                    oq.master_seed, numpy.unique(eids))
+            else:
+                rng = None
+            dparam = Dparam(eids, aggids, rlzs, ci, D, Dc, rng)
+            for aids, d4 in _gen_d4(asset_df, gmf_df, crmodel, dparam):
+                for lti, d3 in enumerate(d4):
+                    if R == 1:
+                        dmgcsq[aids, 0, lti] += d3.sum(axis=1)
+                    else:
+                        for e, rlz in enumerate(dparam.rlzs):
+                            dmgcsq[aids, rlz, lti] += d3[:, e]
+                    tot = d3.sum(axis=0)  # sum on the assets
+                    for e, eid in enumerate(eids):
+                        dddict[eid, oq.K][lti] += tot[e]
+                        if oq.K:
+                            for kids in dparam.aggids:
+                                for a, aid in enumerate(aids):
+                                    dddict[eid, kids[aid]][lti] += d3[a, e]
+
     return _dframe(dddict, ci, oq.loss_types), dmgcsq
 
 
@@ -226,8 +244,14 @@ def execute(self):
         if oq.investigation_time:  # event based
             self.builder = get_loss_builder(self.datastore, oq)  # check
         self.dmgcsq = zero_dmgcsq(len(self.assetcol), self.R, self.crmodel)
-        smap = calc.starmap_from_gmfs(damage_from_gmfs, oq, self.datastore,
-                                      self._monitor)
+        if oq.K:
+            aggids, _ = self.assetcol.build_aggids(
+                oq.aggregate_by, oq.max_aggregations)
+        else:
+            aggids = 0
+        smap = calc.starmap_from_gmfs(damage_from_gmfs, oq,
+                                      self.datastore, self._monitor)
+        smap.monitor.save('aggids', aggids)
         smap.monitor.save('assets', self.assetcol.to_dframe('id'))
         smap.monitor.save('crmodel', self.crmodel)
         return smap.reduce(self.combine)

openquake/calculators/event_based_risk.py

@@ -247,7 +247,7 @@ def gen_outputs(df, crmodel, rng, monitor):
             yield out
 
 
-def check_tot_loss_unit_consistency(units, total_losses, loss_types):
+def _tot_loss_unit_consistency(units, total_losses, loss_types):
     total_losses_units = set()
     for separate_lt in total_losses.split('+'):
         assert separate_lt in loss_types
@@ -279,10 +279,12 @@ def set_oqparam(oq, assetcol, dstore):
                 partial(insurance_losses, policy_df=policy_df))
 
     ideduc = assetcol['ideductible'].any()
-    if oq.total_losses:
-        units = dstore['exposure'].cost_calculator.get_units(oq.loss_types)
-        check_tot_loss_unit_consistency(
-            units.split(), oq.total_losses, oq.loss_types)
+    cc = dstore['exposure'].cost_calculator
+    if oq.total_losses and cc.cost_types:
+        # cc.cost_types is empty in scenario_damage/case_21 (consequences)
+        units = cc.get_units(oq.total_loss_types)
+        _tot_loss_unit_consistency(
+            units.split(), oq.total_losses, oq.total_loss_types)
         sec_losses.append(
             partial(total_losses, kind=oq.total_losses, ideduc=ideduc))
     elif ideduc:

openquake/calculators/tests/scenario_damage_test.py

@@ -24,7 +24,7 @@
 from openquake.qa_tests_data.scenario_damage import (
     case_1, case_1c, case_2, case_3, case_4, case_4b, case_5, case_5a,
     case_6, case_7, case_8, case_9, case_10, case_11, case_12, case_13,
-    case_14, case_16, case_17, case_18, case_19, case_20, case_21)
+    case_14, case_16, case_17, case_18, case_19, case_20, case_21, case_22)
 from openquake.calculators.tests import CalculatorTestCase, strip_calc_id
 from openquake.calculators.extract import extract
 from openquake.calculators.export import export
@@ -298,12 +298,18 @@ def test_case_20(self):
         self.assertEqualFiles('expected/aggrisk.csv', fname)
 
     def test_case_21(self):
-        # infrastructure risk for structural and liquefaction loss types
+        # infrastructure risk for structural, liquefaction and landslides
         out = self.run_calc(case_21.__file__, 'job.ini', exports='csv')
         [agg_csv, aggparent_csv] = out[('aggrisk', 'csv')]
         self.assertEqualFiles('expected/aggrisk.csv', agg_csv)
         self.assertEqualFiles('expected/aggrisk-parent.csv', aggparent_csv)
 
+    def test_case_22(self):
+        # losses with liquefaction and landslides
+        out = self.run_calc(case_22.__file__, 'job.ini', exports='csv')
+        [agg_csv] = out[('aggrisk', 'csv')]
+        self.assertEqualFiles('expected/aggrisk.csv', agg_csv)
+
 
 def losses(aid, alt):
     E = len(alt.event_id.unique())

openquake/commonlib/oqvalidation.py

@@ -1745,6 +1745,18 @@ def ext_loss_types(self):
             etypes = self.loss_types + itypes
         return etypes
 
+    @property
+    def total_loss_types(self):
+        """
+        :returns: the loss types in total_losses or the single loss type
+        """
+        if self.total_losses:
+            return self.total_losses.split('+')
+        elif len(self.loss_types) == 1:
+            return self.loss_types
+        else:
+            self.raise_invalid('please specify total_losses')
+
     def loss_dt(self, dtype=F64):
         """
         :returns: a composite dtype based on the loss types including occupants

openquake/qa_tests_data/scenario_damage/case_21/job.ini

@@ -30,6 +30,7 @@ liquefaction_fragility_file = fragility_model_liquefaction.xml
 landslide_fragility_file = fragility_model_landslide.xml
 
 [consequence]
+total_losses = structural
 consequence_file = {'taxonomy': 'consequence_multiple_loss_types.csv'}
 
 [risk_calculation]

openquake/qa_tests_data/scenario_damage/case_22/consequences.csv

@@ -0,0 +1,7 @@
+taxonomy,consequence,loss_type,slight,moderate,extreme,complete
+Concrete,losses,structural,0.04,0.31,0.6,1
+Wood,losses,structural,0.04,0.31,0.6,1
+Concrete,losses,liquefaction,0,0,0,1
+Wood,losses,liquefaction,0,0,0,1
+Concrete,losses,landslide,0,0,0,1
+Wood,losses,landslide,0,0,0,1

openquake/qa_tests_data/scenario_damage/case_22/expected/aggrisk.csv

@@ -0,0 +1,5 @@
+#,,,,,,,"generated_by='OpenQuake engine 3.22.0-gitff750fe4b4', start_date='2024-10-21T11:16:00', checksum=2624933730, investigation_time=None, risk_investigation_time=None"
+loss_type,no_damage,slight,moderate,extreme,complete,losses_value,losses_ratio
+structural,1.06701E+01,1.41803E+00,1.39786E+00,9.49061E-01,5.64993E-01,1.84217E+04,1.08299E-01
+liquefaction,1.26268E+01,0.00000E+00,0.00000E+00,0.00000E+00,2.37321E+00,0.00000E+00,nan
+landslide,9.00000E+00,0.00000E+00,0.00000E+00,0.00000E+00,6.00000E+00,0.00000E+00,nan

openquake/qa_tests_data/scenario_damage/case_22/exposure_model.csv

@@ -0,0 +1,16 @@
+id,lon,lat,number,structural,night,taxonomy,NAME_1
+a1,83.31382,29.46117,1,11340,5,Wood1,a
+a2,83.31382,29.23617,1,11340,5,Wood1,a
+a3,83.53882,29.08617,1,11340,5,Wood1,a
+a4,80.68882,28.93617,1,11340,5,Wood1,a
+a5,83.53882,29.01117,1,11340,5,Wood1,a
+a6,81.13882,28.78617,1,11340,5,Wood1,a
+a7,83.98882,28.48617,1,11340,5,Wood1,a
+a8,83.23882,29.38617,1,11340,5,Concrete1,a
+a9,83.01382,29.08617,1,11340,5,Concrete1,a
+a10,83.31382,28.71117,1,11340,5,Concrete1,a
+a11,86.91382,27.73617,1,11340,5,Concrete1,a
+a12,83.16382,29.31117,1,11340,5,Concrete1,a
+a13,80.61382,28.93617,1,11340,5,Concrete1,a
+a14,83.91382,29.01117,1,11340,5,Concrete1,a
+a15,82.03882,30.28617,1,11340,5,Concrete1,a

openquake/qa_tests_data/scenario_damage/case_22/exposure_model.xml

@@ -0,0 +1,27 @@
+<?xml version="1.0" encoding="utf-8"?>
+<nrml
+xmlns="http://openquake.org/xmlns/nrml/0.5"
+xmlns:gml="http://www.opengis.net/gml"
+>
+    <exposureModel
+    category="buildings"
+    id="ex1"
+    taxonomySource="GEM taxonomy"
+    >
+        <description />
+        <conversions>
+            <costTypes>
+                <costType name="structural" type="per_asset" unit="USD"/>
+            </costTypes>
+        </conversions>
+        <occupancyPeriods>
+            night
+        </occupancyPeriods>
+        <tagNames>
+            NAME_1
+        </tagNames>
+        <assets>
+            exposure_model.csv
+        </assets>
+    </exposureModel>
+</nrml>

openquake/qa_tests_data/scenario_damage/case_22/fault_rupture.xml

@@ -0,0 +1,19 @@
+<?xml version="1.0" encoding="utf-8"?>
+<nrml xmlns:gml="http://www.opengis.net/gml" xmlns="http://openquake.org/xmlns/nrml/0.4">
+  <simpleFaultRupture>
+    <magnitude>7.0</magnitude>
+    <rake>90</rake>
+    <hypocenter lat="27.6" lon="84.4" depth="30.0"/>
+    <simpleFaultGeometry>
+      <gml:LineString>
+        <gml:posList>
+          85.0 27.3
+          83.8 27.8
+        </gml:posList>
+      </gml:LineString>
+      <dip>30.0</dip>
+      <upperSeismoDepth>20.0</upperSeismoDepth>
+      <lowerSeismoDepth>50.0</lowerSeismoDepth>
+    </simpleFaultGeometry>
+  </simpleFaultRupture>
+</nrml>

openquake/qa_tests_data/scenario_damage/case_22/fragility_model_landslide.xml

@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="UTF-8"?>
+<nrml xmlns="http://openquake.org/xmlns/nrml/0.5">
+	<fragilityModel id="1" assetCategory="RDN" lossCategory="landslide">
+		<description>Road Network</description>
+		<limitStates>slight moderate extreme complete</limitStates>
+		<fragilityFunction id="Wood" format="discrete">
+			<imls imt="DispProb" noDamageLimit="0">0 0.3 0.6 1</imls>
+			<poes ls="slight">0 0 0 1</poes>
+			<poes ls="moderate">0 0 0 1</poes>
+			<poes ls="extreme">0 0 0 1</poes>
+			<poes ls="complete">0 0 0 1</poes>
+		</fragilityFunction>
+		<fragilityFunction id="Concrete" format="discrete">
+			<imls imt="DispProb" noDamageLimit="0">0 0.3 0.6 1</imls>
+			<poes ls="slight">0 1 1 1</poes>
+			<poes ls="moderate">0 1 1 1</poes>
+			<poes ls="extreme">0 1 1 1</poes>
+			<poes ls="complete">0 1 1 1</poes>
+		</fragilityFunction>
+	</fragilityModel>
+</nrml>