Add inter/extrapolate features to ImpactFreqCurve

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]
@@ -2299,15 +2281,208 @@ def plot(self, axis=None, log_frequency=False, **kwargs):
         -------
         matplotlib.axes.Axes
         """
+
+        return self._plot(
+            self.return_per,
+            self.impact,
+            axis,
+            log_frequency,
+            self.frequency_unit,
+            self.unit,
+            self.label,
+            **kwargs,
+        )
+
+    def evaluate(
+        self,
+        return_period,
+        *,
+        method="interpolate",
+        log_frequency=True,
+        log_impact=True,
+        min_impact=0,
+        bin_decimals=None,
+    ):
+        """Evaluate impact frequency curve with different extrapolation options.
+
+        Parameters
+        ----------
+        return_period : Iterable[float]
+            return periods for which to evaluate the impact frequency curve
+            axis to use
+        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.
+
+        Returns
+        -------
+        np.array
+            array of exceeded impacts at given return periods
+        """
+        exceedance_frequency = 1 / np.array(return_period)
+
+        # 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)
+
+        return 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=0.0,
+            bin_decimals=bin_decimals,
+        )
+
+    def plot_extended(
+        self,
+        *,
+        return_period_range=None,
+        axis=None,
+        log_frequency=False,
+        kwargs_interp=None,
+        **kwargs,
+    ):
+        """Plot impact frequency curve with extrapolation options.
+
+        Parameters
+        ----------
+        return_period_range : tuple(float), optional
+            range of return period values to plot, e.g. (5, 500) to plot between 5 and 500 years.
+            If None, the plot range is between 0.5*min(data_return_periods) and
+            1.2*max(data_return_periods), where data_return_periods are the return period values
+            extracted from the impact object's data. Defaults to None.
+        axis : matplotlib.axes.Axes, optional
+            axis to use
+        log_frequency : boolean, optional
+            plot logarithmioc exceedance frequency on x-axis
+        kwargs_interp : dict, optional
+            dict with (key, value) pairs to handle inter- and extrapolation behaviour, e.g.
+            {"method": "extrapolate"}. Default is {"method": "interpolate", "log_frequency": True,
+            "log_impact": True, "min_impact": 0, "bin_decimals": None}. Available options are
+            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.
+        kwargs : dict, optional
+            arguments for plot matplotlib function, e.g. color='b'
+
+        Returns
+        -------
+        matplotlib.axes.Axes
+        """
+        if kwargs_interp is None:
+            kwargs_interp = {}
+        kwargs_interp = {
+            "method": "interpolate",
+            "log_frequency": True,
+            "log_impact": True,
+            "min_impact": 0,
+            "bin_decimals": None,
+        } | kwargs_interp
+
+        if return_period_range is None:
+            return_periods = np.linspace(
+                0.5 * min(self.return_per), 1.2 * max(self.return_per), 500
+            )
+        else:
+            return_periods = np.linspace(
+                min(return_period_range), max(return_period_range), 500
+            )
+
+        impacts = self.evaluate(return_periods, **kwargs_interp)
+
+        axis = self._plot(
+            return_periods,
+            impacts,
+            axis,
+            log_frequency,
+            self.frequency_unit,
+            self.unit,
+            self.label,
+            **kwargs,
+        )
+        # axis = self._plot(self.return_per, self.evaluate(self.return_per, **(kwargs_interp|{"method": "interpolate"})), axis, log_frequency, self.frequency_unit, self.unit, self.label, **kwargs)
+        for rp in [min(self.return_per), max(self.return_per)]:
+            axis.axvline(x=1 / rp if log_frequency else rp, linestyle="--", c="gray")
+
+        return axis
+
+    @classmethod
+    def _plot(
+        cls,
+        return_per,
+        impact,
+        axis,
+        log_frequency,
+        frequency_unit,
+        unit,
+        title,
+        **kwargs,
+    ):
+        """
+        private function to plot an impact's exceedance frequency curve
+        """
+        # check frequency unit
+        return_period_unit = u_dt.convert_frequency_unit_to_time_unit(frequency_unit)
+
         if not axis:
             _, axis = plt.subplots(1, 1)
-        axis.set_title(self.label)
-        axis.set_ylabel("Impact (" + self.unit + ")")
+        axis.set_title(title)
+        axis.set_ylabel("Impact (" + unit + ")")
         if log_frequency:
-            axis.set_xlabel(f"Exceedance frequency ({self.frequency_unit})")
+            axis.set_xlabel(f"Exceedance frequency ({frequency_unit})")
             axis.set_xscale("log")
-            axis.plot(self.return_per**-1, self.impact, **kwargs)
+            axis.plot(return_per**-1, impact, **kwargs)
         else:
-            axis.set_xlabel("Return period (year)")
-            axis.plot(self.return_per, self.impact, **kwargs)
+            axis.set_xlabel(f"Return period ({return_period_unit})")
+            axis.plot(return_per, impact, **kwargs)
         return axis

climada/engine/test/test_impact.py

@@ -321,6 +321,42 @@ def test_ref_value_rp_pass(self):
         self.assertEqual("USD", ifc.unit)
         self.assertEqual("1/week", ifc.frequency_unit)
 
+    def test_evaluate_freq_curve(self):
+        """Test evaluate method of freq curve"""
+        imp = Impact()
+        imp.frequency = np.ones(4) * 0.1
+        imp.at_event = np.zeros(4)
+        imp.at_event[0] = 0.0
+        imp.at_event[1] = 100.0
+        imp.at_event[2] = 50.0
+        imp.at_event[3] = 110.0
+        imp.unit = "USD"
+        imp.frequency_unit = "1/year"
+
+        ifc = imp.calc_freq_curve()
+        npt.assert_array_almost_equal(
+            ifc.evaluate([1, 5, 20], method="stepfunction"), [0.0, 100.0, 110.0]
+        )
+        npt.assert_array_almost_equal(
+            ifc.evaluate([1, 5, 20], method="interpolate"), [np.nan, 100.0, np.nan]
+        )
+        npt.assert_array_almost_equal(
+            ifc.evaluate([1, 5, 20], method="extrapolate_constant"), [0.0, 100.0, 110.0]
+        )
+        npt.assert_array_almost_equal(
+            ifc.evaluate([1, 5, 20], method="extrapolate_constant", bin_decimals=-2),
+            [0.0, 100.0, 100.0],
+        )
+        npt.assert_array_almost_equal(
+            ifc.evaluate(
+                [1.0, 2.5, 4, 20],
+                method="extrapolate",
+                log_frequency=False,
+                log_impact=False,
+            ),
+            [-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,11 +46,11 @@ 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
+        1-D array of test values (e.g., intensities or impacts) for which exceedance frequencies should be
         assigned. If given, test_frequency must be None.
     frequency : array_like
         1-D array of frequencies to be interpolated.
@@ -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,