add at_reg_id_event

climada/engine/impact.py

@@ -191,8 +191,6 @@ def __init__(self,
         else:
             self.imp_mat = sparse.csr_matrix(np.empty((0, 0)))
 
-
-
     def calc(self, exposures, impact_funcs, hazard, save_mat=False, assign_centroids=True):
         """This function is deprecated, use ``ImpactCalc.impact`` instead.
         """
@@ -390,6 +388,44 @@ def impact_per_year(self, all_years=True, year_range=None):
             year_set[year] = sum(self.at_event[orig_year == year])
         return year_set
 
+    def impact_at_reg(self, agg_regions=None):
+        """Aggregate impact on given aggregation regions. This method works 
+        only if Impact.imp_mat was stored during the impact calculation.
+
+        Parameters
+        ----------
+        agg_regions : np.array, list, pd.Series
+            It has the same lenght as the number of coordinates points in exposure 
+            and reports what macro-regions these points belong to. If no aggregation
+            regions are passed, the method aggregates impact at the admin_0 level.
+            Default is None.
+
+        Returns
+        -------
+        pd.DataFrame
+        """
+        if self.imp_mat.nnz == 0:
+            LOGGER.warning("No Impact.imp_mat was stored during the impact calculation,"
+                           "thus the impact per region cannot be computed")
+            return None
+
+        if agg_regions is None:
+            LOGGER.warning("Aggregation regions were not specified,"
+                           "impact is aggregated at the admin_0 level")
+            agg_regions = pd.Series(u_coord.get_country_code(self.coord_exp[:,0],
+                                                             self.coord_exp[:,1]))
+        elif not isinstance(agg_regions, pd.Series):
+            agg_regions = pd.Series(agg_regions)
+
+        at_reg_event = np.hstack([
+                self.imp_mat[:, np.where(agg_regions == reg)[0]].sum(1)
+                for reg in agg_regions.unique()
+                ])
+
+        at_reg_event = pd.DataFrame(at_reg_event, columns=agg_regions.unique())
+
+        return at_reg_event
+
     def calc_impact_year_set(self,all_years=True, year_range=None):
         """This function is deprecated, use Impact.impact_per_year instead."""
         LOGGER.warning("The use of Impact.calc_impact_year_set is deprecated."

climada/engine/test/test_impact.py

@@ -467,6 +467,48 @@ def test_local_exceedance_imp_pass(self):
         self.assertAlmostEqual(np.max(impact_rp), 2916964966.388219, places=5)
         self.assertAlmostEqual(np.min(impact_rp), 444457580.131494, places=5)
 
+class TestImpactReg(unittest.TestCase):
+    """Test impact aggregation per exposure region id or admin0"""
+    def test_impact_at_reg(self):
+        """Test calc local impacts per region"""
+
+        # Read default hazard file
+        hazard = Hazard.from_hdf5(HAZ_TEST_TC)
+
+        # Read an exposure
+        ent = Entity.from_excel(ENT_DEMO_TODAY)
+        ent.check()
+
+        # Calculate impact
+        impact = ImpactCalc(ent.exposures, ent.impact_funcs, hazard).impact(save_mat=True)
+
+        # Aggregate impact at the admin 0 level
+        at_reg_event = impact.impact_at_reg()
+
+        self.assertEqual(at_reg_event.sum().sum(), impact.at_event.sum())
+        self.assertEqual(at_reg_event.shape[0], impact.at_event.shape[0])
+
+        # Aggregate impact at user-defined aggregation regions
+        region_ids = np.hstack([np.repeat(i, 10) for i in range(5)])
+        at_reg_event = impact.impact_at_reg(region_ids)
+
+        self.assertAlmostEqual(at_reg_event.sum().sum(), impact.at_event.sum(), places=2)
+        self.assertEqual(at_reg_event.shape[0], impact.at_event.shape[0])
+        self.assertListEqual(at_reg_event.columns.tolist(), np.unique(region_ids).tolist())
+
+        self.assertEqual(at_reg_event[0].sum(), 2071193014030.06)
+        self.assertEqual(at_reg_event[1].sum(), 2163005244090.206)
+        self.assertEqual(at_reg_event[2].sum(), 2219969197097.4907)
+        self.assertEqual(at_reg_event[3].sum(), 2203363516026.8047)
+        self.assertEqual(at_reg_event[4].sum(), 1827112892434.1558)
+
+        # Do not save Impact.imp_mat in the impact calculation and does do not aggregate
+        impact = ImpactCalc(ent.exposures, ent.impact_funcs, hazard).impact(save_mat=False)
+        at_reg_event = impact.impact_at_reg()
+
+        self.assertIsNone(at_reg_event)
+
+
 class TestRiskTrans(unittest.TestCase):
     """Test risk transfer methods"""
     def test_risk_trans_pass(self):
@@ -816,6 +858,7 @@ def test__exp_build_event(self):
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestImpactPerYear))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestIO))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRPmatrix))
+    TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestImpactReg))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRiskTrans))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestSelect))
     TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestConvertExp))