diff --git a/climada/engine/impact_forecast.py b/climada/engine/impact_forecast.py index 9e8d2af0ba..154f890512 100644 --- a/climada/engine/impact_forecast.py +++ b/climada/engine/impact_forecast.py @@ -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") diff --git a/climada/engine/test/test_impact_forecast.py b/climada/engine/test/test_impact_forecast.py index 257fbbe2ee..e4442748c2 100644 --- a/climada/engine/test/test_impact_forecast.py +++ b/climada/engine/test/test_impact_forecast.py @@ -268,3 +268,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])) diff --git a/climada/hazard/forecast.py b/climada/hazard/forecast.py index abab0d0787..cf6980f5ff 100644 --- a/climada/hazard/forecast.py +++ b/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") diff --git a/climada/hazard/test/test_forecast.py b/climada/hazard/test/test_forecast.py index 99bc1d6d73..26f26de4b1 100644 --- a/climada/hazard/test/test_forecast.py +++ b/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]))