Implement mean min max for impact forecast and hazard forecast (#1187)

Add min/mean/max methods

---------

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

Author: 3 people

climada/engine/impact_forecast.py

@@ -22,6 +22,7 @@
 import logging
 
 import numpy as np
+import scipy.sparse as sparse
 
 from ..util import log_level
 from ..util.checker import size
@@ -185,6 +186,123 @@ def _check_sizes(self):
         size(exp_len=num_entries, var=self.member, var_name="Forecast.member")
         size(exp_len=num_entries, var=self.lead_time, var_name="Forecast.lead_time")
 
+    def _reduce_attrs(self, event_name: str):
+        """
+        Reduce the attributes of an ImpactForecast to a single value.
+
+        Attributes are modified as follows:
+        - lead_time: set to NaT
+        - member: set to -1
+        - event_id: set to 0
+        - event_name: set to the name of the reduction method (default)
+        - date: set to 0
+        - frequency: set to 1
+
+        Parameters
+        ----------
+        event_name : str
+            The event name given to the reduced data.
+        """
+        reduced_attrs = {
+            "lead_time": np.array([np.timedelta64("NaT")]),
+            "member": np.array([-1]),
+            "event_id": np.array([0]),
+            "event_name": np.array([event_name]),
+            "date": np.array([0]),
+            "frequency": np.array([1]),
+        }
+
+        return reduced_attrs
+
+    def min(self):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the minimum
+        value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        ImpactForecast
+            An ImpactForecast object with the min impact matrix and at_event.
+        """
+        red_imp_mat = self.imp_mat.min(axis=0).tocsr()
+        red_at_event = np.array([red_imp_mat.sum()])
+        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("min"),
+        )
+
+    def max(self):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the maximum
+        value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        ImpactForecast
+            An ImpactForecast object with the max impact matrix and at_event.
+        """
+        red_imp_mat = self.imp_mat.max(axis=0).tocsr()
+        red_at_event = np.array([red_imp_mat.sum()])
+        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("max"),
+        )
+
+    def mean(self):
+        """
+        Reduce the impact matrix and at_event of an ImpactForecast to the mean value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        ImpactForecast
+            An ImpactForecast object with the mean impact matrix and at_event.
+        """
+        red_imp_mat = sparse.csr_matrix(self.imp_mat.mean(axis=0))
+        red_at_event = np.array([red_imp_mat.sum()])
+        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("mean"),
+        )
+
     def select(
         self,
         event_ids=None,
@@ -205,10 +323,6 @@ def select(
         lead_time : Sequence of numpy.timedelta64
             Lead times to select
 
-        Returns
-        -------
-        ImpactForecast
-
         See Also
         --------
         :py:meth:`~climada.engine.impact.Impact.select`

climada/engine/test/test_impact_forecast.py

@@ -212,7 +212,7 @@ def test_impact_forecast_concat(impact_forecast, member):
 
 
 def test_impact_forecast_blocked_methods(impact_forecast):
-    """Check if blocked methods raise NotImplementedError"""
+    """Check if ImpactForecast.exceedance_freq_curve raises NotImplementedError"""
     with pytest.raises(NotImplementedError):
         impact_forecast.local_exceedance_impact(np.array([10, 50, 100]))
 
@@ -221,3 +221,40 @@ def test_impact_forecast_blocked_methods(impact_forecast):
 
     with pytest.raises(NotImplementedError):
         impact_forecast.calc_freq_curve(np.array([10, 50, 100]))
+
+
+@pytest.fixture
+def impact_forecast_stats(impact_kwargs, lead_time, member):
+    max_index = 4
+    for key, val in impact_kwargs.items():
+        if isinstance(val, (np.ndarray, list)):
+            impact_kwargs[key] = val[:max_index]
+        elif isinstance(val, csr_matrix):
+            impact_kwargs[key] = val[:max_index, :]
+    impact_kwargs["imp_mat"] = csr_matrix([[1, 0], [0, 1], [3, 2], [2, 3]])
+    impact_kwargs["at_event"] = np.array([1, 1, 5, 5])
+    return ImpactForecast(
+        lead_time=lead_time[:max_index], member=member[:max_index], **impact_kwargs
+    )
+
+
+@pytest.mark.parametrize("attr", ["min", "mean", "max"])
+def test_impact_forecast_min_mean_max(impact_forecast_stats, attr):
+    """Check mean, min, and max methods for ImpactForecast"""
+    imp_fc_reduced = getattr(impact_forecast_stats, attr)()
+
+    # assert imp_mat
+    npt.assert_array_equal(
+        imp_fc_reduced.imp_mat.todense(),
+        getattr(impact_forecast_stats.imp_mat.todense(), attr)(axis=0),
+    )
+    at_event_expected = {"min": [0], "mean": [3], "max": [6]}
+    npt.assert_array_equal(imp_fc_reduced.at_event, at_event_expected[attr])
+
+    # check that attributes where reduced correctly
+    npt.assert_array_equal(np.isnat(imp_fc_reduced.lead_time), [True])
+    npt.assert_array_equal(imp_fc_reduced.member, [-1])
+    npt.assert_array_equal(imp_fc_reduced.event_name, [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])

climada/hazard/forecast.py

@@ -22,6 +22,7 @@
 import logging
 
 import numpy as np
+import scipy.sparse as sparse
 
 from ..util.checker import size
 from ..util.forecast import Forecast
@@ -105,6 +106,115 @@ def _check_sizes(self):
         size(exp_len=num_entries, var=self.member, var_name="Forecast.member")
         size(exp_len=num_entries, var=self.lead_time, var_name="Forecast.lead_time")
 
+    def _reduce_attrs(self, event_name: str):
+        """
+        Reduce the attributes of a HazardForecast to a single value.
+
+        Attributes are modified as follows:
+        - lead_time: set to NaT
+        - member: set to -1
+        - event_id: set to 0
+        - event_name: set to the name of the reduction method (default)
+        - date: set to 0
+        - frequency: set to 1
+
+        Parameters
+        ----------
+        event_name : str
+            The event_name given to the reduced data.
+        """
+        reduced_attrs = {
+            "lead_time": np.array([np.timedelta64("NaT")]),
+            "member": np.array([-1]),
+            "event_id": np.array([0]),
+            "event_name": np.array([event_name]),
+            "date": np.array([0]),
+            "frequency": np.array([1]),
+            "orig": np.array([True]),
+        }
+
+        return reduced_attrs
+
+    def min(self):
+        """
+        Reduce the intensity and fraction of a HazardForecast to the minimum
+        value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        HazardForecast
+            A HazardForecast object with the min intensity and fraction.
+        """
+        red_intensity = self.intensity.min(axis=0).tocsr()
+        red_fraction = self.fraction.min(axis=0).tocsr()
+        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("min"),
+        )
+
+    def max(self):
+        """
+        Reduce the intensity and fraction of a HazardForecast to the maximum
+        value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        HazardForecast
+            A HazardForecast object with the min intensity and fraction.
+        """
+        red_intensity = self.intensity.max(axis=0).tocsr()
+        red_fraction = self.fraction.max(axis=0).tocsr()
+        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("max"),
+        )
+
+    def mean(self):
+        """
+        Reduce the intensity and fraction of a HazardForecast to the mean value.
+
+        Parameters
+        ----------
+        None
+
+        Returns
+        -------
+        HazardForecast
+            A HazardForecast object with the min intensity and fraction.
+        """
+        red_intensity = sparse.csr_matrix(self.intensity.mean(axis=0))
+        red_fraction = sparse.csr_matrix(self.fraction.mean(axis=0))
+        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("mean"),
+        )
+
     @classmethod
     def concat(cls, haz_list: list):
         """Concatenate multiple HazardForecast instances and return a new object"""
@@ -138,10 +248,6 @@ def select(
         lead_time : Sequence of numpy.timedelta64
             Lead times to select
 
-        Returns
-        -------
-        HazardForecast
-
         See Also
         --------
         :py:meth:`~climada.hazard.base.Hazard.select`

climada/hazard/test/test_forecast.py

@@ -233,3 +233,28 @@ def test_write_read_hazard_forecast(haz_fc, tmp_path):
         else:
             # npt.assert_array_equal also works for comparing int, float or list
             npt.assert_array_equal(haz_fc.__dict__[key], haz_fc_read.__dict__[key])
+
+
+@pytest.mark.parametrize("attr", ["min", "mean", "max"])
+def test_hazard_forecast_mean_min_max(haz_fc, attr):
+    """Check mean, min, and max methods for ImpactForecast"""
+    haz_fcst_reduced = getattr(haz_fc, attr)()
+
+    # Assert sparse matrices
+    npt.assert_array_equal(
+        haz_fcst_reduced.intensity.todense(),
+        getattr(haz_fc.intensity.todense(), attr)(axis=0),
+    )
+    npt.assert_array_equal(
+        haz_fcst_reduced.fraction.todense(),
+        getattr(haz_fc.fraction.todense(), attr)(axis=0),
+    )
+
+    # Check that attributes where reduced correctly
+    npt.assert_array_equal(np.isnat(haz_fcst_reduced.lead_time), [True])
+    npt.assert_array_equal(haz_fcst_reduced.member, [-1])
+    npt.assert_array_equal(haz_fcst_reduced.event_name, [attr])
+    npt.assert_array_equal(haz_fcst_reduced.event_id, [0])
+    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])