diff --git a/src/itzi/providers/grass_interface.py b/src/itzi/providers/grass_interface.py index 66e0aff..32b433b 100644 --- a/src/itzi/providers/grass_interface.py +++ b/src/itzi/providers/grass_interface.py @@ -444,9 +444,25 @@ def read_raster_map(self, rast_name): with self.raster_lock: with raster.RasterRow(rast_name, mode="r") as rast: array = np.array(rast, dtype=self.dtype) + array = self._replace_cell_null_sentinel(rast.mtype, array) else: with raster.RasterRow(rast_name, mode="r") as rast: array = np.array(rast, dtype=self.dtype) + array = self._replace_cell_null_sentinel(rast.mtype, array) + return array + + @staticmethod + def _replace_cell_null_sentinel(raster_type: str, array: np.ndarray) -> np.ndarray: + """Normalize GRASS CELL nulls to NaN. + + FCELL/DCELL nulls are already exposed as NaN by pygrass. CELL nulls are + returned as the int32 null sentinel cast to the target dtype. + """ + if raster_type != "CELL" or not np.issubdtype(array.dtype, np.floating): + return array + + null_sentinel = array.dtype.type(np.iinfo(np.int32).min) + array[array == null_sentinel] = np.nan return array def write_raster_map(self, arr, rast_name, mkey, hmin): diff --git a/src/itzi/simulation.py b/src/itzi/simulation.py index 60157db..641358e 100644 --- a/src/itzi/simulation.py +++ b/src/itzi/simulation.py @@ -352,7 +352,9 @@ def update_input_arrays(self): return self # DEM is needed for WSE and rain routing direction if not self.timed_arrays["dem"].is_valid(self.sim_time): - self.set_array("dem", self.timed_arrays["dem"].get(self.sim_time)) + new_dem = self.timed_arrays["dem"].get(self.sim_time) + self._validate_input_array_data("dem", new_dem) + self.set_array("dem", new_dem) # loop through the arrays for arr_key, ta in self.timed_arrays.items(): # DEM done before @@ -381,9 +383,31 @@ def update_input_arrays(self): new_arr = ta.get(self.sim_time) # update array msgr.debug("{}: update input array <{}>".format(self.sim_time, arr_key)) + self._validate_input_array_data(arr_key, new_arr) self.set_array(arr_key, new_arr) return self + def _validate_input_array_data(self, arr_key: str, arr: np.ndarray) -> None: + """Fail or warn when a provider-backed input has no valid cells.""" + active_cells = arr[~self.raster_domain.mask] + active_cell_count = active_cells.size + finite_cell_count = int(np.count_nonzero(np.isfinite(active_cells))) + + if finite_cell_count > 0: + return + + if active_cell_count == 0: + msg = f"{self.sim_time}: active domain contains no cells for input map <{arr_key}>" + else: + msg = ( + f"{self.sim_time}: input map <{arr_key}> contains only NULL/NaN cells " + "inside the active domain" + ) + if arr_key == "dem": + msgr.fatal(msg) + else: + msgr.warning(msg) + def set_array(self, arr_id: str, arr: np.ndarray): """Set an array of the simulation domain.""" if arr_id in ["inflow", "rain"]: diff --git a/tests/core/test_hotstart_timed_inputs.py b/tests/core/test_hotstart_timed_inputs.py index f4e7b1d..5998f31 100644 --- a/tests/core/test_hotstart_timed_inputs.py +++ b/tests/core/test_hotstart_timed_inputs.py @@ -3,6 +3,7 @@ from __future__ import annotations from datetime import datetime, timedelta +import logging from typing import TYPE_CHECKING import numpy as np @@ -14,6 +15,7 @@ from itzi.const import InfiltrationModelType, TemporalType from itzi.data_containers import SimulationConfig, SurfaceFlowParameters +from itzi.itzi_error import ItziFatal from itzi.providers.memory_output import MemoryRasterOutputProvider, MemoryVectorOutputProvider from itzi.providers.xarray_input import XarrayRasterInputProvider from itzi.simulation_builder import SimulationBuilder @@ -252,3 +254,56 @@ def test_resume_with_absolute_time_xarray_inputs(domain_5by5) -> None: domain_5by5, use_relative_time=False, ) + + +def test_build_fails_when_dem_input_has_only_nan_cells(domain_5by5) -> None: + start_time = datetime(2000, 1, 1, 0, 0, 0) + end_time = start_time + timedelta(seconds=25) + dataset, _ = _make_xarray_input_dataset( + domain_5by5, + start_time, + use_relative_time=False, + ) + dataset["dem"] = xr.DataArray( + np.full(domain_5by5.domain_data.shape, np.nan, dtype=np.float32), + dims=("y", "x"), + coords={"y": dataset["y"], "x": dataset["x"]}, + ) + sim_config = _make_simulation_config( + start_time, + end_time, + temporal_type=TemporalType.ABSOLUTE, + ) + + with pytest.raises(ItziFatal, match=r"input map contains only NULL/NaN cells"): + _build_provider_simulation(sim_config, domain_5by5, dataset) + + +def test_build_warns_when_non_dem_input_has_only_nan_cells(domain_5by5, caplog) -> None: + start_time = datetime(2000, 1, 1, 0, 0, 0) + end_time = start_time + timedelta(seconds=25) + dataset, _ = _make_xarray_input_dataset( + domain_5by5, + start_time, + use_relative_time=False, + ) + dataset["rain"] = xr.DataArray( + np.full(domain_5by5.domain_data.shape, np.nan, dtype=np.float32), + dims=("y", "x"), + coords={"y": dataset["y"], "x": dataset["x"]}, + ) + sim_config = _make_simulation_config( + start_time, + end_time, + temporal_type=TemporalType.ABSOLUTE, + ) + + with caplog.at_level(logging.WARNING, logger="itzi"): + simulation = _build_provider_simulation(sim_config, domain_5by5, dataset) + + warning_messages = [record.message for record in caplog.records] + assert any( + "input map contains only NULL/NaN cells inside the active domain" in message + for message in warning_messages + ) + assert np.allclose(simulation.raster_domain.get_array("rain"), 0.0) diff --git a/tests/grass/test_itzi.py b/tests/grass/test_itzi.py index 2650b96..0381c9e 100644 --- a/tests/grass/test_itzi.py +++ b/tests/grass/test_itzi.py @@ -1,5 +1,6 @@ """ """ +from configparser import ConfigParser import os import grass.script as gscript @@ -7,6 +8,7 @@ from itzi import SimulationRunner from itzi.configreader import ConfigReader +from itzi.itzi_error import ItziFatal @pytest.mark.forked @@ -66,3 +68,45 @@ def test_region_mask(test_data_path): # Check tear down assert int(gscript.parse_command("g.region", flags="pg")["cells"]) == init_ncells assert int(gscript.parse_command("r.univar", map="z", flags="g")["null_cells"]) == init_nulls + + +@pytest.mark.forked +@pytest.mark.usefixtures("grass_5by5") +def test_fails_when_region_has_no_dem_data(test_data_temp_path): + gscript.run_command( + "g.region", res=10, s=100, n=150, w=100, e=150, save="outside_dem", flags="o" + ) + + config_dict = { + "input": { + "dem": "z@5by5", + "friction": "n@5by5", + "water_depth": "start_h@5by5", + }, + "time": { + "duration": "00:01:00", + "record_step": "00:00:30", + }, + "output": { + "prefix": "out_5by5_no_dem_overlap", + "values": "water_depth", + }, + "options": { + "hmin": "0.0001", + "dtmax": "0.3", + "cfl": "0.2", + }, + "grass": { + "region": "outside_dem", + }, + } + parser = ConfigParser() + parser.read_dict(config_dict) + config_file = os.path.join(test_data_temp_path, "outside_dem.ini") + with open(config_file, "w") as file_handle: + parser.write(file_handle) + + conf_data = ConfigReader(config_file) + + with pytest.raises(ItziFatal, match=r"input map contains only NULL/NaN cells"): + SimulationRunner(conf_data.get_sim_params(), conf_data.get_grass_params())