Merge branch 'forecast-class' into implement_mean_min_max

Author: luseverin

climada/engine/impact_forecast.py

@@ -25,6 +25,7 @@
 import scipy.sparse as sparse
 
 from ..util import log_level
+from ..util.checker import size
 from ..util.forecast import Forecast
 from .impact import Impact
 
@@ -52,8 +53,8 @@ def __init__(
         impact_kwargs
             Keyword-arguments passed to ~:py:class`climada.engine.impact.Impact`.
         """
-        # TODO: Maybe assert array lengths?
         super().__init__(lead_time=lead_time, member=member, **impact_kwargs)
+        self._check_sizes()
 
     @classmethod
     def from_impact(
@@ -92,13 +93,16 @@ def from_impact(
 
     @property
     def at_event(self):
-        LOGGER.warning("at_event for forecasts is not yet implemented.")
+        """Get the total impact for each member/lead_time combination."""
+        LOGGER.warning(
+            "at_event gives the total impact for one specific combination of member and "
+            "lead_time."
+        )
         return self._at_event
 
     @at_event.setter
     def at_event(self, value):
-        """Set the total exposure value close to a hazard"""
-        LOGGER.warning("at_event for forecasts is not yet implemented.")
+        """Set the total impact for each member/lead_time combination."""
         self._at_event = value
 
     def local_exceedance_impact(
@@ -111,9 +115,14 @@ def local_exceedance_impact(
         bin_decimals=None,
     ):
         """Compution of local exceedance impact for given return periods is not
-        implemented for ImpactForecast. See climada.engine.impact.Impact for details.
-        Returns
-        -------
+        implemented for ImpactForecast.
+
+        See Also
+        --------
+        See :py:meth:`~climada.engine.impact.Impact.local_exceedance_impact`
+
+        Raises
+        ------
         NotImplementedError
         """
 
@@ -132,8 +141,13 @@ def local_return_period(
         bin_decimals=None,
     ):
         """Compution of local return period for given impact thresholds is not
-        implemented for ImpactForecast. See climada.engine.impact.Impact for details.
-        Returns
+        implemented for ImpactForecast.
+
+        See Also
+        --------
+        See :py:meth:`~climada.engine.impact.Impact.local_return_period`
+
+        Raises
         -------
         NotImplementedError
         """
@@ -145,15 +159,33 @@ def local_return_period(
 
     def calc_freq_curve(self, return_per=None):
         """Computation of the impact exceedance frequency curve is not
-        implemented for ImpactForecast. See climada.engine.impact.Impact for details.
-        Returns
-        -------
+        implemented for ImpactForecast.
+
+        See Also
+        --------
+        See :py:meth:`~climada.engine.impact.Impact.calc_freq_curve`
+
+        Raises
+        ------
         NotImplementedError
         """
 
         LOGGER.error("calc_freq_curve is not defined for ImpactForecast")
         raise NotImplementedError("calc_freq_curve is not defined for ImpactForecast")
 
+    def _check_sizes(self):
+        """Check sizes of forecast data vs. impact data.
+
+        Raises
+        ------
+        ValueError
+            If the sizes of the forecast data do not match the
+            :py:attr:`~climada.engine.impact.Impact.event_id`
+        """
+        num_entries = len(self.event_id)
+        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, reduce_method: str):
         """
         Reduce the attributes of an ImpactForecast to a single value.

climada/engine/test/test_impact_forecast.py

@@ -41,13 +41,15 @@ def impact(impact_kwargs):
 
 
 @pytest.fixture
-def lead_time():
-    return pd.timedelta_range(start="1 day", periods=6).to_numpy()
+def lead_time(impact_kwargs):
+    return pd.timedelta_range(
+        start="1 day", periods=len(impact_kwargs["event_id"])
+    ).to_numpy()
 
 
 @pytest.fixture
-def member():
-    return np.arange(6)
+def member(impact_kwargs):
+    return np.arange(len(impact_kwargs["event_id"]))
 
 
 @pytest.fixture
@@ -76,6 +78,12 @@ def test_impact_forecast_init(self, impact_kwargs, lead_time, member):
         npt.assert_array_equal(forecast1.member, member)
         self.assert_impact_kwargs(forecast1, **impact_kwargs)
 
+    def test_impact_forecast_init_error(self, impact, impact_kwargs, lead_time, member):
+        with pytest.raises(ValueError, match="Forecast.lead_time"):
+            ImpactForecast(lead_time=lead_time[:-2], member=member, **impact_kwargs)
+        with pytest.raises(ValueError, match="Forecast.member"):
+            ImpactForecast.from_impact(impact, lead_time=lead_time, member=member[1:])
+
     def test_impact_forecast_from_impact(
         self, impact_forecast, impact_kwargs, lead_time, member
     ):
@@ -84,16 +92,58 @@ def test_impact_forecast_from_impact(
         self.assert_impact_kwargs(impact_forecast, **impact_kwargs)
 
 
-def test_impact_forecast_select(impact_forecast, lead_time, member, impact_kwargs):
+@pytest.mark.parametrize(
+    "var, var_select",
+    [("event_id", "event_ids"), ("event_name", "event_names"), ("date", "dates")],
+)
+def test_impact_forecast_select_events(
+    impact_forecast, lead_time, member, impact_kwargs, var, var_select
+):
     """Check if Impact.select works on the derived class"""
-    event_ids = impact_kwargs["event_id"][np.array([2, 0])]
-    impact_fc = impact_forecast.select(event_ids=event_ids)
+    select_mask = np.array([2, 1])
+    ordered_select_mask = np.array([1, 2])
+    if var == "date":
+        # Date needs to be a valid delta
+        select_mask = np.array([1, 2])
+        ordered_select_mask = np.array([1, 2])
+
+    var_value = np.array(impact_kwargs[var])[select_mask]
+    # event_name is a list, convert to numpy array for indexing
+    impact_fc = impact_forecast.select(**{var_select: var_value})
     # NOTE: Events keep their original order
     npt.assert_array_equal(
-        impact_fc.event_id, impact_forecast.event_id[np.array([0, 2])]
+        impact_fc.event_id,
+        impact_forecast.event_id[ordered_select_mask],
+    )
+    npt.assert_array_equal(
+        impact_fc.event_name,
+        np.array(impact_forecast.event_name)[ordered_select_mask],
+    )
+    npt.assert_array_equal(impact_fc.date, impact_forecast.date[ordered_select_mask])
+    npt.assert_array_equal(
+        impact_fc.frequency, impact_forecast.frequency[ordered_select_mask]
+    )
+    npt.assert_array_equal(impact_fc.member, member[ordered_select_mask])
+    npt.assert_array_equal(impact_fc.lead_time, lead_time[ordered_select_mask])
+    npt.assert_array_equal(
+        impact_fc.imp_mat.todense(),
+        impact_forecast.imp_mat.todense()[ordered_select_mask],
+    )
+
+
+def test_impact_forecast_select_exposure(
+    impact_forecast, lead_time, member, impact_kwargs
+):
+    """Check if Impact.select works on the derived class"""
+    exp_col = 0
+    select_mask = np.array([exp_col])
+    coord_exp = impact_kwargs["coord_exp"][select_mask]
+    impact_fc = impact_forecast.select(coord_exp=coord_exp)
+    npt.assert_array_equal(impact_fc.member, member)
+    npt.assert_array_equal(impact_fc.lead_time, lead_time)
+    npt.assert_array_equal(
+        impact_fc.imp_mat.todense(), impact_forecast.imp_mat.todense()[:, exp_col]
     )
-    npt.assert_array_equal(impact_fc.member, member[np.array([0, 2])])
-    npt.assert_array_equal(impact_fc.lead_time, lead_time[np.array([0, 2])])
 
 
 @pytest.mark.skip("Concat from base class does not work")

climada/hazard/forecast.py

@@ -23,8 +23,9 @@
 
 import numpy as np
 
-from climada.hazard.base import Hazard
-from climada.util.forecast import Forecast
+from ..util.checker import size
+from ..util.forecast import Forecast
+from .base import Hazard
 
 LOGGER = logging.getLogger(__name__)
 
@@ -52,6 +53,7 @@ def __init__(
             py:meth`~climada.hazard.base.Hazard.__init__` for details.
         """
         super().__init__(lead_time=lead_time, member=member, **hazard_kwargs)
+        self._check_sizes()
 
     @classmethod
     def from_hazard(cls, hazard: Hazard, lead_time: np.ndarray, member: np.ndarray):
@@ -89,3 +91,16 @@ def from_hazard(cls, hazard: Hazard, lead_time: np.ndarray, member: np.ndarray):
             intensity=hazard.intensity,
             fraction=hazard.fraction,
         )
+
+    def _check_sizes(self):
+        """Check sizes of forecast data vs. hazard data.
+
+        Raises
+        ------
+        ValueError
+            If the sizes of the forecast data do not match the
+            :py:attr:`~climada.hazard.base.Hazard.event_id`
+        """
+        num_entries = len(self.event_id)
+        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")

climada/hazard/io.py

@@ -917,6 +917,10 @@ def write_hdf5(self, file_name, todense=False):
                     # Centroids have their own write_hdf5 method,
                     # which is invoked at the end of this method (s.b.)
                     continue
+                elif var_name == "lead_time":
+                    hf_data.create_dataset(
+                        var_name, data=var_val.astype("timedelta64[ns]").astype("int64")
+                    )
                 elif isinstance(var_val, sparse.csr_matrix):
                     if todense:
                         hf_data.create_dataset(var_name, data=var_val.toarray())
@@ -987,7 +991,11 @@ def from_hdf5(cls, file_name):
                     continue
                 if var_name == "centroids":
                     continue
-                if isinstance(var_val, np.ndarray) and var_val.ndim == 1:
+                if var_name == "lead_time":
+                    hazard_kwargs[var_name] = np.array(hf_data.get(var_name)).astype(
+                        "timedelta64[ns]"
+                    )
+                elif isinstance(var_val, np.ndarray) and var_val.ndim == 1:
                     hazard_kwargs[var_name] = np.array(hf_data.get(var_name))
                 elif isinstance(var_val, sparse.csr_matrix):
                     hf_csr = hf_data.get(var_name)

climada/hazard/test/test_forecast.py

@@ -41,13 +41,13 @@ def hazard(haz_kwargs):
 
 
 @pytest.fixture
-def lead_time():
-    return pd.timedelta_range("1h", periods=6).to_numpy()
+def lead_time(haz_kwargs):
+    return pd.timedelta_range("1h", periods=len(haz_kwargs["event_id"])).to_numpy()
 
 
 @pytest.fixture
-def member():
-    return np.arange(6)
+def member(haz_kwargs):
+    return np.arange(len(haz_kwargs["event_id"]))
 
 
 @pytest.fixture
@@ -78,6 +78,13 @@ def test_init_hazard_forecast(haz_fc, member, lead_time, haz_kwargs):
     assert_hazard_kwargs(haz_fc, **haz_kwargs)
 
 
+def test_init_hazard_forecast_error(hazard, member, lead_time, haz_kwargs):
+    with pytest.raises(ValueError, match="Forecast.lead_time"):
+        HazardForecast(lead_time=lead_time[:-2], member=member, **haz_kwargs)
+    with pytest.raises(ValueError, match="Forecast.member"):
+        HazardForecast.from_hazard(hazard, lead_time=lead_time, member=member[1:])
+
+
 def test_from_hazard(lead_time, member, hazard, haz_kwargs):
     haz_fc_from_haz = HazardForecast.from_hazard(
         hazard, lead_time=lead_time, member=member
@@ -100,10 +107,60 @@ def test_hazard_forecast_concat(haz_fc, lead_time, member):
     npt.assert_array_equal(haz_fc_concat.member, np.concatenate([member, member]))
 
 
-def test_hazard_forecast_select(haz_fc, lead_time, member):
+@pytest.mark.parametrize(
+    "var, var_select",
+    [("event_id", "event_id"), ("event_name", "event_names"), ("date", "date")],
+)
+def test_hazard_forecast_select(haz_fc, lead_time, member, haz_kwargs, var, var_select):
     """Check if Hazard.select works on the derived class"""
-    haz_fc_select = haz_fc.select(event_id=[4, 1])
-    # NOTE: Events keep their original order
-    npt.assert_array_equal(haz_fc_select.event_id, haz_fc.event_id[np.array([3, 0])])
-    npt.assert_array_equal(haz_fc_select.member, member[np.array([3, 0])])
-    npt.assert_array_equal(haz_fc_select.lead_time, lead_time[np.array([3, 0])])
+
+    select_mask = np.array([3, 2])
+    ordered_select_mask = np.array([3, 2])
+    if var == "date":
+        # Date needs to be a valid delta
+        select_mask = np.array([2, 3])
+        ordered_select_mask = np.array([2, 3])
+
+    var_value = np.array(haz_kwargs[var])[select_mask]
+    # event_name is a list, convert to numpy array for indexing
+    haz_fc_sel = haz_fc.select(**{var_select: var_value})
+    # Note: order is preserved
+    npt.assert_array_equal(
+        haz_fc_sel.event_id,
+        haz_fc.event_id[ordered_select_mask],
+    )
+    npt.assert_array_equal(
+        haz_fc_sel.event_name,
+        np.array(haz_fc.event_name)[ordered_select_mask],
+    )
+    npt.assert_array_equal(haz_fc_sel.date, haz_fc.date[ordered_select_mask])
+    npt.assert_array_equal(haz_fc_sel.frequency, haz_fc.frequency[ordered_select_mask])
+    npt.assert_array_equal(haz_fc_sel.member, member[ordered_select_mask])
+    npt.assert_array_equal(haz_fc_sel.lead_time, lead_time[ordered_select_mask])
+    npt.assert_array_equal(
+        haz_fc_sel.intensity.todense(),
+        haz_fc.intensity.todense()[ordered_select_mask],
+    )
+    npt.assert_array_equal(
+        haz_fc_sel.fraction.todense(),
+        haz_fc.fraction.todense()[ordered_select_mask],
+    )
+
+    assert haz_fc_sel.centroids == haz_fc.centroids
+
+
+def test_write_read_hazard_forecast(haz_fc, tmp_path):
+
+    file_name = tmp_path / "test_hazard_forecast.h5"
+
+    haz_fc.write_hdf5(file_name)
+    haz_fc_read = HazardForecast.from_hdf5(file_name)
+
+    assert haz_fc_read.lead_time.dtype.kind == np.dtype("timedelta64").kind
+
+    for key in haz_fc.__dict__.keys():
+        if key in ["intensity", "fraction"]:
+            (haz_fc.__dict__[key] != haz_fc_read.__dict__[key]).nnz == 0
+        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])