Implement quantiles and median for hazard and impact forecasts #1191

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged

peanutfun merged 44 commits into forecast-class from implement_quantiles

Dec 10, 2025

+214 −0

climada/engine/impact_forecast.py

-Original file line number
+Diff line change
@@ Expand Up / @@ -354,3 +354,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

-Original file line number
+Diff line change
@@ Expand Up @@
         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

-Original file line number
+Diff line change
@@ Expand Up / @@ -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.
-        Reduce the impact matrix and at_event of a HazardForecast to the median value.
+        Reduce the impact matrix and at_event of a HazardForecast to the median value over all events.
+            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

-Original file line number
+Diff line change
@@ Expand Up / @@ -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]))

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implement quantiles and median for hazard and impact forecasts #1191

Diff view

Diff view

There are no files selected for viewing

chahank Dec 10, 2025

Uh oh!

Uh oh!

Implement quantiles and median for hazard and impact forecasts #1191

Implement quantiles and median for hazard and impact forecasts #1191

Uh oh!

Uh oh!

Diff view

Diff view

There are no files selected for viewing

chahank Dec 10, 2025

Choose a reason for hiding this comment

Uh oh!

Uh oh!