Merge pull request #948 from CLIMADA-project/bugfix/imp-calc-mdr-zeros

Bug fix in Hazard.get_mdr and impact calculation in case of a single exposure point and MDR values of 0

Author: bguillod

CHANGELOG.md

@@ -21,6 +21,8 @@ Code freeze date: YYYY-MM-DD
 
 ### Fixed
 
+- Avoids a ValueError in the impact calculation for cases with a single exposure point and MDR values of 0, by explicitly removing zeros in `climada.hazard.Hazard.get_mdr` [#933](https://github.com/CLIMADA-project/climada_python/pull/948)
+
 ### Deprecated
 
 ### Removed

climada/engine/test/test_impact_calc.py

@@ -28,8 +28,8 @@
 
 from climada import CONFIG
 from climada.entity.entity_def import Entity
-from climada.entity import Exposures, ImpactFuncSet, ImpactFunc
-from climada.hazard.base import Hazard
+from climada.entity import Exposures, ImpactFuncSet, ImpactFunc, ImpfTropCyclone
+from climada.hazard.base import Hazard, Centroids
 from climada.engine import ImpactCalc, Impact
 from climada.engine.impact_calc import LOGGER as ILOG
 from climada.util.constants import ENT_DEMO_TODAY, DEMO_DIR
@@ -471,6 +471,34 @@ def test_stitch_risk_metrics(self):
         np.testing.assert_array_equal(eai_exp, [2.25, 1.25, 4.5])
         self.assertEqual(aai_agg, 8.0)  # Sum of eai_exp
 
+    def test_single_exp_zero_mdr(self):
+        """Test for case where exposure has a single value and MDR or fraction contains zeros"""
+        centroids = Centroids.from_lat_lon([-26.16], [28.20])
+        haz = Hazard(
+            intensity=sparse.csr_matrix(np.array([[31.5], [19.0]])),
+            event_id=np.arange(2),
+            event_name=[0,1],
+            frequency=np.ones(2) / 2,
+            fraction=sparse.csr_matrix(np.zeros((2,1))),
+            date=np.array([0, 1]),
+            centroids=centroids,
+            haz_type='TC'
+        )
+        exp = Exposures({'value': [1.],
+                            'longitude': 28.22,
+                            'latitude': -26.17,
+                            'impf_TC': 1},
+                            crs="EPSG:4326")
+        imp_evt = 0.00250988804927603
+        aai_agg = imp_evt/2
+        eai_exp = np.array([aai_agg])
+        at_event = np.array([imp_evt, 0])
+        exp.set_geometry_points()
+        impf_tc = ImpfTropCyclone.from_emanuel_usa()
+        impf_set = ImpactFuncSet([impf_tc])
+        impf_set.check()
+        imp = ImpactCalc(exp, impf_set, haz).impact(save_mat=True)
+        check_impact(self, imp, haz, exp, aai_agg, eai_exp, at_event, at_event)
 
 class TestImpactMatrixCalc(unittest.TestCase):
     """Verify the computation of the impact matrix"""

climada/hazard/base.py

@@ -1103,7 +1103,9 @@ def get_mdr(self, cent_idx, impf):
             impf.id)
             mdr_array = impf.calc_mdr(mdr.toarray().ravel()).reshape(mdr.shape)
             mdr = sparse.csr_matrix(mdr_array)
-        return mdr[:, indices]
+        mdr_out = mdr[:, indices]
+        mdr_out.eliminate_zeros()
+        return mdr_out
 
     def get_paa(self, cent_idx, impf):
         """

climada/hazard/test/test_base.py

@@ -1176,6 +1176,14 @@ def test_get_mdr(self):
             true_mdr = np.digitize(haz.intensity[:, idx].toarray(), [0, 1])
             np.testing.assert_array_almost_equal(mdr.toarray(), true_mdr)
 
+        # #case with zeros everywhere
+        cent_idx = np.array([0, 0, 1])
+        impf.mdd=np.array([0,0,0,1])
+        # how many non-zeros values are expected
+        num_nz_values = 5
+        mdr = haz.get_mdr(cent_idx, impf)
+        self.assertEqual(mdr.nnz, num_nz_values)
+
     def test_get_paa(self):
         haz = dummy_hazard()
         impf = dummy_step_impf(haz)