diff --git a/CHANGELOG.md b/CHANGELOG.md index e301836bd..b124f3599 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/climada/engine/impact.py b/climada/engine/impact.py index d8c944c7b..e8042afef 100644 --- a/climada/engine/impact.py +++ b/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, + ) diff --git a/climada/engine/test/test_impact.py b/climada/engine/test/test_impact.py index 38b3def3d..c79a264a2 100644 --- a/climada/engine/test/test_impact.py +++ b/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""" diff --git a/climada/hazard/base.py b/climada/hazard/base.py index abcbae2e8..b06c22e98 100644 --- a/climada/hazard/base.py +++ b/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 ) diff --git a/climada/util/interpolation.py b/climada/util/interpolation.py index 67bcc3ee8..bf64f37cd 100644 --- a/climada/util/interpolation.py +++ b/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,10 +137,10 @@ 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, @@ -139,7 +148,7 @@ def preprocess_and_interpolate_ev( ) extrapolation = None if method == "interpolate" else method return _interpolate_ev( - test_frequency, + exceedance_frequency, frequency[::-1], values[::-1], logx=log_frequency,