#152 implement setup_2dboundary, update setup_1dboundary

supports:
spatial uniform constant boundary

does not support yet:
spatial uniform timeseries boundary
spatial varying timeseries boundary

hydrolib/hydromt_delft3dfm/dflowfm.py

@@ -1404,7 +1404,7 @@ def setup_1dboundary(
         tstart and tstop entries.
 
         Adds/Updates model layers:
-            * **{boundary_type}bnd_{branch_type}** forcing: 1D boundaries DataArray
+            * **boundary1d_{boundary_type}bnd_{branch_type}** forcing: 1D boundaries DataArray
 
         Parameters
         ----------
@@ -1498,7 +1498,7 @@ def setup_1dboundary(
         )
 
         # 4. set boundaries
-        self.set_forcing(da_out, name=f"{da_out.name}_{branch_type}")
+        self.set_forcing(da_out, name=f"boundary1d_{da_out.name}_{branch_type}") #FIXME: this format cannot be read back due to lack of branch type info from model files
 
     def setup_mesh2d(
         self,
@@ -1521,9 +1521,10 @@ def setup_mesh2d(
         can already be read.
         (2) If no geometry, bbox or existing grid is specified for this setup function, then the self.region is used as
          mesh extent to generate the unstructured mesh.
-        (3) Validation checks have been added to check if the mesh extent is within model region.
-        (4) Only existing meshed with only 2D grid can be read.
-        (5) 1D2D network files are not supported as mesh2d_fn.
+        (3) At mesh border, cells that intersect with geometry border will be kept.
+        (4) Validation checks have been added to check if the mesh extent is within model region.
+        (5) Only existing meshed with only 2D grid can be read.
+        (6) 1D2D network files are not supported as mesh2d_fn.
 
         Adds/Updates model layers:
 
@@ -1781,6 +1782,132 @@ def setup_maps_from_raster_reclass(
                     )
                     self._MAPS[var]["averagingrelsize"] = relsize
 
+    def setup_2dboundary(
+        self,
+        boundaries_fn:str = None,
+        boundaries_timeseries_fn: str = None,
+        boundaries_geodataset_fn: str = None,
+        boundary_value: float = 0.0,
+        boundary_type: str = "waterlevel",
+        boundary_unit: str = "m",
+    ):
+        """
+        Prepares the 2D boundaries from line geometry.
+
+        The values can either be a spatially-uniform constant using ``boundaries_fn`` and ``boundary_value`` (default),
+        or spatially-uniform timeseries read from ``boundaries_timeseries_fn`` (TODO: Not Implemented, check the example for meteo)
+        or spatially-varying timeseries read from ``boundaries_geodataset_fn`` (TODO: Not Implemented, check if hydromt support line geometry time series)
+        The ``boundary_type`` can either be "waterlevel" or "discharge".
+
+        If ``boundaries_timeseries_fn``  or ``boundaries_geodataset_fn`` has missing values,
+        the constant ``boundary_value`` will be used.
+
+        The dataset/timeseries are clipped to the model region (see below note), and model time based on the model config
+        tstart and tstop entries.
+
+        Note that:
+        (1) Only line geometry that are contained within the region of 0-3 grid cell sizes outside of the mesh2d will be kept.
+        (2) Because of the above, this function must be called before the mesh refinement.
+        (3) when using constant boundary, the output forcing will be written to time series with constant values.
+
+        Adds/Updates model layers:
+            * **boundary2d_{boundary_name}** forcing: 2D boundaries DataArray
+
+        Parameters
+        ----------
+        boundaries_fn: str Path
+            Path or data source name for line geometry file.
+            * Required variables if a combined time series data csv file: ['index'] with type int
+        boundaries_timeseries_fn: str, Path
+            Path to tabulated timeseries csv file with time index in first column
+            and location index with type int in the first row,
+            see :py:meth:`hydromt.open_timeseries_from_table`, for details.
+            NOTE: Require equidistant time series
+        boundaries_geodataset_fn : str, Path
+            Path or data source name for geospatial point timeseries file.
+            This can either be a netcdf file with geospatial coordinates
+            or a combined point location file with a timeseries data csv file
+            which can be setup through the data_catalog yml file, see hydromt User Manual (https://deltares.github.io/hydromt/latest/_examples/prep_data_catalog.html#GeoDataset-from-vector-files)
+            * Required variables if netcdf: ['discharge', 'waterlevel'] depending on ``boundary_type``
+            * Required coordinates if netcdf: ['time', 'index', 'y', 'x']
+            * Required variables if a combined point location file: ['index'] with type int
+            * Required index types if a time series data csv file: int
+            NOTE: Require equidistant time series
+        boundary_value : float, optional
+            Constant value to use for all boundaries, and to
+            fill in missing data. By default 0.0 m.
+        boundary_type : str, optional
+            Type of boundary tu use. One of ["waterlevel", "discharge"].
+            By default "waterlevel".
+        boundary_unit : str, optional.
+            Unit corresponding to [boundary_type].
+            If ``boundary_type`` = "waterlevel"
+                Allowed unit is [m]
+            if ''boundary_type`` = "discharge":
+               Allowed unit is [m3/s]
+            By default m.
+
+        Raises:
+        -------
+        NotImplementedError:
+            The use of boundaries_timeseries_fn and boundaries_geodataset_fn is not yet implemented.
+
+        """
+        self.logger.info(f"Preparing 2D boundaries.")
+
+        if boundary_type == "waterlevel": assert boundary_unit in ["m"]
+        if boundary_type == "discharge": assert boundary_unit in ["m3/s"]
+
+        _mesh_region = self._mesh.ugrid.to_geodataframe().unary_union
+        _boundary_region = _mesh_region.buffer(3*self.res).difference(_mesh_region) # region where 2d boundary is allowed
+        _boundary_region = gpd.GeoDataFrame({'geometry':[_boundary_region]}, crs = self.crs)
+
+        refdate, tstart, tstop = self.get_model_time()  # time slice
+
+        # 1. read constant boundaries
+        if boundaries_fn is not None:
+            gdf_bnd = self.data_catalog.get_geodataframe(
+                boundaries_fn,
+                geom=_boundary_region,
+                crs = self.crs,
+                predicate='contains',
+            )
+            # preprocess
+            gdf_bnd = gdf_bnd.explode()
+            # set index
+            if "boundary_id" not in gdf_bnd:
+                gdf_bnd["boundary_id"] = [f"2dboundary_{i}" for i in range(len(gdf_bnd))]
+        else:
+            gdf_bnd = None
+        # 2. read timeseries boundaries
+        if boundaries_timeseries_fn is not None:
+            raise NotImplementedError("Does not support reading timeseries boundaries yet.")
+        else:
+            df_bnd = pd.DataFrame({"time":
+                                       pd.date_range(start=pd.to_datetime(tstart), end = pd.to_datetime(tstop), freq='D'),
+                                   "global":
+                                       np.nan})
+        # 3. read spatially varying timeseries boundaries
+        if boundaries_geodataset_fn is not None:
+            raise NotImplementedError("Does not support reading geodataset boundaries yet.")
+        else:
+            da_bnd = None
+
+        # 4. Derive DataArray with boundary values at boundary locations in boundaries_branch_type
+        da_out = workflows.compute_2dboundary_values(
+            boundaries=gdf_bnd,
+            da_bnd=da_bnd,
+            df_bnd=df_bnd,
+            boundary_value=boundary_value,
+            boundary_type=boundary_type,
+            boundary_unit=boundary_unit,
+            logger=self.logger,
+        )
+
+        # 5. set boundaries
+        self.set_forcing(da_out, name=f"boundary2d_{da_out.name}")
+
+
     # ## I/O
     # TODO: remove after hydromt 0.6.1 release
     @property
@@ -2009,7 +2136,7 @@ def _prepare_inifields(da_dict, da):
         # save filepath in the config
         self.set_config("geometry.inifieldfile", inifield_model_filename)
 
-    def read_geoms(self) -> None:
+    def read_geoms(self) -> None: #FIXME: gives an error when only 2D model.
         """
         Read model geometries files at <root>/<geoms> and add to geoms property.
 
@@ -2070,26 +2197,37 @@ def write_geoms(self) -> None:
             )
             self.set_config("geometry.storagenodefile", storage_fn)
 
-    def read_forcing(self) -> None:
+    def read_forcing(self) -> None: #FIXME reading of forcing should include boundary, lateral and meteo
         """Read forcing at <root/?/> and parse to dict of xr.DataArray"""
         self._assert_read_mode
         # Read external forcing
         ext_model = self.dfmmodel.external_forcing.extforcefilenew
         if ext_model is not None:
+            # read boundary blocks #FIXME: there might be better options
             df_ext = pd.DataFrame([f.__dict__ for f in ext_model.boundary])
-            # Forcing dataarrays to prepare for each quantity
-            forcing_names = np.unique(df_ext.quantity).tolist()
-            # Loop over forcing names to build data arrays
-            for name in forcing_names:
-                # Get the dataframe corresponding to the current variable
-                df = df_ext[df_ext.quantity == name]
-                # Get the corresponding nodes gdf
-                node_geoms = self.network1d_nodes[
-                    np.isin(self.network1d_nodes["nodeId"], df.nodeid.values)
-                ]
-                da_out = utils.read_1dboundary(df, quantity=name, nodes=node_geoms)
-                # Add to forcing
-                self.set_forcing(da_out)
+            # 1d boundary
+            df_ext_1d = df_ext.loc[~df_ext.nodeid.isna(),:]
+            if len(df_ext_1d) > 0:
+                # Forcing data arrays to prepare for each quantity
+                forcing_names = np.unique(df_ext_1d.quantity).tolist()
+                # Loop over forcing names to build data arrays
+                for name in forcing_names:
+                    # Get the dataframe corresponding to the current variable
+                    df = df_ext_1d[df_ext_1d.quantity == name]
+                    # Get the corresponding nodes gdf
+                    node_geoms = self.network1d_nodes[
+                        np.isin(self.network1d_nodes["nodeId"], df.nodeid.values)
+                    ]
+                    da_out = utils.read_1dboundary(df, quantity=name, nodes=node_geoms)
+                    # Add to forcing
+                    self.set_forcing(da_out)
+            # 2d boundary
+            df_ext_2d = df_ext.loc[df_ext.nodeid.isna(), :]
+            if len(df_ext_2d) > 0:
+                for _,df in df_ext_2d.iterrows():
+                    da_out = utils.read_2dboundary(df, workdir = self.dfmmodel.filepath.parent)
+                    # Add to forcing
+                    self.set_forcing(da_out)
 
     def write_forcing(self) -> None:
         """write forcing into hydrolib-core ext and forcing models"""
@@ -2099,7 +2237,12 @@ def write_forcing(self) -> None:
             self._assert_write_mode
             self.logger.info("Writting forcing files.")
             savedir = dirname(join(self.root, self._config_fn))
-            _, ext_fn = utils.write_1dboundary(self.forcing, savedir)
+            # create new external forcing file
+            ext_fn = "bnd.ext"
+            Path(join(savedir, ext_fn)).unlink(missing_ok=True)
+            # populate external forcing file
+            utils.write_1dboundary(self.forcing, savedir, ext_fn=ext_fn)
+            utils.write_2dboundary(self.forcing, savedir, ext_fn=ext_fn)
             self.set_config("external_forcing.extforcefilenew", ext_fn)
 
     def read_mesh(self):