Improve efficiency of local exceedance intensitry/impact and local return periods functions

CHANGELOG.md

@@ -13,12 +13,13 @@ Code freeze date: YYYY-MM-DD
 ### Added
 
 ### Changed
-
+- `Hazard.local_exceedance_intensity`, `Hazard.local_return_period` and `Impact.local_exceedance_impact`, `Impact.local_return_period`, using the `climada.util.interpolation` module: New default (no binning), binning on decimals, and faster implementation [#1012](https://github.com/CLIMADA-project/climada_python/pull/1012)
 ### Fixed
 
 ### Deprecated
 
 ### Removed
+- `climada.util.interpolation.round_to_sig_digits` [#1012](https://github.com/CLIMADA-project/climada_python/pull/1012)
 
 ## 6.0.1
 

climada/engine/test/test_impact.py

@@ -556,6 +556,7 @@ def test_local_exceedance_impact(self):
             method="extrapolate",
             log_frequency=False,
             log_impact=False,
+            bin_decimals=2,
         )
         np.testing.assert_allclose(
             impact_stats.values[:, 1:].astype(float),
@@ -594,7 +595,7 @@ def test_local_exceedance_impact_methods(self):
 
         # test log log extrapolation
         impact_stats, _, _ = impact.local_exceedance_impact(
-            return_periods=(1000, 30, 0.1), method="extrapolate"
+            return_periods=(1000, 30, 0.1), method="extrapolate", bin_decimals=-1
         )
         np.testing.assert_allclose(
             impact_stats.values[:, 1:].astype(float),
@@ -661,6 +662,7 @@ def test_local_return_period(self):
             method="extrapolate",
             log_frequency=False,
             log_impact=False,
+            bin_decimals=2,
         )
 
         np.testing.assert_allclose(
@@ -707,7 +709,7 @@ def test_local_return_period_methods(self):
 
         # test log log extrapolation
         return_stats, _, _ = impact.local_return_period(
-            threshold_impact=(0.1, 300, 1e5), method="extrapolate"
+            threshold_impact=(0.1, 300, 1e5), method="extrapolate", bin_decimals=2
         )
         np.testing.assert_allclose(
             return_stats.values[:, 1:].astype(float),

climada/hazard/base.py

@@ -488,6 +488,7 @@
         min_intensity=None,
         log_frequency=True,
         log_intensity=True,
+        bin_decimals=None,
     ):
         """Compute local exceedance intensity for given return periods. The default method
         is fitting the ordered intensitites per centroid to the corresponding cummulated
@@ -506,19 +507,23 @@
             periods larger than the Hazard object's observed local return periods will be assigned
             the largest local intensity, and return periods smaller than the Hazard object's
             observed local return periods will be assigned 0. If set to "extrapolate", local
-            exceedance intensities will be extrapolated (and interpolated).
-            Defauls to "interpolate".
+            exceedance intensities will be extrapolated (and interpolated). The extrapolation to
+            large return periods uses the two highest intensites 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_intensity : float, optional
             Minimum threshold to filter the hazard intensity. If set to None, self.intensity_thres
             will be used. Defaults to None.
         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.
+            If set to True, (cummulative) frequency values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
         log_intensity : bool, optional
-            This parameter is only used if method is set to "interpolate". If set to True,
-            intensity values are converted to log scale before inter- and extrapolation.
-            Defaults to True.
+            If set to True, intensity values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
+        bin_decimals : int, optional
+            Number of decimals to group and bin intensity 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
         -------
@@ -531,6 +536,23 @@
             GeoDataFrame label, for reporting and plotting
         column_label : function
             Column-label-generating function, for reporting and plotting
+
+        See Also
+        --------
+        util.interpolation.preprocess_and_interpolate_ev :
+            inter- and extrapolation method
+
+        Notes
+        -------
+        If an integer bin_decimals is given, the intensity values are binned according to their
+        bin_decimals decimals, and their corresponding frequencies are summed. This binning leads
+        to a smoother (and coarser) interpolation, and a more stable extrapolation. For instance,
+        if bin_decimals=1, the two values 12.01 and 11.97 with corresponding frequencies 0.1 and
+        0.2 are combined to a value 12.0 with frequency 0.3. The default bin_decimals=None results
+        in not binning the values.
+        E.g., if your intensities range from 1 to 100, you could use bin_decimals=1, if your
+        intensities range from 1e6 to 1e9, you could use bin_decimals=-5, if your intensities
+        range from 0.0001 to .01, you could use bin_decimals=5.
         """
         if not min_intensity and min_intensity != 0:
             min_intensity = self.intensity_thres
@@ -550,23 +572,32 @@
 
         # calculate local exceedance intensity
         test_frequency = 1 / np.array(return_periods)
-        exceedance_intensity = np.array(
-            [
-                u_interp.preprocess_and_interpolate_ev(
-                    test_frequency,
-                    None,
-                    self.frequency,
-                    self.intensity.getcol(i_centroid).toarray().flatten(),
-                    log_frequency=log_frequency,
-                    log_values=log_intensity,
-                    value_threshold=min_intensity,
-                    method=method,
-                    y_asymptotic=0.0,
-                )
-                for i_centroid in range(self.intensity.shape[1])
-            ]
+
+        exceedance_intensity = np.full(
+            (self.intensity.shape[1], len(test_frequency)),
+            np.nan if method == "interpolate" else 0.0,
         )
 
+        nonzero_centroids = np.where(self.intensity.getnnz(axis=0) > 0)[0]
+        if not len(nonzero_centroids) == 0:
+            exceedance_intensity[nonzero_centroids, :] = np.array(
+                [
+                    u_interp.preprocess_and_interpolate_ev(
+                        test_frequency,
+                        None,
+                        self.frequency,
+                        self.intensity.getcol(i_centroid).toarray().flatten(),
+                        log_frequency=log_frequency,
+                        log_values=log_intensity,
+                        value_threshold=min_intensity,
+                        method=method,
+                        y_asymptotic=0.0,
+                        bin_decimals=bin_decimals,
+                    )
+                    for i_centroid in nonzero_centroids
+                ]
+            )
+
         # create the output GeoDataFrame
         gdf = gpd.GeoDataFrame(
             geometry=self.centroids.gdf["geometry"], crs=self.centroids.gdf.crs
@@ -576,9 +607,11 @@
 
         # create label and column_label
         label = f"Intensity ({self.units})"
-        column_label = lambda column_names: [
-            f"Return Period: {col} {return_period_unit}" for col in column_names
-        ]
+
+        def column_label(column_names):
+            return [
+                f"Return Period: {col} {return_period_unit}" for col in column_names
+            ]
 
         return gdf, label, column_label
 
@@ -609,6 +642,7 @@
         min_intensity=None,
         log_frequency=True,
         log_intensity=True,
+        bin_decimals=None,
     ):
         """Compute local return periods for given hazard intensities. The default method
         is fitting the ordered intensitites per centroid to the corresponding cummulated
@@ -628,18 +662,24 @@
             intensities will be assigned NaN, and threshold intensities smaller than the Hazard
             object's local intensities will be assigned the smallest observed local return period.
             If set to "extrapolate", local return periods will be extrapolated (and interpolated).
+            The extrapolation to large threshold intensities uses the two highest intensites of
+            the centroid and their return periods and extends the interpolation between these
+            points to the given threshold intensity (similar for small threshold intensites).
             Defaults to "interpolate".
         min_intensity : float, optional
             Minimum threshold to filter the hazard intensity. If set to None, self.intensity_thres
             will be used. Defaults to None.
         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.
+            If set to True, (cummulative) frequency values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
         log_intensity : bool, optional
-            This parameter is only used if method is set to "interpolate". If set to True,
-            intensity values are converted to log scale before inter- and extrapolation.
-            Defaults to True.
+            If set to True, intensity values are converted to log scale before
+            inter- and extrapolation. Defaults to True.
+        bin_decimals : int, optional
+            Number of decimals to group and bin intensity 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
         -------
@@ -652,6 +692,23 @@
             GeoDataFrame label, for reporting and plotting
         column_label : function
             Column-label-generating function, for reporting and plotting
+
+        See Also
+        --------
+        util.interpolation.preprocess_and_interpolate_ev :
+            inter- and extrapolation method
+
+        Notes
+        -------
+        If an integer bin_decimals is given, the intensity values are binned according to their
+        bin_decimals decimals, and their corresponding frequencies are summed. This binning leads
+        to a smoother (and coarser) interpolation, and a more stable extrapolation. For instance,
+        if bin_decimals=1, the two values 12.01 and 11.97 with corresponding frequencies 0.1 and
+        0.2 are combined to a value 12.0 with frequency 0.3. The default bin_decimals=None results
+        in not binning the values.
+        E.g., if your intensities range from 1 to 100, you could use bin_decimals=1, if your
+        intensities range from 1e6 to 1e9, you could use bin_decimals=-5, if your intensities
+        range from 0.0001 to .01, you could use bin_decimals=5.
         """
         if not min_intensity and min_intensity != 0:
             min_intensity = self.intensity_thres
@@ -669,23 +726,31 @@
         ]:
             raise ValueError(f"Unknown method: {method}")
 
-        # calculate local return periods
-        return_periods = np.array(
-            [
-                u_interp.preprocess_and_interpolate_ev(
-                    None,
-                    np.array(threshold_intensities),
-                    self.frequency,
-                    self.intensity.getcol(i_centroid).toarray().flatten(),
-                    log_frequency=log_frequency,
-                    log_values=log_intensity,
-                    value_threshold=min_intensity,
-                    method=method,
-                    y_asymptotic=np.nan,
-                )
-                for i_centroid in range(self.intensity.shape[1])
-            ]
+        return_periods = np.full(
+            (self.intensity.shape[1], len(threshold_intensities)), np.nan
         )
+
+        nonzero_centroids = np.where(self.intensity.getnnz(axis=0) > 0)[0]
+
+        if not len(nonzero_centroids) == 0:
+            return_periods[nonzero_centroids, :] = np.array(
+                [
+                    u_interp.preprocess_and_interpolate_ev(
+                        None,
+                        np.array(threshold_intensities),
+                        self.frequency,
+                        self.intensity.getcol(i_centroid).toarray().flatten(),
+                        log_frequency=log_frequency,
+                        log_values=log_intensity,
+                        value_threshold=min_intensity,
+                        method=method,
+                        y_asymptotic=np.nan,
+                        bin_decimals=bin_decimals,
+                    )
+                    for i_centroid in nonzero_centroids
+                ]
+            )
+
         return_periods = safe_divide(1.0, return_periods)
 
         # create the output GeoDataFrame
@@ -697,9 +762,9 @@
 
         # create label and column_label
         label = f"Return Periods ({return_period_unit})"
-        column_label = lambda column_names: [
-            f"Threshold Intensity: {col} {self.units}" for col in column_names
-        ]
+
+        def column_label(column_names):
+            return [f"Threshold Intensity: {col} {self.units}" for col in column_names]
 
         return gdf, label, column_label
 

climada/hazard/plot.py

@@ -23,7 +23,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
 from deprecation import deprecated
 
 import climada.util.plot as u_plot
 
@@ -37,62 +37,62 @@
     Contains all plotting methods of the Hazard class
     """
 
-    @deprecated(
-        details="The use of Hazard.plot_rp_intensity is deprecated."
-        "Use Hazard.local_exceedance_intensity and util.plot.plot_from_gdf instead."
-    )
     def plot_rp_intensity(
         self,
         return_periods=(25, 50, 100, 250),
-        smooth=True,
         axis=None,
-        figsize=(9, 13),
-        adapt_fontsize=True,
+        kwargs_local_exceedance_intensity=None,
         **kwargs,
     ):
         """
-        This function is deprecated,
-        use Impact.local_exceedance_impact and util.plot.plot_from_gdf instead.
-
         Compute and plot hazard exceedance intensity maps for different
-        return periods. Calls local_exceedance_inten.
+        return periods. Calls local_exceedance_intensity. For handling large data sets and for
+        further options, see Notes.
 
         Parameters
         ----------
         return_periods: tuple(int), optional
             return periods to consider
-        smooth: bool, optional
-            smooth plot to plot.RESOLUTIONxplot.RESOLUTION
         axis: matplotlib.axes._subplots.AxesSubplot, optional
             axis to use
-        figsize: tuple, optional
-            figure size for plt.subplots
+        kwargs_local_exceedance_intensity: dict
+            Dictionary of keyword arguments for the method hazard.local_exceedance_intensity.
         kwargs: optional
             arguments for pcolormesh matplotlib function used in event plots
 
         Returns
         -------
         axis, inten_stats:  matplotlib.axes._subplots.AxesSubplot, np.ndarray
             intenstats is return_periods.size x num_centroids
+
+        See Also
+        ---------
+        hazard.local_exceedance_intensity: method to calculate local exceedance frequencies.
+
+         Notes
+         -----
+         For handling large data, and for more fleixble options in the exceedance
+         intensity computation and in the plotting, we recommend to use
+         gdf, title, labels = hazard.local_exceedance_intensity() and
+         util.plot.plot_from_gdf(gdf, title, labels) instead.
         """
-        inten_stats = self.local_exceedance_intensity(return_periods)[0].values[:, 1:].T
-        inten_stats = inten_stats.astype(float)
-        colbar_name = "Intensity (" + self.units + ")"
-        title = list()
-        for ret in return_periods:
-            title.append("Return period: " + str(ret) + " years")
-        axis = u_plot.geo_im_from_array(
-            inten_stats,
-            self.centroids.coord,
-            colbar_name,
-            title,
-            smooth=smooth,
-            axes=axis,
-            figsize=figsize,
-            adapt_fontsize=adapt_fontsize,
-            **kwargs,
+        LOGGER.info(
+            "Some errors in the previous calculation of local exceedance intensities have been corrected,"
+            " see Hazard.local_exceedance_intensity. To reproduce data with the "
+            "previous calculation, use CLIMADA v5.0.0 or less."
+        )
+
+        if kwargs_local_exceedance_intensity is None:
+            kwargs_local_exceedance_intensity = {}
+
+        inten_stats, title, column_labels = self.local_exceedance_intensity(
+            return_periods, **kwargs_local_exceedance_intensity
+        )
+
+        axis = u_plot.plot_from_gdf(
+            inten_stats, title, column_labels, axis=axis, **kwargs
         )
-        return axis, inten_stats
+        return axis, inten_stats.values[:, 1:].T.astype(float)
 
     def plot_intensity(
         self,