Bugfixes for the methods related to maps

docs/changelog.rst

@@ -34,6 +34,7 @@ Changed
 - Change default spacing in setup_channels from ``None`` to ``np.inf``. (PR #133)
 - Added ``maxdist`` variable to setup_rivers and setup_channels. (PR #153)
 - Renamed ``manhole_defaults_fn`` to ``manholes_defaults_fn`` in ``setup_manholes`` for consistency. (PR #187)
+- No data values in maps from ``setup_maps_from_rasterdataset`` are now handled as -999.0. (PR #161)
 
 Fixed
 -----

hydromt_delft3dfm/dflowfm.py

@@ -2389,10 +2389,10 @@
             'lanczos', 'average', 'mode', 'gauss', 'max', 'min', 'med', 'q1', 'q3',
             'sum', 'rms']
         interpolation_method : str, optional
-            Interpolation method for DFlow-FM. By default triangulation. Except for
-            waterlevel and waterdepth then the default is mean.
-            When methods other than 'triangulation', the relative search cell size will
-            be estimated based on resolution of the raster.
+            Interpolation method for DFlow-FM. By default mean for waterlevel and
+            waterdepth, and triangulation for all other variables. When methods other
+            than 'triangulation' are used, the relative search cell size will be
+            estimated based on resolution of the raster.
             Available methods: ['triangulation', 'mean', 'nearestNb', 'max', 'min',
             'invDist', 'minAbs', 'median']
         locationtype : str, optional
@@ -2401,15 +2401,15 @@
             Variable name, only in case data is of type DataArray or if a Dataset is
             added as is (split_dataset=False).
         split_dataset: bool, optional
-            If data is a xarray.Dataset, either add it as is to maps or split it into
-            several xarray.DataArrays.
+            If data is a xarray.Dataset, either add it as a Dataset to maps or split it
+            into a xarray.DataArrays per variable.
             Default to True.
         """
         # check for name when split_dataset is False
         if split_dataset is False and name is None:
             self.logger.error("name must be specified when split_dataset = False")
 
-        # Call super method
+        # Call super method from HydroMT Core
         variables = super().setup_maps_from_rasterdataset(
             raster_fn=raster_fn,
             variables=variables,
@@ -2419,6 +2419,11 @@
             split_dataset=split_dataset,
         )
 
+        for var in variables:
+            da = self.maps[var]
+            da.where(da == da.raster.nodata, -999.0)
+            self.set_maps(da, var)
+
         allowed_methods = [
             "triangulation",
             "mean",
@@ -2965,17 +2970,21 @@
         self._assert_read_mode()
         # Read initial fields
         inifield_model = self.dfmmodel.geometry.inifieldfile
-        # seperate 1d and 2d
-        # inifield_model_1d = [
-        #     i for i in inifield_model.initial if "1d" in i.locationtype
-        # ] # not supported yet
-        inifield_model_2dinitial = [
-            i for i in inifield_model.initial if "2d" in i.locationtype
-        ]
-        inifield_model_2dparameter = [
-            i for i in inifield_model.parameter if "2d" in i.locationtype
-        ]
-        inifield_model_2d = inifield_model_2dinitial + inifield_model_2dparameter
+        if inifield_model:
+            # seperate 1d and 2d
+            # inifield_model_1d = [
+            #     i for i in inifield_model.initial if "1d" in i.locationtype
+            # ] # not supported yet
+            inifield_model_2dinitial = [
+                i for i in inifield_model.initial if "2d" in i.locationtype
+            ]
+            inifield_model_2dparameter = [
+                i for i in inifield_model.parameter if "2d" in i.locationtype
+            ]
+            inifield_model_2d = inifield_model_2dinitial + inifield_model_2dparameter
+        else:
+            inifield_model_2d = []
+
         if any(inifield_model_2d):
             # Loop over initial / parameter to read the geotif
             inilist = inifield_model_2d
@@ -3043,6 +3052,7 @@
             if da.raster.nodata is None or np.isnan(da.raster.nodata):
                 da.raster.set_nodata(-999)
             da.raster.to_raster(_fn)
+            self.logger.info(f"Writing file {mapsroot}/{name}.tif")
             # Prepare dict
             if interp_method == "triangulation":
                 inidict = {
@@ -3083,20 +3093,26 @@
                     # update config if infiltration
                     if name == "infiltcap":
                         self.set_config("grw.infiltrationmodel", 2)
-
+                else:
+                    self.logger.error(
+                        f"Could not write map to model: {name} not recognized"
+                    )
             elif isinstance(ds, xr.Dataset):
                 for v in ds.data_vars:
                     if v in self._MAPS:
                         _prepare_inifields(self._MAPS[v], ds[v])
                         # update config if frcition
-                        if "frictype" in self._MAPS[name]:
+                        if self._MAPS[v] == "frictype":
                             self.set_config(
                                 "physics.uniffricttype", self._MAPS[name]["frictype"]
                             )
                         # update config if infiltration
-                        if name == "infiltcap":
+                        if v == "infiltcap":
                             self.set_config("grw.infiltrationmodel", 2)
-
+                    else:
+                        self.logger.error(
+                            f"Could not write map to model: {v} not found in map {name}"
+                        )
         # Assign initial fields to model and write
         inifield_model = IniFieldModel(initial=inilist, parameter=paramlist)
         # Bug: when initialising IniFieldModel hydrolib-core does not parse correclty

tests/test_dflowfm.py

@@ -78,3 +78,20 @@ def test_write_structures(tmpdir):
 
     # indirectly call hidden write_structures() method
     model.write_geoms(write_mesh_gdf=False)
+
+
+def test_setup_maps_from_rasterdataset(tmpdir):
+    model = DFlowFMModel(root=join(EXAMPLEDIR, "dflowfm_local"), mode="r")
+    model.read()
+    model.set_root(tmpdir, mode="w")
+    raster_fn = ""
+    variables = []
+    #TODO expand this test
+    #TODO assert if there are no NaN values, should be -999
+    # model.setup_maps_from_rasterdataset(raster_fn, variables)
+
+
+def test_read_maps():
+    model = DFlowFMModel(root=join(EXAMPLEDIR, "dflowfm_local"), mode="r")
+    #TODO assert if initialfields are read correctly
+    #TODO check if NaN values are read correctly