re-added two functions in dflowfm.py

hydromt_delft3dfm/dflowfm.py

@@ -1610,6 +1610,63 @@ def setup_1dboundary(
         # FIXME: this format cannot be read back due to lack of branch type info
         # from model files
 
+    def _read_forcing_geodataset(
+        self,
+        forcing_geodataset_fn: Union[str, Path],
+        forcing_name: str = "discharge",
+        buffer: float = 1000.0,
+    ):
+        """Read forcing geodataset."""
+
+        refdate, tstart, tstop = self.get_model_time()  # time slice
+
+        if (
+            forcing_geodataset_fn is not None
+            and self.data_catalog[forcing_geodataset_fn].data_type == "GeoDataset"
+        ):
+            da = self.data_catalog.get_geodataset(
+                forcing_geodataset_fn,
+                geom=self.region.buffer(
+                    buffer
+                ),  # only select data within region of interest (large region for forcing)
+                variables=[forcing_name],
+                time_tuple=(
+                    tstart,
+                    tstop
+                    + timedelta(
+                        seconds=1
+                    ),  # extend with 1 seconds to include the last timestep
+                ),
+            ).rename(forcing_name)
+            # error if time mismatch
+            if np.logical_and(
+                pd.to_datetime(da.time.values[0]) == pd.to_datetime(tstart),
+                pd.to_datetime(da.time.values[-1]) == pd.to_datetime(tstop),
+            ):
+                pass
+            else:
+                self.logger.error(
+                    "forcing has different start and end time. Please check the forcing file. support yyyy-mm-dd HH:MM:SS. "
+                )
+            # reproject if needed and convert to location
+            if da.vector.crs != self.crs:
+                da = da.vector.to_crs(self.crs)
+            # get geom
+            gdf = da.vector.to_gdf(reducer=np.mean)
+        elif (
+            forcing_geodataset_fn is not None
+            and self.data_catalog[forcing_geodataset_fn].data_type == "GeoDataFrame"
+        ):
+            gdf = self.data_catalog.get_geodataframe(
+                forcing_geodataset_fn,
+                geom=self.region.buffer(self._network_snap_offset),
+            )
+            da = None
+        else:
+            gdf = None
+            da = None
+        return gdf, da
+
     def _snap_geom_to_branches_and_drop_nonsnapped(
         self, branches: gpd.GeoDataFrame, geoms: gpd.GeoDataFrame, snap_offset=0.0
     ):
@@ -2602,6 +2659,83 @@ def __set_map_parameters_based_on_variable(
                     relsize = 1.01
                 self._MAPS[var]["averagingrelsize"] = relsize
 
+    def setup_2dboundary_from_lines(
+        self,
+        boundaries_geodataset_fn: str = None,
+        boundary_value: float = 0.0,
+        boundary_type: str = "waterlevel",
+        tolerance: float = 3.0,
+    ):
+        """
+        Prepares the 2D boundaries from geodataset of line geometries.
+        E.g. `boundary_value` m3/s `boundary_type` boundary for all lines in `boundaries_geodataset_fn`.
+
+        Use ``boundaries_geodataset_fn`` to set the boundary timeseries from a geodataset
+        of line geometries.
+        Support also geodataframe of line geometries in combination of ``boundary_value`` and ``boundary_type``.
+
+        Only lines that are within a max distance defined in ``tolerance`` are used.
+
+        The timeseries are clipped to the model time based on the model config
+        tstart and tstop entries.
+        If the timeseries has missing values, the constant ``boundary_value`` will be used.
+
+        Adds/Updates model layers:
+            * **boundary2d_{boundary_name}** forcing: 2D boundaries DataArray
+
+        Parameters
+        ----------
+        boundaries_geodataset_fn : str, Path
+            Path or data source name for geospatial point location file.
+            * Required variables if geodataset is provided [``boundary_type``]
+            NOTE: Use universal datetime format e.g. yyyy-mm-dd to avoid ambiguity when using a csv timeseries.
+        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".
+        tolerance: float, optional
+            Search tolerance factor between boundary polyline and grid cells.
+            Unit: in cell size units (i.e., not meters)
+            By default, 3.0
+
+        """
+        self.logger.info("Preparing 2D boundaries.")
+
+        if boundary_type == "waterlevel":
+            boundary_unit = "m"
+        if boundary_type == "discharge":
+            boundary_unit = "m3/s"
+
+        # 1. read boundary geodataset
+        gdf_bnd, da_bnd = self._read_forcing_geodataset(
+            boundaries_geodataset_fn,
+            boundary_type,
+            buffer=self.res * tolerance,
+        )
+        if len(gdf_bnd) == 0:
+            return None
+
+        # set boundary index -> name of the pli line
+        gdf_bnd.index = [f"{boundary_type}bnd_{i}" for i in gdf_bnd.index]
+
+        # 2. Compute boundary dataarray
+        da_out = workflows.compute_forcing_values_lines(
+            gdf=gdf_bnd,
+            da=da_bnd,
+            forcing_value=boundary_value,
+            forcing_type=boundary_type + "bnd",
+            forcing_unit=boundary_unit,
+            logger=self.logger,
+        )
+
+        # 3. set laterals
+        self.set_forcing(da_out, name=f"boundary2d_lines_{boundary_type+'bnd'}")
+
+        # adjust parameters
+        self.set_config("geometry.openboundarytolerance", tolerance)
+
     def setup_2dboundary(
         self,
         boundaries_fn: str = None,
@@ -3393,8 +3527,8 @@ def read_mesh(self):
     def write_mesh(self, write_gui=True):
         """Write 1D branches and 2D mesh at <root/dflowfm/fm_net.nc> in model ready format."""
         self._assert_write_mode()
-        savedir = dirname(join(self.root, self._config_fn)) # from this branch
-        savedir = join(self.root, "dflowfm") # from main
+        savedir = dirname(join(self.root, self._config_fn))  # from this branch
+        savedir = join(self.root, "dflowfm")  # from main
         mesh_filename = "fm_net.nc"
 
         # write mesh