add impact.select reset_frequency option

AUTHORS.md

@@ -31,3 +31,4 @@
 * Leonie Villiger
 * Kam Lam Yeung
 * Sarah Hülsen
+* Timo Schmid

CHANGELOG.md

@@ -14,6 +14,7 @@ Code freeze date: YYYY-MM-DD
 
 - Convenience method `api_client.Client.get_dataset_file`, combining `get_dataset_info` and `download_dataset`, returning a single file objet. [#821](https://github.com/CLIMADA-project/climada_python/pull/821)
 - Read and Write methods to and from csv files for the `DiscRates` class. [#818](ttps://github.com/CLIMADA-project/climada_python/pull/818)
+- Add reset_frequency option for the impact.select() function. [#847](https://github.com/CLIMADA-project/climada_python/pull/847)
 
 ### Changed
 

climada/engine/impact.py

@@ -1475,9 +1475,9 @@
 
         return imp_fit
 
     def select(self,
                event_ids=None, event_names=None, dates=None,
-               coord_exp=None):
+               coord_exp=None,reset_frequency=False):
         """
         Select a subset of events and/or exposure points from the impact.
         If multiple input variables are not None, it returns all the impacts
@@ -1509,6 +1509,9 @@
         coord_exp : np.array, optional
             Selection of exposures coordinates [lat, lon] (in degrees)
             The default is None.
+        reset_frequency : bool, optional
+            Change frequency of events proportional to difference between first and last
+            year (old and new). Default: False.
 
         Raises
         ------
@@ -1580,6 +1583,19 @@
             LOGGER.info("The total value cannot be re-computed for a "
                         "subset of exposures and is set to None.")
 
+        # reset frequency if date span has changed (optional):
+        if reset_frequency:
+            if self.frequency_unit not in ['1/year', 'annual', '1/y', '1/a']:
+                LOGGER.warning("Resetting the frequency is based on the calendar year of given"
+                    " dates but the frequency unit here is %s. Consider setting the frequency"
+                    " manually for the selection or changing the frequency unit to %s.",
+                    self.frequency_unit, DEF_FREQ_UNIT)
+            year_span_old = np.abs(dt.datetime.fromordinal(self.date.max()).year -
+                                   dt.datetime.fromordinal(self.date.min()).year) + 1
+            year_span_new = np.abs(dt.datetime.fromordinal(imp.date.max()).year -
+                                   dt.datetime.fromordinal(imp.date.min()).year) + 1
+            imp.frequency = imp.frequency * year_span_old / year_span_new
+
         # cast frequency vector into 2d array for sparse matrix multiplication
         freq_mat = imp.frequency.reshape(len(imp.frequency), 1)
         # .A1 reduce 1d matrix to 1d array

climada/engine/test/test_impact.py

@@ -27,6 +27,7 @@
 import h5py
 from pyproj import CRS
 from rasterio.crs import CRS as rCRS
+import datetime as dt
 
 from climada.entity.entity_def import Entity
 from climada.hazard.base import Hazard
@@ -67,6 +68,24 @@ def dummy_impact():
         haz_type="TC",
     )
 
+def dummy_impact_yearly():
+    """Return an impact containing events in multiple years"""
+    imp = dummy_impact()
+
+    years = np.arange(2010,2010+len(imp.date))
+
+    # Edit the date and frequency
+    imp.date = np.array([dt.date(year,1,1).toordinal() for year in years])
+    imp.frequency_unit = "1/year"
+    imp.frequency = np.ones(len(years))/len(years)
+
+    # Calculate the correct expected annual impact
+    freq_mat = imp.frequency.reshape(len(imp.frequency), 1)
+    imp.eai_exp = imp.imp_mat.multiply(freq_mat).sum(axis=0).A1
+    imp.aai_agg = imp.eai_exp.sum()
+
+    return imp
+
 
 class TestImpact(unittest.TestCase):
     """"Test initialization and more"""
@@ -868,6 +887,37 @@ def test_select_imp_map_fail(self):
         with self.assertRaises(ValueError):
             imp.select(event_ids=[0], event_names=[1, 'two'], dates=(0, 2))
 
+    def test_select_reset_frequency(self):
+        """Test that reset_frequency option works correctly"""
+
+        imp = dummy_impact_yearly() # 6 events, 1 per year
+
+        # select first 4 events
+        n_yr = 4
+        sel_imp = imp.select(dates=(imp.date[0],imp.date[n_yr-1]), reset_frequency=True)
+
+        # check unchanged attributes
+        self.assertTrue(u_coord.equal_crs(sel_imp.crs, imp.crs))
+        self.assertEqual(sel_imp.unit, imp.unit)
+        self.assertEqual(sel_imp.frequency_unit, imp.frequency_unit)
+        np.testing.assert_array_equal(sel_imp.coord_exp, imp.coord_exp)
+
+        # check sub-selected attributes
+        np.testing.assert_array_equal(sel_imp.event_id, imp.event_id[0:n_yr])
+        self.assertEqual(sel_imp.event_name, imp.event_name[0:n_yr])
+        np.testing.assert_array_equal(sel_imp.date, imp.date[0:n_yr])
+        np.testing.assert_array_equal(sel_imp.at_event, imp.at_event[0:n_yr])
+        np.testing.assert_array_equal(sel_imp.imp_mat.todense(),
+                                      imp.imp_mat[0:n_yr,:].todense())
+
+        # check frequency-related attributes
+        np.testing.assert_array_equal(sel_imp.frequency, [1/n_yr]*n_yr)
+        self.assertEqual(sel_imp.aai_agg,imp.at_event[0:n_yr].sum()/n_yr)
+
+        np.testing.assert_array_equal(sel_imp.eai_exp,
+                                      imp.imp_mat[0:n_yr,:].todense().sum(axis=0).A1/n_yr)
+
+
 class TestConvertExp(unittest.TestCase):
     def test__build_exp(self):
         """Test that an impact set can be converted to an exposure"""