Merge branch 'develop' into feature/option-appraisal-split-make-legacy

Author: spjuhel

CHANGELOG.md

@@ -13,6 +13,7 @@ Code freeze date: YYYY-MM-DD
 ### Added
 
 - Better type hints and overloads signatures for ImpactFuncSet [#1250](https://github.com/CLIMADA-project/climada_python/pull/1250)
+- 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).
@@ -21,8 +22,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()

climada/entity/impact_funcs/trop_cyclone.py

@@ -23,6 +23,7 @@
 
 import logging
 from enum import Enum
+from typing import Literal, overload
 
 import numpy as np
 import pandas as pd
@@ -92,7 +93,7 @@ class CountryCode(Enum):
             "NAM", "NER", "NGA", "NLD", "NOR", "POL", "PRK", "PRT", "PSE", "REU", "ROU",
             "RUS", "RWA", "SDN", "SEN", "SGP", "SGS", "SJM", "SLE", "SMR", "SPM", "SRB",
             "SSD", "STP", "SVK", "SVN", "SWE", "SYC", "TCD", "TGO", "TUN", "TUR", "UKR",
-            "UMI", "VAT", "XKX", "ZMB",
+            "UMI", "VAT", "XKO", "ZMB",
         ],
     }
 
@@ -363,10 +364,47 @@ def calibrated_regional_vhalf(
             reg_v_half[regions_short[-1]] = np.round(df_reg["v_half"].values[0], 5)
         return reg_v_half
 
+    @overload
+    @staticmethod
+    def get_countries_per_region(
+        region: Literal["all"] = "all",
+    ) -> tuple[
+        dict[str, str],  # region_name
+        dict[str, int],  # impf_id
+        dict[str, list[int]],  # numeric
+        dict[str, list[str]],  # alpha3
+    ]: ...
+
+    @overload
+    @staticmethod
+    def get_countries_per_region(
+        region: None,
+    ) -> tuple[
+        dict[str, str], dict[str, int], dict[str, list[int]], dict[str, list[str]]
+    ]: ...
+
+    @overload
+    @staticmethod
+    def get_countries_per_region(
+        region: str,
+    ) -> tuple[
+        str, int, list[int], list[str]  # region_name  # impf_id  # numeric  # alpha3
+    ]: ...
+
     @staticmethod
     def get_countries_per_region(region=None):
-        """Returns dictionaries with numerical (numeric) and alphabetical (alpha3) ISO3 codes
-        of all countries associated to a calibration region.
+        """Returns countries within a TC calibration region and associated impact functions.
+
+        This method returns a tuple with numerical (numeric) and alphabetical (alpha3)
+        ISO3 codes of all countries associated to a calibration region.
+
+        If no region or "all" is provided as argument, the method return a tuple of
+        dictionaries with short name of the tropical cyclone calibration regions as
+        keys and the values for each of those.
+
+        Notes
+        -----
+
         Only contains countries that were affected by tropical cyclones
         between 1980 and 2017 according to EM-DAT.
 
@@ -395,9 +433,12 @@ def get_countries_per_region(region=None):
             return (
                 CountryCode.REGION_NAME.value,
                 CountryCode.IMPF_ID.value,
-                coordinates.country_to_iso(
-                    CountryCode.ALPHA3.value, representation="numeric"
-                ),
+                {
+                    reg: coordinates.country_to_iso(
+                        CountryCode.ALPHA3.value[reg], representation="numeric"
+                    )
+                    for reg in CountryCode.REGION_NAME.value
+                },
                 CountryCode.ALPHA3.value,
             )
 

climada/hazard/base.py

@@ -554,15 +554,6 @@ def local_exceedance_intensity(
             self.frequency_unit
         )
 
-        # check method
-        if method not in [
-            "interpolate",
-            "extrapolate",
-            "extrapolate_constant",
-            "stepfunction",
-        ]:
-            raise ValueError(f"Unknown method: {method}")
-
         # calculate local exceedance intensity
         test_frequency = 1 / np.array(return_periods)
 
@@ -707,15 +698,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.intensity.shape[1], len(threshold_intensities)), np.nan
         )