Merge pull request #971 from CLIMADA-project/feature/rp_impact

Add local_return_period method to Impact objects and new tutorial

Author: ValentinGebhart

CHANGELOG.md

@@ -12,6 +12,8 @@ Code freeze date: YYYY-MM-DD
 
 ### Added
 
+- `climada.engine.impact.Impact.local_return_period` method [#971](https://github.com/CLIMADA-project/climada_python/pull/971)
+- `doc.tutorial.climada_util_local_exceedance_values.ipynb` tutorial explaining `Hazard.local_exceedance_intensity`, `Hazard.local_return_period`, `Impact.local_exceedance_impact`, and `Impact.local_return_period` methods [#971](https://github.com/CLIMADA-project/climada_python/pull/971)
 - `Hazard.local_exceedance_intensity`, `Hazard.local_return_period` and `Impact.local_exceedance_impact`, that all use the `climada.util.interpolation` module [#918](https://github.com/CLIMADA-project/climada_python/pull/918)
 - `climada.util.interpolation` module for inter- and extrapolation util functions used in local exceedance intensity and return period functions [#930](https://github.com/CLIMADA-project/climada_python/pull/930)
 - `climada.exposures.exposures.Exposures.geometry` property

climada/engine/impact.py

@@ -54,6 +54,7 @@
 from climada.entity import Exposures
 from climada.util.constants import CMAP_IMPACT, DEF_CRS, DEF_FREQ_UNIT
 from climada.util.select import get_attributes_with_matching_dimension
+from climada.util.value_representation import safe_divide
 
 LOGGER = logging.getLogger(__name__)
 
@@ -499,7 +500,7 @@ def local_exceedance_impact(
     ):
         """Compute local exceedance impact for given return periods. The default method
         is fitting the ordered impacts per centroid to the corresponding cummulated
-        frequency with by linear interpolation on log-log scale.
+        frequency with linear interpolation on log-log scale.
 
         Parameters
         ----------
@@ -610,6 +611,113 @@ def local_exceedance_imp(self, return_periods=(25, 50, 100, 250)):
             .T.astype(float)
         )
 
+    def local_return_period(
+        self,
+        threshold_impact=(1000.0, 10000.0),
+        method="interpolate",
+        min_impact=0,
+        log_frequency=True,
+        log_impact=True,
+    ):
+        """Compute local return periods for given threshold impacts. The default method
+        is fitting the ordered impacts per centroid to the corresponding cummulated
+        frequency with linear interpolation on log-log scale.
+
+        Parameters
+        ----------
+        threshold_impact : array_like
+            User-specified impact values for which the return period should be calculated
+            locally (at each centroid). Defaults to (1000, 10000)
+        method : str
+            Method to interpolate to new threshold impacts. Currently available are
+            "interpolate", "extrapolate", "extrapolate_constant" and "stepfunction". If set to
+            "interpolate", threshold impacts outside the range of the Impact object's local
+            impacts will be assigned NaN. If set to "extrapolate_constant" or
+            "stepfunction", threshold impacts larger than the Impacts object's local
+            impacts will be assigned NaN, and threshold impacts smaller than the Impact
+            object's local impacts will be assigned the smallest observed local return period.
+            If set to "extrapolate", local return periods will be extrapolated (and interpolated).
+            Defaults to "interpolate".
+        min_impacts : float, optional
+            Minimum threshold to filter the impact. Defaults to 0.
+        log_frequency : bool, optional
+            This parameter is only used if method is set to "interpolate". If set to True,
+            (cummulative) frequency values are converted to log scale before inter- and
+            extrapolation. Defaults to True.
+        log_impact : bool, optional
+            This parameter is only used if method is set to "interpolate". If set to True,
+            impact values are converted to log scale before inter- and extrapolation.
+            Defaults to True.
+
+        Returns
+        -------
+        gdf : gpd.GeoDataFrame
+            GeoDataFrame containing return periods for given threshold impacts. Each column
+            corresponds to a threshold_impact value, each row corresponds to a centroid. Values
+            in the gdf correspond to the return period for the given centroid and
+            threshold_impact value
+        label : str
+            GeoDataFrame label, for reporting and plotting
+        column_label : function
+            Column-label-generating function, for reporting and plotting
+        """
+
+        LOGGER.info("Computing return period map for impacts: %s", threshold_impact)
+        if self.imp_mat.size == 0:
+            raise ValueError(
+                "Attribute imp_mat is empty. Recalculate Impact"
+                "instance with parameter save_mat=True"
+            )
+
+        # check frequency unit
+        return_period_unit = u_dt.convert_frequency_unit_to_time_unit(
+            self.frequency_unit
+        )
+
+        # check method
+        if method not in [
+            "interpolate",
+            "extrapolate",
+            "extrapolate_constant",
+            "stepfunction",
+        ]:
+            raise ValueError(f"Unknown method: {method}")
+
+        # calculate local return periods
+        return_periods = np.array(
+            [
+                u_interp.preprocess_and_interpolate_ev(
+                    None,
+                    np.array(threshold_impact),
+                    self.frequency,
+                    self.imp_mat.getcol(i_centroid).toarray().flatten(),
+                    log_frequency=log_frequency,
+                    log_values=log_impact,
+                    value_threshold=min_impact,
+                    method=method,
+                    y_asymptotic=np.nan,
+                )
+                for i_centroid in range(self.imp_mat.shape[1])
+            ]
+        )
+        return_periods = safe_divide(1.0, return_periods)
+
+        # create the output GeoDataFrame
+        gdf = gpd.GeoDataFrame(
+            geometry=gpd.points_from_xy(self.coord_exp[:, 1], self.coord_exp[:, 0]),
+            crs=self.crs,
+        )
+        col_names = [f"{thresh_impact}" for thresh_impact in threshold_impact]
+        gdf[col_names] = return_periods
+
+        # create label and column_label
+        label = f"Return Periods ({return_period_unit})"
+        column_label = lambda column_names: [
+            f"Impact: {col} {self.unit}" for col in column_names
+        ]
+
+        return gdf, label, column_label
+
     def calc_freq_curve(self, return_per=None):
         """Compute impact exceedance frequency curve.
 

climada/engine/test/test_impact.py

@@ -642,6 +642,120 @@ def test_local_exceedance_impact_methods(self):
             rtol=0.8,
         )
 
+    def test_local_return_period(self):
+        """Test local return periods with lin lin interpolation"""
+
+        impact = dummy_impact()
+        impact.coord_exp = np.array([np.arange(4), np.arange(4)]).T
+        impact.imp_mat = sparse.csr_matrix(
+            np.array([[0, 0, 1, 2], [0, 0, 4, 4], [0, 0, 1, 1], [0, 1, 1, 3]])
+        )
+        impact.frequency = np.ones(4)
+        # first centroid has impacts None with cum frequencies None
+        # second centroid has impacts 1 with cum frequencies 1
+        # third centroid has impacts 1, 4 with cum frequencies 4, 1
+        # fourth centroid has impacts 1,2,3,4 with cum frequencies 4,3,2,1
+        # testing at threshold impacts (0.5, 2.5, 5)
+        return_stats, _, _ = impact.local_return_period(
+            threshold_impact=(0.5, 2.5, 5),
+            method="extrapolate",
+            log_frequency=False,
+            log_impact=False,
+        )
+
+        np.testing.assert_allclose(
+            return_stats[return_stats.columns[1:]].values,
+            np.array(
+                [
+                    [np.nan, np.nan, np.nan],
+                    [1.0, np.nan, np.nan],
+                    [1 / 4.5, 1 / 2.5, np.nan],
+                    [1 / 4.5, 1 / 2.5, np.nan],
+                ]
+            ),
+        )
+
+    def test_local_return_period_methods(self):
+        """Test local return periods different methods"""
+        impact = dummy_impact()
+        impact.coord_exp = np.array([np.arange(4), np.arange(4)]).T
+        impact.imp_mat = sparse.csr_matrix(
+            np.array([[0, 0, 0, 1e1], [0, 0, 1e1, 1e2], [0, 1e3, 1e3, 1e3]])
+        )
+        impact.frequency = np.array([1.0, 0.1, 0.01])
+        # first centroid has impacts None with cum frequencies None
+        # second centroid has impacts 1e3 with frequencies .01, cum freq .01
+        # third centroid has impacts 1e1, 1e3 with cum frequencies .1, .01, cum freq .11, .01
+        # fourth centroid has impacts 1e1, 1e2, 1e3 with cum frequencies 1., .1, .01, cum freq 1.11, .11, .01
+        # testing at threshold impacts .1, 300, 1e5
+
+        # test stepfunction
+        return_stats, _, _ = impact.local_return_period(
+            threshold_impact=(0.1, 300, 1e5), method="stepfunction"
+        )
+        np.testing.assert_allclose(
+            return_stats.values[:, 1:].astype(float),
+            np.array(
+                [
+                    [np.nan, np.nan, np.nan],
+                    [100, 100, np.nan],
+                    [1 / 0.11, 100, np.nan],
+                    [1 / 1.11, 100, np.nan],
+                ]
+            ),
+        )
+
+        # test log log extrapolation
+        return_stats, _, _ = impact.local_return_period(
+            threshold_impact=(0.1, 300, 1e5), method="extrapolate"
+        )
+        np.testing.assert_allclose(
+            return_stats.values[:, 1:].astype(float),
+            np.array(
+                [
+                    [np.nan, np.nan, np.nan],
+                    [100, 100, np.nan],
+                    [1.0, 30, 1e3],
+                    [0.01, 30, 1e4],
+                ]
+            ),
+            rtol=0.8,
+        )
+
+        # test log log interpolation and extrapolation with constant
+        return_stats, _, _ = impact.local_return_period(
+            threshold_impact=(0.1, 300, 1e5), method="extrapolate_constant"
+        )
+        np.testing.assert_allclose(
+            return_stats.values[:, 1:].astype(float),
+            np.array(
+                [
+                    [np.nan, np.nan, np.nan],
+                    [100, 100, np.nan],
+                    [1 / 0.11, 30, np.nan],
+                    [1 / 1.11, 30, np.nan],
+                ]
+            ),
+            rtol=0.8,
+        )
+
+        # test log log interpolation and no extrapolation
+        return_stats, _, _ = impact.local_return_period(
+            threshold_impact=(0.1, 300, 1e5)
+        )
+        np.testing.assert_allclose(
+            return_stats.values[:, 1:].astype(float),
+            np.array(
+                [
+                    [np.nan, np.nan, np.nan],
+                    [np.nan, np.nan, np.nan],
+                    [np.nan, 30, np.nan],
+                    [np.nan, 30, np.nan],
+                ]
+            ),
+            rtol=0.8,
+        )
+
 
 class TestImpactReg(unittest.TestCase):
     """Test impact aggregation per aggregation region or admin 0"""

climada/hazard/base.py

@@ -494,7 +494,7 @@ def local_exceedance_intensity(
     ):
         """Compute local exceedance intensity for given return periods. The default method
         is fitting the ordered intensitites per centroid to the corresponding cummulated
-        frequency with by linear interpolation on log-log scale.
+        frequency with linear interpolation on log-log scale.
 
         Parameters
         ----------
@@ -615,7 +615,7 @@ def local_return_period(
     ):
         """Compute local return periods for given hazard intensities. The default method
         is fitting the ordered intensitites per centroid to the corresponding cummulated
-        frequency with by linear interpolation on log-log scale.
+        frequency with linear interpolation on log-log scale.
 
         Parameters
         ----------