Add inter/extrapolate features to ImpactFreqCurve

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).
@@ -22,6 +23,7 @@ Code freeze date: YYYY-MM-DD
 - Fixed asset count in impact logging message [#1195](https://github.com/CLIMADA-project/climada_python/pull/1195).
 
 ### 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/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
         )

climada/util/interpolation.py

@@ -29,7 +29,7 @@
 
 
 def preprocess_and_interpolate_ev(
-    test_frequency,
+    exceedance_frequency,
     test_values,
     frequency,
     values,
@@ -46,12 +46,12 @@ def preprocess_and_interpolate_ev(
 
     Parameters
     ----------
-    test_frequency : array_like
-        1-D array of test frequencies for which values (e.g., intensities or impacts) should be
+    exceedance_frequency : array_like
+        1-D array of test exceedance frequencies for which values (e.g., intensities or impacts) should be
         assigned. If given, test_values must be None.
     test_values : array_like
-        1-D array of test values (e.g., intensities or impacts) for which frequencies should be
-        assigned. If given, test_frequency must be None.
+        1-D array of test values (e.g., intensities or impacts) for which exceedance frequencies should be
+        assigned. If given, exceedance_frequency must be None.
     frequency : array_like
         1-D array of frequencies to be interpolated.
     values : array_like
@@ -106,13 +106,22 @@ def preprocess_and_interpolate_ev(
     could use bin_decimals=5.
     """
 
+    # check method
+    if method not in [
+        "interpolate",
+        "extrapolate",
+        "extrapolate_constant",
+        "stepfunction",
+    ]:
+        raise ValueError(f"Unknown method: {method}")
+
     # check that only test frequencies or only test values are given
-    if test_frequency is not None and test_values is not None:
+    if exceedance_frequency is not None and test_values is not None:
         raise ValueError(
             "Both test frequencies and test values are given. This method only handles one of "
             "the two. To use this method, please only use one of them."
         )
-    if test_frequency is None and test_values is None:
+    if exceedance_frequency is None and test_values is None:
         raise ValueError("No test values or test frequencies are given.")
 
     # sort values and frequencies
@@ -128,18 +137,18 @@ def preprocess_and_interpolate_ev(
     frequency = np.cumsum(frequency[::-1])[::-1]
 
     # if test frequencies are provided
-    if test_frequency is not None:
+    if exceedance_frequency is not None:
         if method == "stepfunction":
             return _stepfunction_ev(
-                test_frequency,
+                exceedance_frequency,
                 frequency[::-1],
                 values[::-1],
                 y_threshold=value_threshold,
                 y_asymptotic=y_asymptotic,
             )
         extrapolation = None if method == "interpolate" else method
         return _interpolate_ev(
-            test_frequency,
+            exceedance_frequency,
             frequency[::-1],
             values[::-1],
             logx=log_frequency,