Implement quantiles and median for hazard and impact forecasts (#1191)

* Add quantile and median for hazard and impact forecast
* Review

---------

Co-authored-by: luseverin <[email protected]>
Co-authored-by: Chahan Kropf <[email protected]>
Co-authored-by: Lukas Riedel <[email protected]>

Author: 4 people

climada/engine/impact_forecast.py

@@ -412,3 +412,56 @@ def select(
             coord_exp=coord_exp,
             reset_frequency=reset_frequency,
         )
+
+    def _quantile(self, q: float, event_name: str | None = None):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the quantile value.
+        """
+        red_imp_mat = sparse.csr_matrix(np.quantile(self.imp_mat.toarray(), q, axis=0))
+        red_at_event = np.array([red_imp_mat.sum()])
+        if event_name is None:
+            event_name = f"quantile_{q}"
+        return ImpactForecast(
+            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,
+            **self._reduce_attrs(event_name),
+        )
+
+    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.
+        """
+        return self._quantile(q=q)
+
+    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(q=0.5, event_name="median")

climada/engine/test/test_impact_forecast.py

@@ -298,3 +298,49 @@ def test_impact_forecast_min_mean_max(impact_forecast_stats, attr):
     npt.assert_array_equal(imp_fc_reduced.event_id, [0])
     npt.assert_array_equal(imp_fc_reduced.frequency, [1])
     npt.assert_array_equal(imp_fc_reduced.date, [0])
+
+
+@pytest.mark.parametrize("quantile", [0.3, 0.6, 0.8])
+def test_impact_forecast_quantile(impact_forecast, quantile):
+    """Check quantile method for ImpactForecast"""
+    imp_fcst_quantile = impact_forecast.quantile(q=quantile)
+
+    # assert imp_mat
+    npt.assert_array_equal(
+        imp_fcst_quantile.imp_mat.toarray().squeeze(),
+        np.quantile(impact_forecast.imp_mat.toarray(), quantile, axis=0),
+    )
+    # assert at_event
+    npt.assert_array_equal(
+        imp_fcst_quantile.at_event,
+        np.quantile(impact_forecast.at_event, quantile, 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_{quantile}"])
+    )
+    npt.assert_array_equal(imp_fcst_quantile.frequency, np.array([1]))
+    npt.assert_array_equal(imp_fcst_quantile.date, np.array([0]))
+
+
+def test_median(impact_forecast):
+    imp_fcst_median = impact_forecast.median()
+    imp_fcst_quantile = impact_forecast.quantile(q=0.5)
+    npt.assert_array_equal(
+        imp_fcst_median.imp_mat.toarray(), imp_fcst_quantile.imp_mat.toarray()
+    )
+    npt.assert_array_equal(imp_fcst_median.imp_mat.toarray(), [[2.5, 2.5]])
+
+    # check that attributes where reduced correctly
+    npt.assert_array_equal(imp_fcst_median.member, np.array([-1]))
+    npt.assert_array_equal(imp_fcst_median.lead_time, np.array([np.timedelta64("NaT")]))
+    npt.assert_array_equal(imp_fcst_median.event_id, np.array([0]))
+    npt.assert_array_equal(imp_fcst_median.event_name, np.array(["median"]))
+    npt.assert_array_equal(imp_fcst_median.frequency, np.array([1]))
+    npt.assert_array_equal(imp_fcst_median.date, np.array([0]))

climada/hazard/forecast.py

@@ -281,3 +281,60 @@ def select(
             extent=extent,
             reset_frequency=reset_frequency,
         )
+
+    def _quantile(self, q: float, event_name: str | None = None):
+        """
+        Reduce the impact matrix and at_event of a HazardForecast to the quantile value.
+        """
+        red_intensity = sparse.csr_matrix(
+            np.quantile(self.intensity.toarray(), q, axis=0)
+        )
+        red_fraction = sparse.csr_matrix(
+            np.quantile(self.fraction.toarray(), q, axis=0)
+        )
+        if event_name is None:
+            event_name = f"quantile_{q}"
+        return HazardForecast(
+            haz_type=self.haz_type,
+            pool=self.pool,
+            units=self.units,
+            centroids=self.centroids,
+            frequency_unit=self.frequency_unit,
+            intensity=red_intensity,
+            fraction=red_fraction,
+            **self._reduce_attrs(event_name),
+        )
+
+    def quantile(self, q: float):
+        """
+        Reduce the impact matrix and at_event of a HazardForecast to the quantile value.
+
+        The quantile value is computed by taking the quantile of the impact matrix
+        along the event dimension axis (axis=0) and then taking the quantile of the
+        resulting array.
+
+        Parameters
+        ----------
+        q : float
+            The quantile to compute, between 0 and 1.
+
+        Returns
+        -------
+        HazardForecast
+            A HazardForecast object with the quantile intensity and fraction.
+        """
+        return self._quantile(q=q)
+
+    def median(self):
+        """
+        Reduce the impact matrix and at_event of a HazardForecast to the median value.
+
+        The median value is computed by taking the median of the impact matrix along the
+        event dimension axis (axis=0) and then taking the median of the resulting array.
+
+        Returns
+        -------
+        HazardForecast
+            A HazardForecast object with the median intensity and fraction.
+        """
+        return self._quantile(q=0.5, event_name="median")

climada/hazard/test/test_forecast.py

@@ -258,3 +258,61 @@ def test_hazard_forecast_mean_min_max(haz_fc, attr):
     npt.assert_array_equal(haz_fcst_reduced.frequency, [1])
     npt.assert_array_equal(haz_fcst_reduced.date, [0])
     npt.assert_array_equal(haz_fcst_reduced.orig, [True])
+
+
+@pytest.mark.parametrize("quantile", [0.3, 0.6, 0.8])
+def test_hazard_forecast_quantile(haz_fc, quantile):
+    """Check quantile method for HazardForecast"""
+    haz_fcst_quantile = haz_fc.quantile(q=quantile)
+
+    # assert intensity
+    npt.assert_array_equal(
+        haz_fcst_quantile.intensity.toarray().squeeze(),
+        np.quantile(haz_fc.intensity.toarray(), quantile, axis=0),
+    )
+    # assert fraction
+    npt.assert_array_equal(
+        haz_fcst_quantile.fraction.toarray().squeeze(),
+        np.quantile(haz_fc.fraction.toarray(), quantile, axis=0),
+    )
+
+    # check that attributes where reduced correctly
+    npt.assert_array_equal(
+        haz_fcst_quantile.lead_time, np.array([np.timedelta64("NaT")])
+    )
+    npt.assert_array_equal(haz_fcst_quantile.member, np.array([-1]))
+    npt.assert_array_equal(
+        haz_fcst_quantile.event_name, np.array([f"quantile_{quantile}"])
+    )
+    npt.assert_array_equal(haz_fcst_quantile.event_id, np.array([0]))
+    npt.assert_array_equal(haz_fcst_quantile.frequency, np.array([1]))
+    npt.assert_array_equal(haz_fcst_quantile.date, np.array([0]))
+    npt.assert_array_equal(haz_fcst_quantile.orig, np.array([True]))
+
+
+def test_median(haz_fc):
+    haz_fcst_median = haz_fc.median()
+    haz_fcst_quantile = haz_fc.quantile(q=0.5)
+    npt.assert_array_equal(
+        haz_fcst_median.intensity.todense(), haz_fcst_quantile.intensity.todense()
+    )
+    npt.assert_array_equal(
+        haz_fcst_median.intensity.todense(),
+        np.median(haz_fc.intensity.todense(), axis=0),
+    )
+    npt.assert_array_equal(
+        haz_fcst_median.fraction.todense(), haz_fcst_quantile.fraction.todense()
+    )
+    npt.assert_array_equal(
+        haz_fcst_median.fraction.todense(),
+        np.median(haz_fc.fraction.todense(), axis=0),
+    )
+
+    # check that attributes where reduced correctly
+    npt.assert_array_equal(haz_fcst_median.member, np.array([-1]))
+    npt.assert_array_equal(haz_fcst_median.lead_time, np.array([np.timedelta64("NaT")]))
+    npt.assert_array_equal(haz_fcst_median.event_id, np.array([0]))
+    npt.assert_array_equal(haz_fcst_median.event_name, np.array(["median"]))
+    npt.assert_array_equal(haz_fcst_median.frequency, np.array([1]))
+    npt.assert_array_equal(haz_fcst_median.date, np.array([0]))
+    npt.assert_array_equal(haz_fcst_median.orig, np.array([True]))