Add quantile and median for impact forecast

Author: luseverin

climada/engine/impact_forecast.py

@@ -376,3 +376,54 @@ def mean(self):
             imp_mat=red_imp_mat,
             haz_type=self.haz_type,
         )
+
+    def quantile(self, q: float):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the quantile value.
+
+        Parameters
+        ----------
+        q : float
+            The quantile to compute, which must be between 0 and 1.
+
+        Returns
+        -------
+        ImpactForecast
+            An ImpactForecast object with the quantile impact matrix and at_event.
+        """
+        red_imp_mat = sparse.csr_matrix(np.quantile(self.imp_mat.toarray(), q, axis=0))
+        red_at_event = np.array([red_imp_mat.sum()])
+        reduced_attrs = self._reduce_attrs(f"quantile_{q}")
+        return ImpactForecast(
+            lead_time=reduced_attrs["lead_time"],
+            member=reduced_attrs["member"],
+            event_id=reduced_attrs["event_id"],
+            event_name=reduced_attrs["event_name"],
+            date=reduced_attrs["date"],
+            frequency=reduced_attrs["frequency"],
+            frequency_unit=self.frequency_unit,
+            coord_exp=self.coord_exp,
+            crs=self.crs,
+            eai_exp=self.eai_exp,
+            at_event=red_at_event,
+            tot_value=self.tot_value,
+            aai_agg=self.aai_agg,
+            unit=self.unit,
+            imp_mat=red_imp_mat,
+            haz_type=self.haz_type,
+        )
+
+    def median(self):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the median value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        ImpactForecast
+            An ImpactForecast object with the median impact matrix and at_event.
+        """
+        return self.quantile(0.5)

climada/engine/test/test_impact_forecast.py

@@ -261,3 +261,32 @@ def test_impact_forecast_mean_min_max(impact_forecast):
     assert imp_fcst_mean.date == 0
     assert imp_fcst_min.date == 0
     assert imp_fcst_max.date == 0
+
+
+def test_impact_forecast_quantile(impact_forecast):
+    """Check quantile method for ImpactForecast"""
+    for q in [0.0, 0.5, 0.8]:
+        imp_fcst_quantile = impact_forecast.quantile(q)
+
+        # assert imp_mat
+        npt.assert_array_equal(
+            imp_fcst_quantile.imp_mat.toarray().squeeze(),
+            np.quantile(impact_forecast.imp_mat.toarray(), q, axis=0),
+        )
+        # assert at_event
+        npt.assert_array_equal(
+            imp_fcst_quantile.at_event,
+            np.quantile(impact_forecast.at_event, q, axis=0).sum(),
+        )
+
+        # check that attributes where reduced correctly
+        npt.assert_array_equal(imp_fcst_quantile.member, np.array([-1]))
+        npt.assert_array_equal(
+            imp_fcst_quantile.lead_time, np.array([np.timedelta64("NaT")])
+        )
+        npt.assert_array_equal(imp_fcst_quantile.event_id, np.array([0]))
+        npt.assert_array_equal(
+            imp_fcst_quantile.event_name, np.array([f"quantile_{q}"])
+        )
+        npt.assert_array_equal(imp_fcst_quantile.frequency, np.array([1]))
+        npt.assert_array_equal(imp_fcst_quantile.date, np.array([0]))