Add polynomial s shape impact function

Author: Chahan Kropf

climada/entity/impact_funcs/base.py

@@ -273,3 +273,70 @@ def set_sigmoid_impf(self, *args, **kwargs):
         LOGGER.warning("The use of ImpactFunc.set_sigmoid_impf is deprecated."
                        "Use ImpactFunc.from_sigmoid_impf instead.")
         self.__dict__ = ImpactFunc.from_sigmoid_impf(*args, **kwargs).__dict__
+
+    @classmethod
+    def from_poly_s_shape(
+        cls,
+        intensity: tuple[float, float, float],
+        threshold: float,
+        half_point: float,
+        upper_limit: float,
+        exponent: float,
+        haz_type: str,
+        impf_id: int = 1,
+        **kwargs):
+        """S-shape polynomial impact function hinging on four parameter.
+
+        f(I) = D(I)**exponent / (1 + D(I)**exponent) * upper_limit
+
+        D(I) = max[I - threshold, 0] / (half_point - threshold)
+
+        This function is inspired by Emanuel et al. (2011)
+        https://doi.org/10.1175/WCAS-D-11-00007.1
+
+        Parameters
+        ----------
+        intensity: tuple(float, float, float)
+            tuple of 3 intensity numbers along np.arange(min, max, step)
+        threshold : float
+            Intensity threshold below which there is no impact.
+            Should in general be larger than 0 for computational efficiency
+            of impacts.
+        half_point : float
+            Intensity at which 50% of maxixmum impact is expected.
+            Must be larger than thres.
+        upper_limit : float
+            Maximum impact value. Must be larger than 0.
+        exponent: float
+            Exponent of the polynomial. Must be larger than 0.
+        haz_type: str
+            Reference string for the hazard (e.g., 'TC', 'RF', 'WS', ...)
+        impf_id : int, optional, default=1
+            Impact function id
+        kwargs :
+            keyword arguments passed to ImpactFunc()
+
+        Returns
+        -------
+        impf : climada.entity.impact_funcs.ImpactFunc
+            s-shapep polynomial impact function
+        """
+        intensity = np.arange(*intensity)
+
+        if threshold > half_point:
+            mdd = np.zeros_like(intensity)
+        else:
+            D = (intensity - threshold) / (half_point - threshold)
+            D[D < 0] = 0
+            mdd = D**exponent / (1 + D**exponent) * upper_limit
+        paa = np.ones_like(intensity)
+
+        impf = cls(
+            haz_type=haz_type,
+            id=impf_id,
+            intensity=intensity,
+            paa=paa,
+            mdd=mdd,
+            **kwargs
+        )
+        return impf

climada/entity/impact_funcs/test/test_base.py

@@ -47,7 +47,6 @@ def test_from_step(self):
         self.assertEqual(imp_fun.haz_type, 'TC')
         self.assertEqual(imp_fun.id, 2)
 
-
     def test_from_sigmoid(self):
         """Check default impact function: sigmoid function"""
         inten = (0, 100, 5)
@@ -60,6 +59,34 @@ def test_from_sigmoid(self):
         self.assertEqual(imp_fun.haz_type, 'RF')
         self.assertEqual(imp_fun.id, 2)
 
+    def test_from_poly_s_shape(self):
+        """Check default impact function: polynomial s-shape"""
+        inten = (0, 5, 1)
+        impf = ImpactFunc.from_poly_s_shape(
+            inten, threshold=0.2, half_point=1, upper_limit=0.8, exponent=4,
+            haz_type='RF', impf_id=2, intensity_unit='m'
+            )
+        correct_mdd = np.array([0, 0.4, 0.76995746, 0.79470418, 0.79843158])
+        self.assertTrue(np.array_equal(impf.paa, np.ones(5)))
+        np.testing.assert_array_almost_equal(impf.mdd, correct_mdd)
+        self.assertTrue(np.array_equal(impf.intensity, np.arange(0, 5, 1)))
+        self.assertEqual(impf.haz_type, 'RF')
+        self.assertEqual(impf.id, 2)
+        self.assertEqual(impf.intensity_unit, 'm')
+
+        # If threshold > half_point, mdd should all be 0
+        impf = ImpactFunc.from_poly_s_shape(
+            inten, threshold=2, half_point=1, upper_limit=0.8, exponent=4,
+            haz_type='RF', impf_id=2, intensity_unit='m'
+            )
+        correct_mdd = np.array([0, 0.4, 0.76995746, 0.79470418, 0.79843158])
+        self.assertTrue(np.array_equal(impf.paa, np.ones(5)))
+        np.testing.assert_array_almost_equal(impf.mdd, np.zeros(5))
+        self.assertTrue(np.array_equal(impf.intensity, np.arange(0, 5, 1)))
+        self.assertEqual(impf.haz_type, 'RF')
+        self.assertEqual(impf.id, 2)
+        self.assertEqual(impf.intensity_unit, 'm')
+
 # Execute Tests
 if __name__ == "__main__":
     TESTS = unittest.TestLoader().loadTestsFromTestCase(TestInterpolation)