Merge remote-tracking branch 'origin/develop' into feature/option-appraisal-dataclasses

Author: spjuhel

.gitignore

@@ -51,22 +51,11 @@ coverage.xml
 .hypothesis/
 
 # Temporary data generated by unit tests
-climada/engine/test/data/supplychain/*
-climada/engine/test/data/test.csv
-climada/engine/test/data/test_imp_mat.npz
-climada/engine/test/data/test.xlsx
-climada/entity/impact_funcs/test/data/test_write.xlsx
-climada/hazard/test/data/tc_tracks_nc/*
-climada/hazard/test/data/tc_tracks.h5
-climada/hazard/test/data/test_centr.h5
-climada/hazard/test/data/test_haz.h5
-climada/hazard/test/data/test_unsupported.h5
-climada/test/data/1988234N13299.nc
-climada/test/data/exposure_buildings_default_47_7.h5
-climada/test/data/exposure_high_47_8.h5
-climada/test/data/test_write_hazard.tif
-climada/util/test/data/save_test.pkl
-climada/util/test/data/test_write_raster.tif
+climada/engine/test/data/
+climada/entity/impact_funcs/test/data/
+climada/hazard/test/data/
+climada/test/data/
+climada/util/test/data/
 
 # Translations
 *.mo
@@ -128,51 +117,6 @@ venv.bak/
 # mac finder files
 .DS_Store
 
-# climada system data files
-# they get downloaded separately from the repo
-data/system/global_coast*
-data/system/global_country*
-data/system/tmp_elevation.tif
-
-# climada system data files: Nightlight and population data
-data/system/BlackMarble*tif
-data/system/*stable_lights*.*
-data/system/gpw_v4_population_count_rev11_2015_30_sec.tif
-
-# climada system data files: IBTrACS data
-data/system/IBTrACS.ALL.v04r00.nc
-data/system/data_master_sl_short/
-
-# climada system data files: SPEI data
-data/system/spei06.nc
-
-# climada system data files: World Bank data
-data/system/OGHIST.xls
-data/system/Wealth-Accounts_CSV/
-
-# climada system data files: Credit Suisse Research institute data:
-data/system/GDP2Asset_factors_CRI_2016
-
-# climada system data files: SPAM data
-data/system/cell5m_allockey_xy.csv
-data/system/LicenseV*.pdf
-data/system/readme_global_v3r2.txt
-data/system/spam*csv
-
-# climada system data files: GPW data
-data/system/gpw-v*
-data/system/gpw_v*.tif
-
-# climada system data files: NASA distance to coast
-data/system/GMT_intermediate_coast*
-
-# climada system data files: folders for hazard and exposure data
-data/system/hazard/
-data/system/litpop/
-data/system/litpop_2014/
-
-# climada data: ISIMIP crop data folder:
-data/ISIMIP_crop/
 
 # climada data results folder:
 results/
@@ -183,5 +127,3 @@ results/
 !.pylintrc
 !.readthedocs.yml
 !.github
-climada_petals/
-climada_python/

CHANGELOG.md

@@ -14,6 +14,7 @@ Code freeze date: YYYY-MM-DD
 
 - Better type hints and overloads signatures for ImpactFuncSet [#1250](https://github.com/CLIMADA-project/climada_python/pull/1250)
 - Adds `MeasureConfig` and related dataclasses for new `Measure` object retrocompatibility and (de)serialization capabilities [#1276](https://github.com/CLIMADA-project/climada_python/pull/1276)
+- Add inter- and extrapolation options to `ImpactFreqCurve` with method `interpolate` [#1252](https://github.com/CLIMADA-project/climada_python/pull/1252)
 
 ### Changed
 - Updated Impact Calculation Tutorial (`doc.climada_engine_Impact.ipynb`) [#1095](https://github.com/CLIMADA-project/climada_python/pull/1095).
@@ -22,8 +23,10 @@ Code freeze date: YYYY-MM-DD
 ### Fixed
 
 - Fixed asset count in impact logging message [#1195](https://github.com/CLIMADA-project/climada_python/pull/1195).
+- `Hazard.from_raster_xarray` now returns a sparse matrix instead of a sparse array [#1261](https://github.com/CLIMADA-project/climada_python/pull/1261).
 
 ### Deprecated
+- `Impact.calc_freq_curve()` should not be given the parameter `return_per`. Use the parameter `return_periods` in `Impact.calc_freq_curve().interpolate()` instead.
 
 ### Removed
 - `climada.util.earth_engine.py` Google Earth Engine methods did not facilitate direct use of GEE data in CLIMADA. Code was relocated to [climada-snippets](https://github.com/CLIMADA-project/climada-snippets). [#1109](https://github.com/CLIMADA-project/climada_python/pull/1109)

climada/engine/impact.py

@@ -576,15 +576,6 @@ def local_exceedance_impact(
             self.frequency_unit
         )
 
-        # check method
-        if method not in [
-            "interpolate",
-            "extrapolate",
-            "extrapolate_constant",
-            "stepfunction",
-        ]:
-            raise ValueError(f"Unknown method: {method}")
-
         # calculate local exceedance impact
         test_frequency = 1 / np.array(return_periods)
 
@@ -732,15 +723,6 @@ def local_return_period(
             self.frequency_unit
         )
 
-        # check method
-        if method not in [
-            "interpolate",
-            "extrapolate",
-            "extrapolate_constant",
-            "stepfunction",
-        ]:
-            raise ValueError(f"Unknown method: {method}")
-
         return_periods = np.full((self.imp_mat.shape[1], len(threshold_impact)), np.nan)
 
         nonzero_centroids = np.where(self.imp_mat.getnnz(axis=0) > 0)[0]
@@ -804,6 +786,12 @@ def calc_freq_curve(self, return_per=None):
         ifc_impact = self.at_event[sort_idxs][::-1]
 
         if return_per is not None:
+            warnings.warn(
+                "Calculating the frequency curve on user-specified return periods is deprecated. "
+                "Use ImpactFreqCurve.calc_freq_curve().interpolate() instead.",
+                DeprecationWarning,
+                stacklevel=2,
+            )
             interp_imp = np.interp(return_per, ifc_return_per, ifc_impact)
             ifc_return_per = return_per
             ifc_impact = interp_imp
@@ -2299,15 +2287,108 @@ def plot(self, axis=None, log_frequency=False, **kwargs):
         -------
         matplotlib.axes.Axes
         """
+        # check frequency unit
+        return_period_unit = u_dt.convert_frequency_unit_to_time_unit(
+            self.frequency_unit
+        )
+
         if not axis:
             _, axis = plt.subplots(1, 1)
         axis.set_title(self.label)
         axis.set_ylabel("Impact (" + self.unit + ")")
+
         if log_frequency:
             axis.set_xlabel(f"Exceedance frequency ({self.frequency_unit})")
             axis.set_xscale("log")
             axis.plot(self.return_per**-1, self.impact, **kwargs)
+
         else:
-            axis.set_xlabel("Return period (year)")
+            axis.set_xlabel(f"Return period ({return_period_unit})")
             axis.plot(self.return_per, self.impact, **kwargs)
+
         return axis
+
+    def interpolate(
+        self,
+        return_periods,
+        *,
+        method="interpolate",
+        log_frequency=True,
+        log_impact=True,
+        min_impact=0,
+        bin_decimals=None,
+        y_asymptotic=0.0,
+    ):
+        """Interpolate and extrapolate impact frequency curve using different methods.
+
+        Parameters
+        ----------
+        return_periods : Iterable[float]
+            return periods for which to evaluate the impact frequency curve
+
+        method : str, optional
+            Method to interpolate to new return periods. Currently available are "interpolate",
+            "extrapolate", "extrapolate_constant" and "stepfunction". If set to "interpolate",
+            return periods outside the range of the Impact object's observed return periods
+            will be assigned NaN. If set to "extrapolate_constant" or "stepfunction",
+            return periods larger than the Impact object's observed return periods will be
+            assigned the largest impact, and return periods smaller than the Impact object's
+            observed return periods will be assigned 0. If set to "extrapolate",
+            exceedance impacts will be extrapolated (and interpolated). The extrapolation to
+            large return periods uses the two highest impacts of the centroid and their return
+            periods and extends the interpolation between these points to the given return period
+            (similar for small return periods). Defauls to "interpolate".
+        min_impact : float, optional
+            Minimum threshold to filter the impact. Defaults to 0.
+        log_frequency : bool, optional
+            If set to True, (cummulative) frequency values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
+        log_impact : bool, optional
+            If set to True, impact values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
+        bin_decimals : int, optional
+            Number of decimals to group and bin impact values. Binning results in smoother (and
+            coarser) interpolation and more stable extrapolation. For more details and sensible
+            values for bin_decimals, see Notes. If None, values are not binned. Defaults to None.
+        y_asymptotic : float, optional
+            Has no effect if method is "interpolate". Else, if data size < 2 or if method
+            is set to "extrapolate_constant" or "stepfunction", it provides return value for
+            exceeded impact for return periods smaller than the data range. Defaults to 0.
+
+        Returns
+        -------
+        ImpactFreqCurve
+            impact frequency curve with inter- and extrapolated values
+
+        See Also
+        --------
+        util.interpolation.preprocess_and_interpolate_ev :
+            inter- and extrapolation method
+        """
+        exceedance_frequency = 1 / np.array(return_periods)
+
+        # sort return periods of ImpactFreqCurve
+        sorted_idxs = np.argsort(self.return_per)
+        impacts = np.squeeze(np.array(self.impact)[sorted_idxs])
+        rps = np.asarray(self.return_per)[sorted_idxs]
+        frequency = np.diff(1 / np.array(rps)[::-1], prepend=0)[::-1]
+        impact_interpolated = u_interp.preprocess_and_interpolate_ev(
+            exceedance_frequency,
+            None,
+            frequency,
+            impacts,
+            log_frequency=log_frequency,
+            log_values=log_impact,
+            value_threshold=min_impact,
+            method=method,
+            y_asymptotic=y_asymptotic,
+            bin_decimals=bin_decimals,
+        )
+
+        return ImpactFreqCurve(
+            return_per=return_periods,
+            impact=impact_interpolated,
+            unit=self.unit,
+            frequency_unit=self.frequency_unit,
+            label=self.label,
+        )

climada/engine/test/test_impact.py

@@ -321,6 +321,61 @@ def test_ref_value_rp_pass(self):
         self.assertEqual("USD", ifc.unit)
         self.assertEqual("1/week", ifc.frequency_unit)
 
+    def test_interpolate_freq_curve(self):
+        """Test inter- and extrapolate method of freq curve"""
+        imp = Impact()
+        imp.frequency = np.array([0.2, 0.1, 0.1, 0.1])
+        imp.at_event = np.array([0.0, 100.0, 50.0, 110.0])
+        imp.unit = "USD"
+        imp.frequency_unit = "1/year"
+
+        ifc = imp.calc_freq_curve()
+        # the ifc has values
+        # impacts [0, 50, 100, 110]
+        # exceedance frequencies [.5, .3, .2, .1]
+        # return periods [2, 3.3, 5, 10]
+
+        # stepfunction assigns zero return periods below data and max(impact) for those above
+        npt.assert_array_almost_equal(
+            ifc.interpolate([1, 5, 20], method="stepfunction").impact,
+            [0.0, 100.0, 110.0],
+        )
+
+        # interpolate assigns nan to return periods outside of data
+        npt.assert_array_almost_equal(
+            ifc.interpolate([1, 5, 20], method="interpolate").impact,
+            [np.nan, 100.0, np.nan],
+        )
+
+        # extrapolate_constant assigns zero return periods below data and max(impact) for those above
+        npt.assert_array_almost_equal(
+            ifc.interpolate([1, 5, 20], method="extrapolate_constant").impact,
+            [0.0, 100.0, 110.0],
+        )
+
+        # by binning the last two digits, 100 and 110 are rounded to 100
+        npt.assert_array_almost_equal(
+            ifc.interpolate(
+                [1, 5, 20], method="extrapolate_constant", bin_decimals=-2
+            ).impact,
+            [0.0, 100.0, 100.0],
+        )
+
+        # extrapolation is done by neglecting 0 impacts (min_impact=0)
+        # rp=1: extrapolate impacts [50, 100] and ex_freqs [.3, .2] to ex_freq=1 --> 0
+        # rp=2.5: extrapolate impacts [50, 100] and ex_freqs [.3, .2] to ex_freq=0.4 --> 0
+        # rp=4: extrapolate impacts [50, 100] and ex_freqs [.3, .2] to ex_freq=0.25 --> 75
+        # rp=1: extrapolate impacts [100, 110] and ex_freqs [.2, .1] to ex_freq=0.05 --> 115
+        npt.assert_array_almost_equal(
+            ifc.interpolate(
+                [1.0, 2.5, 4, 20],
+                method="extrapolate",
+                log_frequency=False,
+                log_impact=False,
+            ).impact,
+            [-300.0, 0.0, 75.0, 115.0],
+        )
+
 
 class TestImpactPerYear(unittest.TestCase):
     """Test calc_impact_year_set method"""

climada/entity/impact_funcs/test/test_tc.py

@@ -24,7 +24,11 @@
 import numpy as np
 import pandas as pd
 
-from climada.entity.impact_funcs.trop_cyclone import ImpfSetTropCyclone, ImpfTropCyclone
+from climada.entity.impact_funcs.trop_cyclone import (
+    CountryCode,
+    ImpfSetTropCyclone,
+    ImpfTropCyclone,
+)
 
 
 class TestEmanuelFormula(unittest.TestCase):
@@ -168,6 +172,16 @@ def test_get_countries_per_region(self):
         self.assertListEqual(out[2], [124, 840])
         self.assertListEqual(out[3], ["CAN", "USA"])
 
+    def test_get_countries_per_region_all_or_none(self):
+        ifs = ImpfSetTropCyclone()
+        out = ifs.get_countries_per_region()
+        out2 = ifs.get_countries_per_region("all")
+        self.assertEqual(out, out2)
+        for reg in CountryCode.REGION_NAME.value:
+            out_reg = ifs.get_countries_per_region(reg)
+            for i in range(4):
+                self.assertEqual(out[i][reg], out_reg[i])
+
     def test_get_imf_id_regions_per_countries(self):
         """Test get_impf_id_regions_per_countries()"""
         ifs = ImpfSetTropCyclone()