Add forecast classes (#1115)

* Add HazardForecast and ImpactForecast classes.
* Add tests.
* Add new forecast tutorial.

---------

Co-authored-by: Evelyn Mühlhofer <42063186+Evelyn-M@users.noreply.github.com>
Co-authored-by: Luca Severino <91593121+luseverin@users.noreply.github.com>
Co-authored-by: Valentin Gebhart <60438839+ValentinGebhart@users.noreply.github.com>
Co-authored-by: Chahan M. Kropf <chahan.kropf@usys.ethz.ch>
Co-authored-by: Chahan Kropf <chahan.kropf@posteo.com>
Co-authored-by: elianekobler <138141646+elianekobler@users.noreply.github.com>

Author: 7 people

CHANGELOG.md

@@ -12,12 +12,14 @@ Code freeze date: YYYY-MM-DD
 
 ### Added
 
+- Support for ensemble forecast data as hazard via `HazardForecast` and `ImpactForecast`, see `doc.user_guide.climada_engine_ImpactForecast.ipynb` [#1115](https://github.com/CLIMADA-project/climada_python/pull/)
 - Better type hints and overloads signatures for ImpactFuncSet [#1250](https://github.com/CLIMADA-project/climada_python/pull/1250)
 - Adds `MeasureConfig` and related dataclasses for new `Measure` object retrocompatibility and (de)serialization capabilities [#1276](https://github.com/CLIMADA-project/climada_python/pull/1276)
 - Add inter- and extrapolation options to `ImpactFreqCurve` with method `interpolate` [#1252](https://github.com/CLIMADA-project/climada_python/pull/1252)
 - Add `trajectory` module to ease the handling of multiple risk assessment at different point in time, and enable interpolation in between. [#1197](https://github.com/CLIMADA-project/climada_python/pull/1197), [#1198](https://github.com/CLIMADA-project/climada_python/pull/1198), [#1199](https://github.com/CLIMADA-project/climada_python/pull/1199), [#1200](https://github.com/CLIMADA-project/climada_python/pull/1197), [#1201](https://github.com/CLIMADA-project/climada_python/pull/1197)
 
 ### Changed
+
 - Updated Impact Calculation Tutorial (`doc.climada_engine_Impact.ipynb`) [#1095](https://github.com/CLIMADA-project/climada_python/pull/1095).
 - Makes current `measure` module a legacy module, moving it to `_legacy_measure`, to retain compatibility with `CostBenefit` class and various tests. [#1274](https://github.com/CLIMADA-project/climada_python/pull/1274)
 

climada/engine/__init__.py

@@ -22,3 +22,4 @@
 from .cost_benefit import *
 from .impact import *
 from .impact_calc import *
+from .impact_forecast import ImpactForecast

climada/engine/impact.py

@@ -188,8 +188,7 @@ def __init__(
             )
         if len(self.coord_exp) != len(self.eai_exp):
             raise AttributeError(
-                "Number of exposures points is different from"
-                "number of eai_exp values"
+                "Number of exposures points is different fromnumber of eai_exp values"
             )
         if imp_mat is not None:
             self.imp_mat = imp_mat
@@ -1226,8 +1225,9 @@ def plot_rp_imp(
 
         impacts_stats_vals = impacts_stats.values[:, 1:].T.astype(float)
         if not log10_scale:
-            min_impact, max_impact = np.nanmin(impacts_stats_vals), np.nanmax(
-                impacts_stats_vals
+            min_impact, max_impact = (
+                np.nanmin(impacts_stats_vals),
+                np.nanmax(impacts_stats_vals),
             )
             kwargs.update(
                 {
@@ -1419,6 +1419,8 @@ def write_attribute(group, name, value):
 
         def write_dataset(group, name, value):
             """Write a dataset"""
+            if name == "lead_time":
+                value = value.astype("timedelta64[ns]").astype("int64")
             group.create_dataset(name, data=value, dtype=_str_type_helper(value))
 
         def write_dict(group, name, value):
@@ -1459,7 +1461,6 @@ def write_csr(group, name, value):
 
         # Open file in write mode
         with h5py.File(file_path, "w") as file:
-
             # Now write all attributes
             # NOTE: Remove leading underscore to write '_tot_value' as regular attribute
             for name, value in self.__dict__.items():
@@ -1600,13 +1601,18 @@ def from_excel(cls, file_name):
     def read_excel(self, *args, **kwargs):
         """This function is deprecated, use Impact.from_excel instead."""
         LOGGER.warning(
-            "The use of Impact.read_excel is deprecated."
-            "Use Impact.from_excel instead."
+            "The use of Impact.read_excel is deprecated.Use Impact.from_excel instead."
         )
         self.__dict__ = Impact.from_excel(*args, **kwargs).__dict__
 
     @classmethod
-    def from_hdf5(cls, file_path: Union[str, Path]):
+    def from_hdf5(
+        cls,
+        file_path: Union[str, Path],
+        *,
+        add_scalar_attrs: Iterable[str] | None = None,
+        add_array_attrs: Iterable[str] | None = None,
+    ):
         """Create an impact object from an H5 file.
 
         This assumes a specific layout of the file. If values are not found in the
@@ -1651,6 +1657,10 @@ def from_hdf5(cls, file_path: Union[str, Path]):
         ----------
         file_path : str or Path
             The file path of the file to read.
+        add_scalar_attrs : Iterable of str, optional
+            Additional scalar attributes to read from file. Defaults to None.
+        add_array_attrs : Iterable of str, optional
+            Additional array attributes to read from file. Defaults to None.
 
         Returns
         -------
@@ -1659,7 +1669,6 @@ def from_hdf5(cls, file_path: Union[str, Path]):
         """
         kwargs = dict()
         with h5py.File(file_path, "r") as file:
-
             # Impact matrix
             if "imp_mat" in file:
                 impact_matrix = file["imp_mat"]
@@ -1679,17 +1688,27 @@ def from_hdf5(cls, file_path: Union[str, Path]):
             # Scalar attributes
             scalar_attrs = set(
                 ("crs", "tot_value", "unit", "aai_agg", "frequency_unit", "haz_type")
-            ).intersection(file.attrs.keys())
+            )
+            if add_scalar_attrs is not None:
+                scalar_attrs = scalar_attrs.union(add_scalar_attrs)
+            scalar_attrs = scalar_attrs.intersection(file.attrs.keys())
             kwargs.update({attr: file.attrs[attr] for attr in scalar_attrs})
 
             # Array attributes
             # NOTE: Need [:] to copy array data. Otherwise, it would be a view that is
             #       invalidated once we close the file.
             array_attrs = set(
                 ("event_id", "date", "coord_exp", "eai_exp", "at_event", "frequency")
-            ).intersection(file.keys())
+            )
+            if add_array_attrs is not None:
+                array_attrs = array_attrs.union(add_array_attrs)
+            array_attrs = array_attrs.intersection(file.keys())
             kwargs.update({attr: file[attr][:] for attr in array_attrs})
-
+            # correct lead_time attribut to timedelta
+            if "lead_time" in kwargs:
+                kwargs["lead_time"] = np.array(file["lead_time"][:]).astype(
+                    "timedelta64[ns]"
+                )
             # Special handling for 'event_name' because it should be a list of strings
             if "event_name" in file:
                 # pylint: disable=no-member
@@ -2196,9 +2215,12 @@ def stack_attribute(attr_name: str) -> np.ndarray:
         imp_mat = sparse.vstack(imp_mats)
 
         # Concatenate other attributes
-        kwargs = {
-            attr: stack_attribute(attr) for attr in ("date", "frequency", "at_event")
-        }
+        concat_attrs = {
+            name.lstrip("_")  # Private attributes with getter/setter
+            for name, value in first_imp.__dict__.items()
+            if isinstance(value, np.ndarray)
+        }.difference(("event_id", "coord_exp", "eai_exp", "aai_agg"))
+        kwargs = {attr: stack_attribute(attr) for attr in concat_attrs}
 
         # Get remaining attributes from first impact object in list
         return cls(

climada/engine/impact_calc.py

@@ -29,6 +29,8 @@
 
 from climada import CONFIG
 from climada.engine.impact import Impact
+from climada.engine.impact_forecast import ImpactForecast
+from climada.hazard.forecast import HazardForecast
 
 LOGGER = logging.getLogger(__name__)
 
@@ -217,7 +219,7 @@ def _return_impact(self, imp_mat_gen, save_mat):
 
         Returns
         -------
-        Impact
+        Impact or ImpactForecast
             Impact Object initialize from the impact matrix
 
         See Also
@@ -230,12 +232,31 @@ def _return_impact(self, imp_mat_gen, save_mat):
             at_event, eai_exp, aai_agg = self.risk_metrics(
                 imp_mat, self.hazard.frequency
             )
+            if isinstance(self.hazard, HazardForecast):
+                eai_exp = np.full_like(eai_exp, np.nan, dtype=eai_exp.dtype)
+                aai_agg = np.nan
+                LOGGER.warning(
+                    "eai_exp and aai_agg are undefined with forecasts. "
+                    "Setting them to NaN arrays."
+                )
+
         else:
+            if isinstance(self.hazard, HazardForecast):
+                raise ValueError(
+                    "Saving impact matrix is required when using HazardForecast."
+                    "Please set save_mat=True."
+                )
             imp_mat = None
             at_event, eai_exp, aai_agg = self.stitch_risk_metrics(imp_mat_gen)
-        return Impact.from_eih(
+
+        impact = Impact.from_eih(
             self.exposures, self.hazard, at_event, eai_exp, aai_agg, imp_mat
         )
+        if isinstance(self.hazard, HazardForecast):
+            return ImpactForecast.from_impact(
+                impact, self.hazard.lead_time, self.hazard.member
+            )
+        return impact
 
     def _return_empty(self, save_mat):
         """
@@ -248,21 +269,37 @@ def _return_empty(self, save_mat):
 
         Returns
         -------
-        Impact
+        Impact or ImpactForecast
             Empty impact object with correct array sizes.
         """
         at_event = np.zeros(self.n_events)
-        eai_exp = np.zeros(self.n_exp_pnt)
-        aai_agg = 0.0
+        if isinstance(self.hazard, HazardForecast):
+            eai_exp = np.full(self.n_exp_pnt, np.nan)
+            aai_agg = np.nan
+        else:
+            eai_exp = np.zeros(self.n_exp_pnt)
+            aai_agg = 0.0
+
         if save_mat:
             imp_mat = sparse.csr_matrix(
                 (self.n_events, self.n_exp_pnt), dtype=np.float64
             )
         else:
+            if isinstance(self.hazard, HazardForecast):
+                raise ValueError(
+                    "Saving impact matrix is required when using HazardForecast. "
+                    "Please set save_mat=True."
+                )
             imp_mat = None
-        return Impact.from_eih(
+
+        impact = Impact.from_eih(
             self.exposures, self.hazard, at_event, eai_exp, aai_agg, imp_mat
         )
+        if isinstance(self.hazard, HazardForecast):
+            return ImpactForecast.from_impact(
+                impact, self.hazard.lead_time, self.hazard.member
+            )
+        return impact
 
     def minimal_exp_gdf(
         self, impf_col, assign_centroids, ignore_cover, ignore_deductible