added method "extrapolate_constant" and changed behaviour of "interpolate"

ValentinGebhart · ValentinGebhart · commit 55c033f50c12 · 2024-09-28T12:49:39.000+02:00
diff --git a/climada/engine/impact.py b/climada/engine/impact.py
@@ -488,7 +488,9 @@ def local_exceedance_impact(
             locally (at each centroid). Defaults to (25, 50, 100, 250).
         method : str
             Method to interpolate to new return periods. Currently available are "interpolate",
-            "extrapolate" and "stepfunction". If set to "interpolate" or "stepfunction",
+            "extrapolate", "extrapolate_constant" and "stepfunction". If set to "interpolate",
+            return periods outside the range of the Impact object's observed local return periods
+            will be assigned NaN. If set to "extrapolate_constant" or "stepfunction",
             return periods larger than the Impact object's observed local return periods will be
             assigned the largest local impact, and return periods smaller than the Impact object's
             observed local return periods will be assigned 0. If set to "extrapolate", local
@@ -551,14 +553,15 @@ def local_exceedance_impact(
             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.
+                    1/np.array(return_periods), frequency[::-1], impact[::-1],
+                    y_threshold=min_impact, y_asymptotic=0.
                 )
-            elif method == 'extrapolate' or method == 'interpolate':
-                extrapolation = (method == 'extrapolate')
+            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.
+                    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}")
diff --git a/climada/hazard/base.py b/climada/hazard/base.py
@@ -467,7 +467,9 @@ def local_exceedance_intensity(
             locally (at each centroid). Defaults to (25, 50, 100, 250).
         method : str
             Method to interpolate to new return periods. Currently available are "interpolate",
-            "extrapolate" and "stepfunction". If set to "interpolate" or "stepfunction", return
+            "extrapolate", "extrapolate_constant" and "stepfunction". If set to "interpolate",
+            return periods outside the range of the Hazard object's observed local return periods
+            will be assigned NaN. If set to "extrapolate_constant" or "stepfunction", return
             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
@@ -528,14 +530,15 @@ def local_exceedance_intensity(
             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.
+                    1/np.array(return_periods), frequency[::-1], intensity[::-1],
+                    y_threshold=min_intensity, y_asymptotic=0.
                 )
-            elif method == 'interpolate' or method == 'extrapolate':
-                extrapolation = (method == 'extrapolate')
+            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, y_asymptotic=0., extrapolation=extrapolation
+                    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}")
@@ -587,7 +590,9 @@ def local_return_period(
             locally (at each centroid). Defaults to (10, 20)
         method : str
             Method to interpolate to new threshold intensities. Currently available are
-            "interpolate", "extrapolate" and "stepfunction". If set to "interpolate" or
+            "interpolate", "extrapolate", "extrapolate_constant" and "stepfunction". If set to
+            "interpolate", threshold intensities outside the range of the Hazard object's local
+            intensities will be assigned NaN. If set to "extrapolate_constant" or
             "stepfunction", threshold intensities larger than the Hazard object's local
             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.
@@ -653,11 +658,12 @@ def local_return_period(
                 return_periods[:,i] = u_interp.stepfunction_ev(
                     threshold_intensities, intensity, frequency, x_threshold=min_intensity
                 )
-            elif method == 'interpolate' or method == "extrapolate":
-                extrapolation = (method == "extrapolate")
+            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
+                    logy=log_frequency, x_threshold=min_intensity, extrapolation=extrapolation,
+                    y_asymptotic=np.nan
                 )
             else:
                 raise ValueError(f"Unknown method: {method}")
diff --git a/climada/hazard/test/test_base.py b/climada/hazard/test/test_base.py
@@ -1032,7 +1032,7 @@ def test_degenerate_pass(self):
         return_period = np.array([25, 50, 100, 250])
         haz.intensity = sparse.csr_matrix(np.zeros(haz.intensity.shape))
         inten_stats = haz.local_exceedance_intensity(return_period)[0].values[:, 1:].T.astype(float)
-        self.assertTrue(np.array_equal(inten_stats, np.zeros((4, 100))))
+        np.testing.assert_allclose(inten_stats, np.full((4, 100), np.nan))
 
     def test_local_exceedance_intensity(self):
         """Test local exceedance frequencies with lin lin interpolation"""
@@ -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)
+            return_period, log_frequeny=False, log_intensity=False, method='extrapolate_constant')
         np.testing.assert_allclose(
             inten_stats[inten_stats.columns[1:]].values,
             np.array([
@@ -1073,7 +1073,8 @@ def test_local_return_period(self):
         # third centroid has intensities 1 with cum frequencies 1 (0 intensity is neglected)
         # testing at intensities 1, 2, 3
         return_stats, _, _ = haz.local_return_period(
-            threshold_intensities, log_frequency=False, log_intensity=False, min_intensity=0)
+            threshold_intensities, log_frequency=False, log_intensity=False, min_intensity=0,
+            method='extrapolate_constant')
         np.testing.assert_allclose(
             return_stats[return_stats.columns[1:]].values,
             np.array([
diff --git a/climada/test/test_engine.py b/climada/test/test_engine.py
@@ -335,9 +335,9 @@ def test_local_exceedance_impact_methods(self):
                 ]),
             rtol=0.8)
 
-        # test log log interpolation and no extrapolation
+        # test log log interpolation and extrapolation with constant
         impact_stats, _, _ = impact.local_exceedance_impact(
-                return_periods=(1000, 30, .1))
+                return_periods=(1000, 30, .1), method = "extrapolate_constant")
         np.testing.assert_allclose(
             impact_stats.values[:,1:].astype(float),
             np.array([
@@ -348,9 +348,22 @@ def test_local_exceedance_impact_methods(self):
                 ]),
             rtol=0.8)
 
+        # test log log interpolation and no extrapolation
+        impact_stats, _, _ = impact.local_exceedance_impact(
+                return_periods=(1000, 30, .1))
+        np.testing.assert_allclose(
+            impact_stats.values[:,1:].astype(float),
+            np.array([
+                [np.nan, np.nan, np.nan],
+                [np.nan, np.nan, np.nan],
+                [np.nan, 1e2, np.nan],
+                [np.nan, 300, np.nan]
+                ]),
+            rtol=0.8)
+
         # test lin lin interpolation with no extrapolation
         impact_stats, _, _ = impact.local_exceedance_impact(
-                return_periods=(1000, 30, .1),
+                return_periods=(1000, 30, .1), method = "extrapolate_constant",
                 log_frequency=False, log_impact=False)
         np.testing.assert_allclose(
             impact_stats.values[:,1:].astype(float),
diff --git a/climada/test/test_hazard.py b/climada/test/test_hazard.py
@@ -432,9 +432,9 @@ def test_local_exceedance_intensity_methods(self):
                 ]),
             rtol=0.8)
 
-        # test log log interpolation and no extrapolation
+        # test log log interpolation and extrapolation with constant
         inten_stats, _, _ = haz.local_exceedance_intensity(
-                return_periods=(1000, 30, .1))
+                return_periods=(1000, 30, .1), method='extrapolate_constant')
         np.testing.assert_allclose(
             inten_stats.values[:,1:].astype(float),
             np.array([
@@ -444,9 +444,22 @@ def test_local_exceedance_intensity_methods(self):
                 ]),
             rtol=0.8)
 
+        # test log log interpolation and no extrapolation
+        inten_stats, _, _ = haz.local_exceedance_intensity(
+                return_periods=(1000, 30, .1))
+        np.testing.assert_allclose(
+            inten_stats.values[:,1:].astype(float),
+            np.array([
+                [np.nan, np.nan, np.nan],
+                [np.nan, 1e2, np.nan],
+                [np.nan, 300, np.nan]
+                ]),
+            rtol=0.8)
+
         # 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_frequeny=False, log_intensity=False,
+                method='extrapolate_constant')
         np.testing.assert_allclose(
             inten_stats.values[:,1:].astype(float),
             np.array([
@@ -495,9 +508,9 @@ def test_local_return_period_methods(self):
                 ]),
             rtol=0.8)
 
-        # test log log interpolation and no extrapolation
+        # test log log interpolation and extrapolation with constant
         return_stats, _, _ = haz.local_return_period(
-                threshold_intensities=(.1, 300, 1e5))
+                threshold_intensities=(.1, 300, 1e5), method='extrapolate_constant')
         np.testing.assert_allclose(
             return_stats.values[:,1:].astype(float),
             np.array([
@@ -507,6 +520,18 @@ def test_local_return_period_methods(self):
                 ]),
         rtol=0.8)
 
+        # test log log interpolation and no extrapolation
+        return_stats, _, _ = haz.local_return_period(
+                threshold_intensities=(.1, 300, 1e5))
+        np.testing.assert_allclose(
+            return_stats.values[:,1:].astype(float),
+            np.array([
+                [np.nan, np.nan, np.nan],
+                [np.nan, 30, np.nan],
+                [np.nan, 30, np.nan]
+                ]),
+        rtol=0.8)
+
 # Execute Tests
 if __name__ == "__main__":
     TESTS = unittest.TestLoader().loadTestsFromTestCase(TestCentroids)
