add risk functions for 100 and 250 return periods

Author: gabrielaznar

climada/engine/cost_benefit.py

@@ -19,7 +19,7 @@
 Define CostBenefit class.
 """
 
-__all__ = ['CostBenefit', 'risk_aai_agg']
+__all__ = ['CostBenefit', 'risk_aai_agg', 'risk_rp_100', 'risk_rp_250']
 
 import logging
 import numpy as np
@@ -52,6 +52,30 @@ def risk_aai_agg(impact):
     """
     return impact.aai_agg
 
+def risk_rp_100(impact):
+    """Risk measurement as exceedance impact at 100 years return period.
+
+    Parameters:
+        impact (Impact): an Impact instance
+
+    Returns:
+        float
+    """
+    efc = impact.calc_freq_curve([100])
+    return efc.impact[0]
+
+def risk_rp_250(impact):
+    """Risk measurement as exceedance impact at 250 years return period.
+
+    Parameters:
+        impact (Impact): an Impact instance
+
+    Returns:
+        float
+    """
+    efc = impact.calc_freq_curve([250])
+    return efc.impact[0]
+
 class CostBenefit():
     """Impact definition. Compute from an entity (exposures and impact
     functions) and hazard.
@@ -140,7 +164,7 @@ def calc(self, hazard, entity, haz_future=None, ent_future=None, \
         self.unit = entity.exposures.value_unit
 
         # save measure colors
-        for meas in entity.measures.get_measure():
+        for meas in entity.measures.get_measure(hazard.tag.haz_type):
             self.color_rgb[meas.name] = meas.color_rgb
 
         if not haz_future and not ent_future:
@@ -295,6 +319,8 @@ def plot_waterfall(self, hazard, entity, haz_future, ent_future,
         time_dep = self._time_dependency_array()
         risk_curr = self._npv_unaverted_impact(risk_future, entity.disc_rates,
                                                time_dep)
+        LOGGER.info('Risk function at {:d}: {:.3e}'.format(self.present_year,
+                    risk_future))
 
         # changing future
         time_dep = self._time_dependency_array(imp_time_depen)
@@ -308,6 +334,8 @@ def plot_waterfall(self, hazard, entity, haz_future, ent_future,
         risk_future = fut_risk
         risk_tot = self._npv_unaverted_impact(risk_future, entity.disc_rates,
                                               time_dep, curr_risk)
+        LOGGER.info('Risk function at {:d}: {:.3e}'.format(self.future_year,
+                    risk_future))
 
         axis.bar(1, risk_curr/norm_fact)
         axis.text(1, risk_curr/norm_fact, str(int(round(risk_curr/norm_fact))), \
@@ -328,6 +356,8 @@ def plot_waterfall(self, hazard, entity, haz_future, ent_future,
         plt.xticks(np.arange(4)+1, ['Risk ' + str(self.present_year), \
             'Economic \ndevelopment', 'Climate \nchange', 'Risk ' + str(self.future_year)])
         axis.set_ylabel('Impact (' + self.unit + ' ' + norm_name + ')')
+        axis.set_title('Total accumulated damage from {:d} to {:d}'.format(
+                self.present_year, self.future_year))
 
         return fig, axis
 
@@ -362,7 +392,7 @@ def _calc_impact_measures(self, hazard, exposures, meas_set, imp_fun_set, \
             impact_meas['no measure']['impact'] = imp_tmp
 
         # compute impact for each measure
-        for measure in meas_set.get_measure():
+        for measure in meas_set.get_measure(hazard.tag.haz_type):
             imp_tmp, risk_transf = measure.calc_impact(exposures, imp_fun_set, hazard)
             impact_meas[measure.name] = dict()
             impact_meas[measure.name]['cost'] = measure.cost
@@ -443,6 +473,9 @@ def _time_dependency_array(self, imp_time_depen=None):
         Parameters:
             imp_time_depen (float, optional): parameter which represent time
                 evolution of impact. Time array is all ones if not provided
+
+        Returns:
+            np.array
         """
         n_years = self.future_year - self.present_year + 1
         if imp_time_depen:

climada/engine/test/test_cost_benefit.py

@@ -26,7 +26,9 @@
 from climada.entity.entity_def import Entity
 from climada.entity.disc_rates import DiscRates
 from climada.hazard.base import Hazard
-from climada.engine.cost_benefit import CostBenefit, risk_aai_agg, DEF_RP
+from climada.engine.cost_benefit import CostBenefit, risk_aai_agg, DEF_RP, \
+risk_rp_100, risk_rp_250
+from climada.engine import Impact
 from climada.util.constants import ENT_DEMO_FUTURE, ENT_DEMO_TODAY
 
 HAZ_DATA_DIR = os.path.join(os.path.dirname(__file__), '../../hazard/test/data')
@@ -45,10 +47,11 @@ def test_calc_impact_measures_pass(self):
         hazard.read_mat(HAZ_TEST_MAT)
         entity = Entity()
         entity.read_mat(ENT_TEST_MAT)
-        entity.exposures.rename(columns={'if_': 'if_TC'}, inplace=True)
         entity.check()
-        for meas in entity.measures.get_measure():
+        entity.measures._data['TC'] = entity.measures._data.pop('XX')
+        for meas in entity.measures.get_measure('TC'):
             meas.haz_type = 'TC'
+        entity.check()
 
         cost_ben = CostBenefit()
         cost_ben._calc_impact_measures(hazard, entity.exposures, entity.measures,
@@ -112,10 +115,10 @@ def test_calc_cb_no_change_pass(self):
         hazard.read_mat(HAZ_TEST_MAT)
         entity = Entity()
         entity.read_mat(ENT_TEST_MAT)
-        entity.exposures.rename(columns={'if_': 'if_TC'}, inplace=True)
-        entity.check()
-        for meas in entity.measures.get_measure():
+        entity.measures._data['TC'] = entity.measures._data.pop('XX')
+        for meas in entity.measures.get_measure('TC'):
             meas.haz_type = 'TC'
+        entity.check()
 
         cost_ben = CostBenefit()
         cost_ben._calc_impact_measures(hazard, entity.exposures, entity.measures,
@@ -148,10 +151,10 @@ def test_calc_cb_change_pass(self):
         hazard.read_mat(HAZ_TEST_MAT)
         entity = Entity()
         entity.read_mat(ENT_TEST_MAT)
-        entity.exposures.rename(columns={'if_': 'if_TC'}, inplace=True)
-        entity.check()
-        for meas in entity.measures.get_measure():
+        entity.measures._data['TC'] = entity.measures._data.pop('XX')
+        for meas in entity.measures.get_measure('TC'):
             meas.haz_type = 'TC'
+        entity.check()
 
         cost_ben = CostBenefit()
         cost_ben._calc_impact_measures(hazard, entity.exposures, entity.measures,
@@ -160,8 +163,6 @@ def test_calc_cb_change_pass(self):
         ent_future = Entity()
         ent_future.read_excel(ENT_DEMO_FUTURE)
         ent_future.check()
-        for meas in ent_future.measures.get_measure():
-            meas.haz_type = 'TC'
 
         haz_future = copy.deepcopy(hazard)
         haz_future.intensity.data += 25
@@ -369,7 +370,46 @@ def test_calc_no_change_pass(self):
 
         self.assertEqual(cost_ben.tot_climate_risk, 1.2150496306913972e+11)
 
+class TestRiskFuncs(unittest.TestCase):
+    '''Test risk functions definitions'''
+
+    def test_impact():
+        ent = Entity()
+        ent.read_excel(ENT_DEMO_TODAY)
+        ent.check()
+        hazard = Hazard('TC')
+        hazard.read_mat(HAZ_TEST_MAT)
+        impact = Impact()
+        ent.exposures.assign_centroids(hazard)
+        impact.calc(ent.exposures, ent.impact_funcs, hazard)
+        return impact        
+
+    def test_risk_aai_agg_pass(self):
+        """Test risk_aai_agg"""
+        impact = self.test_impact()
+        risk = risk_aai_agg(impact)
+        self.assertAlmostEqual(6.512201157564421e+09, risk, 5)
+        self.assertTrue(np.isclose(6.512201157564421e+09, risk))
+
+    def test_risk_rp_100_pass(self):
+        """Test risk_rp_100"""
+        impact = self.test_impact()
+        exc_freq = impact.calc_freq_curve([100])
+
+        risk = risk_rp_100(impact)
+        self.assertAlmostEqual(exc_freq.impact[0], risk)
+
+    def test_risk_rp_200_pass(self):
+        """Test risk_rp_200"""
+        impact = self.test_impact()
+        exc_freq = impact.calc_freq_curve([250])
+
+        risk = risk_rp_250(impact)
+        self.assertAlmostEqual(exc_freq.impact[0], risk)
+
 # Execute Tests
+TESTS = unittest.TestLoader().loadTestsFromTestCase(TestRiskFuncs)
 TESTS = unittest.TestLoader().loadTestsFromTestCase(TestCalc)
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestSteps))
+TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRiskFuncs))
 unittest.TextTestRunner(verbosity=2).run(TESTS)