refactorize methods by introducing new interpolate wrapper function and new util function

Author: ValentinGebhart

climada/engine/impact.py

@@ -53,6 +53,7 @@
 import climada.util.plot as u_plot
 from climada.util.select import get_attributes_with_matching_dimension
 import climada.util.interpolation as u_interp
+import climada.util.checker as u_check
 
 LOGGER = logging.getLogger(__name__)
 
@@ -523,55 +524,32 @@ def local_exceedance_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
-        if self.frequency_unit in ['1/year', 'annual', '1/y', '1/a']:
-            return_period_unit = 'years'
-        elif self.frequency_unit in ['1/month', 'monthly', '1/m']:
-            return_period_unit = 'months'
-        elif self.frequency_unit in ['1/week', 'weekly', '1/w']:
-            return_period_unit = 'weeks'
-        else:
-            LOGGER.warning("Hazard's frequency unit %s is not known, "
-                            "years will be used as return period unit.",
-                            self.frequency_unit)
-            return_period_unit = 'years'
-
-        num_centroids = self.imp_mat.shape[1]
-        imp_stats = np.zeros((len(return_periods), num_centroids))
-
-        for i in range(num_centroids):
-            # sort intensties and frequencies at given centroid
-            impact = np.squeeze(self.imp_mat[:,i].toarray())
-            sorted_idxs = np.argsort(impact)
-            impact = np.squeeze(impact[sorted_idxs])
-            frequency = self.frequency[sorted_idxs]
-
-            # group values with same impact
-            frequency, impact = u_interp.group_frequency(frequency, impact)
-
-            # fit intensities to cummulative frequencies
-            frequency = np.cumsum(frequency[::-1])[::-1]
-            if method == 'stepfunction':
-                imp_stats[:,i] = u_interp.stepfunction_ev(
-                    1/np.array(return_periods), frequency[::-1], impact[::-1],
-                    y_threshold=min_impact, y_asymptotic=0.
-                )
-            elif method in ['interpolate', 'extrapolate', 'extrapolate_constant']:
-                extrapolation = None if method == 'interpolate' else method
-                imp_stats[:,i] = u_interp.interpolate_ev(
-                    1/np.array(return_periods), frequency[::-1], impact[::-1],
-                    logx=log_frequency, logy=log_impact, y_threshold=min_impact,
-                    extrapolation=extrapolation, y_asymptotic=0.
-                )
-            else:
-                raise ValueError(f"Unknown method: {method}")
+        return_period_unit = u_check.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 exceedance impact
+        test_frequency = 1/np.array(return_periods)
+        exceedance_impact = np.array([
+            u_interp.preprocess_and_interpolate_ev(
+                test_frequency, None, 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=0.
+            )
+            for i_centroid in range(self.imp_mat.shape[1])
+        ])
 
         # 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'{ret_per}' for ret_per in return_periods]
-        gdf[col_names] = imp_stats.T
+        gdf[col_names] = exceedance_impact
         # create label and column_label
         label = f'Impact ({self.unit})'
         column_label = lambda column_names: [f'Return Period: {col} {return_period_unit}'

climada/hazard/base.py

@@ -453,7 +453,7 @@ def local_exceedance_intensity(
             return_periods=(25, 50, 100, 250),
             method='interpolate',
             min_intensity=None,
-            log_frequeny=True,
+            log_frequency=True,
             log_intensity=True
     ):
         """Compute local exceedance intensity for given return periods. The default method
@@ -502,51 +502,28 @@ def local_exceedance_intensity(
         if not min_intensity and min_intensity != 0:
             min_intensity = self.intensity_thres
         #check frequency unit
-        if self.frequency_unit in ['1/year', 'annual', '1/y', '1/a']:
-            return_period_unit = 'years'
-        elif self.frequency_unit in ['1/month', 'monthly', '1/m']:
-            return_period_unit = 'months'
-        elif self.frequency_unit in ['1/week', 'weekly', '1/w']:
-            return_period_unit = 'weeks'
-        else:
-            LOGGER.warning(f"Hazard's frequency unit {self.frequency_unit} is not known, "
-                            "years will be used as return period unit.")
-            return_period_unit = 'years'
-        
-        num_centroids = self.intensity.shape[1]
-        inten_stats = np.zeros((len(return_periods), num_centroids))
-
-        for i in range(num_centroids):
-            # sort intensties and frequencies at given centroid
-            intensity = np.squeeze(self.intensity[:,i].toarray())
-            sorted_idxs = np.argsort(intensity)
-            intensity = np.squeeze(intensity[sorted_idxs])
-            frequency = self.frequency[sorted_idxs]
-
-            # group values with same intensity
-            frequency, intensity = u_interp.group_frequency(frequency, intensity)
-
-            # fit intensities to cummulative frequencies
-            frequency = np.cumsum(frequency[::-1])[::-1]
-            if method == 'stepfunction':
-                inten_stats[:,i] = u_interp.stepfunction_ev(
-                    1/np.array(return_periods), frequency[::-1], intensity[::-1],
-                    y_threshold=min_intensity, y_asymptotic=0.
-                )
-            elif method in ['interpolate', 'extrapolate', 'extrapolate_constant']:
-                extrapolation = None if method == 'interpolate' else method
-                inten_stats[:,i] = u_interp.interpolate_ev(
-                    1/np.array(return_periods), frequency[::-1], intensity[::-1],
-                    logx=log_frequeny, logy=log_intensity, y_threshold=min_intensity,
-                    extrapolation=extrapolation, y_asymptotic=0.
-                )
-            else:
-                raise ValueError(f"Unknown method: {method}")
+        return_period_unit = u_check.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 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.
+            )
+            for i_centroid in range(self.intensity.shape[1])
+        ])
 
         # create the output GeoDataFrame
         gdf = gpd.GeoDataFrame(geometry = self.centroids.gdf['geometry'], crs = self.centroids.gdf.crs)
         column_names = [f'{rp}' for rp in return_periods]
-        gdf[column_names] = inten_stats.T
+        gdf[column_names] = exceedance_intensity
 
         # create label and column_label
         label = f'Intensity ({self.units})'
@@ -624,55 +601,29 @@ def local_return_period(
         if not min_intensity and min_intensity != 0:
             min_intensity = self.intensity_thres
         #check frequency unit
-        if self.frequency_unit in ['1/year', 'annual', '1/y', '1/a']:
-            return_period_unit = 'Years'
-        elif self.frequency_unit in ['1/month', 'monthly', '1/m']:
-            return_period_unit = 'Months'
-        elif self.frequency_unit in ['1/week', 'weekly', '1/w']:
-            return_period_unit = 'Weeks'
-        else:
-            LOGGER.warning("Hazard's frequency unit %s is not known, "
-                           "years will be used as return period unit.", self.frequency_unit)
-            return_period_unit = 'Years'
-
-        # Ensure threshold_intensities is a numpy array
-        threshold_intensities = np.array(threshold_intensities)
-        
-        num_centroids = self.intensity.shape[1]
-        return_periods = np.zeros((len(threshold_intensities), num_centroids))
-        
-        for i in range(num_centroids):
-            # sort intensties and frequencies at given centroid
-            intensity = np.squeeze(self.intensity[:,i].toarray())
-            sorted_idxs = np.argsort(intensity)
-            intensity = np.squeeze(intensity[sorted_idxs])
-            frequency = self.frequency[sorted_idxs]
-
-            # group values with same intensity
-            frequency, intensity = u_interp.group_frequency(frequency, intensity)
-
-            # fit intensities to cummulative frequencies
-            frequency = np.cumsum(frequency[::-1])[::-1]
-            if method == 'stepfunction':
-                return_periods[:,i] = u_interp.stepfunction_ev(
-                    threshold_intensities, intensity, frequency, x_threshold=min_intensity
-                )
-            elif method in ['interpolate', 'extrapolate', 'extrapolate_constant']:
-                extrapolation = None if method == 'interpolate' else method
-                return_periods[:,i] = u_interp.interpolate_ev(
-                    threshold_intensities, intensity, frequency, logx=log_intensity,
-                    logy=log_frequency, x_threshold=min_intensity, extrapolation=extrapolation,
-                    y_asymptotic=np.nan
-                )
-            else:
-                raise ValueError(f"Unknown method: {method}")
+        return_period_unit = u_check.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_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 = safe_divide(1., return_periods)
-        
+
         # create the output GeoDataFrame
         gdf = gpd.GeoDataFrame(geometry = self.centroids.gdf['geometry'],
                                crs = self.centroids.gdf.crs)
         col_names = [f'{tresh_inten}' for tresh_inten in threshold_intensities]
-        gdf[col_names] = return_periods.T
+        gdf[col_names] = return_periods
 
         # create label and column_label
         label = f'Return Periods ({return_period_unit})'

climada/hazard/test/test_base.py

@@ -1049,7 +1049,7 @@ def test_local_exceedance_intensity(self):
         # third centroid has intensities 1 with cum frequencies 1
         # testing at frequencies 2, 1.5, 1
         inten_stats, _, _ = haz.local_exceedance_intensity(
-            return_period, log_frequeny=False, log_intensity=False, method='extrapolate_constant')
+            return_period, log_frequency=False, log_intensity=False, method='extrapolate_constant')
         np.testing.assert_allclose(
             inten_stats[inten_stats.columns[1:]].values,
             np.array([

climada/test/test_hazard.py

@@ -458,7 +458,7 @@ def test_local_exceedance_intensity_methods(self):
 
         # test lin lin interpolation without extrapolation
         inten_stats, _, _ = haz.local_exceedance_intensity(
-                return_periods=(1000, 30, .1), log_frequeny=False, log_intensity=False,
+                return_periods=(1000, 30, .1), log_frequency=False, log_intensity=False,
                 method='extrapolate_constant')
         np.testing.assert_allclose(
             inten_stats.values[:,1:].astype(float),

climada/util/checker.py

@@ -214,3 +214,29 @@ def prune_csr_matrix(matrix: sparse.csr_matrix):
     matrix.check_format()
     matrix.eliminate_zeros()
     matrix.sum_duplicates()
+
+def convert_frequency_unit_to_time_unit(frequency_unit):
+    """Converts common frequency units to corresponding time units. Unknown frequency
+    units are converted to "years".
+
+    Parameters
+    ----------
+    frequency_unit : str
+        unit of frequency
+
+    Returns
+    -------
+    str
+        corresponding time unit.
+    """
+    if frequency_unit in ['1/year', 'annual', '1/y', '1/a']:
+        time_unit = 'years'
+    elif frequency_unit in ['1/month', 'monthly', '1/m']:
+        time_unit = 'months'
+    elif frequency_unit in ['1/week', 'weekly', '1/w']:
+        time_unit = 'weeks'
+    else:
+        LOGGER.warning(f"Frequency unit {frequency_unit} is not known, "
+                        "years will be used as time unit.")
+        time_unit = 'years'
+    return time_unit