diff --git a/docs/api.rst b/docs/api.rst
index a0cc5346..9eb27f7d 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -44,12 +44,16 @@ Setup components
    SfincsModel.setup_river_outflow
    SfincsModel.setup_observation_points
    SfincsModel.setup_structures
+   SfincsModel.setup_drainage_structures
    SfincsModel.setup_waterlevel_forcing
    SfincsModel.setup_waterlevel_bnd_from_mask
    SfincsModel.setup_discharge_forcing
    SfincsModel.setup_discharge_forcing_from_grid
    SfincsModel.setup_precip_forcing
    SfincsModel.setup_precip_forcing_from_grid
+   SfincsModel.setup_pressure_forcing_from_grid
+   SfincsModel.setup_wind_forcing
+   SfincsModel.setup_wind_forcing_from_grid
    SfincsModel.setup_tiles
 
 Plot methods
@@ -182,6 +186,8 @@ Input/Output methods
    utils.write_xyn
    utils.read_geoms
    utils.write_geoms
+   utils.read_drn
+   utils.write_drn
    utils.read_sfincs_map_results
    utils.read_sfincs_his_results
 
diff --git a/docs/changelog.rst b/docs/changelog.rst
index 1f36aad0..57249db2 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -8,17 +8,24 @@ The format is based on `Keep a Changelog`_, and this project adheres to
 v1.0.1 (unreleased)
 ===================
 
-
 Bugfix
 ------
 - bugfix in `SfincsModel.write_forcing` to ensure all NetCDF files are written instead of only the first one. PR #86
 - bugfix in `SfincsModel.read_config` & `SfincsInput.read` for relative paths in inp file. PR #88
 - bugfix in `SfincsModel.setup_subgrid` to ensure that a 'big geotiff' will be written by default when 'write_dep_tif' or 'write_man_tif' are True
 - fix memory issues caused by rasterizing the model region and reprojecting before clipping of rasters. PR #94 
+- bugfix in `Sfincs.read_forcing` when combining attributes from the locations stored in the gis folder with the actual forcing locations. PR #99
+- bugfix in `SfincsModel.setup_discharge_from_grid` when snapping based on upstream area in case a src points is outside of the uparea grid domain. PR #99
 
 New
 -----------
-- `SfincsModel.setup_cn_infiltration_with_kr` to setup three layers related to the curve number (maximum and effective infiltration capacity; seff and smax) and recovery rate (kr)
+- `SfincsModel.setup_cn_infiltration_with_kr` to setup three layers related to the curve number (maximum and effective infiltration capacity; seff and smax) and recovery rate (kr). PR#87
+- `SfincsModelsetup_drainage_structures` to setup drainage structures (pumps,culverts) from a geodataframe. PR#90
+- Added `SfincsModel.setup_wind_forcing`, `SfincsModel.setup_wind_forcing_from_grid` and `SfincsModel.setup_pressure_forcing_from_grid` methods to easily add wind and pressure forcing.  PR #92
+- `SfincsModel.read_config` allows to use a template input file from a directory different than the model root. PR #102
+- Added the option to use landuse/landcover data combined with a reclass table to `SfincsModel.setup_constant_infiltration`.  PR #103
+- Enabled to provide locations only (so no timeseries) for `SfincsModel.setup_waterlevel_forcing` and `SfincsModel.setup_discharge_forcing` PR #104
+- New optional buffer argument in  `SfincsModel.setup_discharge_forcing` to select gauges around boundary only. PR #104
 
 Changed
 -------
@@ -206,4 +213,4 @@ Documentation
 
 
 .. _Keep a Changelog: https://keepachangelog.com/en/1.0.0/
-.. _Semantic Versioning: https://semver.org/spec/v2.0.0.html
\ No newline at end of file
+.. _Semantic Versioning: https://semver.org/spec/v2.0.0.html
diff --git a/docs/user_guide/sfincs_model_setup.rst b/docs/user_guide/sfincs_model_setup.rst
index f3c5023e..e5067c5b 100644
--- a/docs/user_guide/sfincs_model_setup.rst
+++ b/docs/user_guide/sfincs_model_setup.rst
@@ -36,7 +36,7 @@ For more information about each file, see the `SFINCS documentation <https://sfi
    * - :py:attr:`~hydromt_sfincs.SfincsModel.subgrid`
      - sbgfile
    * - :py:attr:`~hydromt_sfincs.SfincsModel.geoms`
-     - obsfile, thdfile, weirfile
+     - obsfile, thdfile, weirfile, drnfile
    * - :py:attr:`~hydromt_sfincs.SfincsModel.forcing`
      - bndfile, bzsfile, srcfile, disfile, precipfile, netbndbzsbzifile, netsrcdisfile, netamprfile, netampfile, netamuamvfile
    * - :py:attr:`~hydromt_sfincs.SfincsModel.states`
@@ -122,6 +122,8 @@ Geoms setup methods
      - This component adds observation points to the model (obsfile). 
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_structures`
      - This component adds line element structures to the model (thdfile, weirfile).
+   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_drainage_structures`
+     - This component adds drainage structures (pump, culvert) to the model (drnfile).     
 
 Forcing setup methods
 ---------------------
@@ -145,9 +147,15 @@ Forcing setup methods
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_discharge_forcing_from_grid`
      - This component adds discharge forcing (srcfile, disfile) based on a gridded discharge dataset.
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_precip_forcing`
-     - This component adds spatially uniform precipitation forcing (precipfile).
+     - This component adds spatially uniform precipitation forcing from timeseries/constants (precipfile).
    * - :py:func:`~hydromt_sfincs.SfincsModel.setup_precip_forcing_from_grid`
      - This component adds precipitation forcing from a gridded spatially varying data source (netamprfile).
+   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_pressure_forcing_from_grid`
+     - This component adds pressure forcing from a gridded spatially varying data source (netampfile).
+   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_wind_forcing`
+     - This component adds spatially uniform wind forcing from timeseries/constants (wndfile).
+   * - :py:func:`~hydromt_sfincs.SfincsModel.setup_wind_forcing_from_grid`
+     - This component adds wind forcing from a gridded spatially varying data source (netamuamvfile).               
 
 Other setup methods
 -------------------
diff --git a/hydromt_sfincs/plots.py b/hydromt_sfincs/plots.py
index b439f559..085c19dd 100644
--- a/hydromt_sfincs/plots.py
+++ b/hydromt_sfincs/plots.py
@@ -52,8 +52,15 @@ def plot_forcing(forcing: Dict, **kwargs):
         unit = da.attrs.get("unit", "")
         prefix = ""
         if da.ndim == 3:
-            da = da.mean(dim=[da.raster.x_dim, da.raster.y_dim])
-            prefix = "mean "
+            if name.startswith("press"):
+                da = da.min(dim=[da.raster.x_dim, da.raster.y_dim])
+                prefix = "min "
+            elif name.startswith("wind_u") or name.startswith("wind_v"):
+                da = da.max(dim=[da.raster.x_dim, da.raster.y_dim])
+                prefix = "max "
+            else:
+                da = da.mean(dim=[da.raster.x_dim, da.raster.y_dim])
+                prefix = "mean "
         # convert to Single index dataframe (bar plots don't work with xarray)
         df = da.to_pandas()
         if isinstance(df.index, pd.MultiIndex):
@@ -62,6 +69,20 @@ def plot_forcing(forcing: Dict, **kwargs):
         df.index = mdates.date2num(df.index)
         if name.startswith("precip"):
             axes[i].bar(df.index, df.values, facecolor="darkblue")
+        elif (
+            name.startswith("press")
+            or name.startswith("wind_u")
+            or name.startswith("wind_v")
+        ):
+            df.plot.line(ax=axes[i])
+        elif name.startswith("wnd"):
+            df.plot(ax=axes[i], kind="line", secondary_y="dir", legend=False)
+            # set tick color for y-axis of variable 1
+            axes[i].tick_params(axis="y", labelcolor="C0")
+            axes[i].right_ax.set_ylabel("Wind direction [degrees]")
+            # set tick color and label for y-axis of variable 2
+            axes[i].right_ax.tick_params(axis="y", labelcolor="C1")
+
         else:
             df.plot.line(ax=axes[i]).legend(
                 title="index",
diff --git a/hydromt_sfincs/sfincs.py b/hydromt_sfincs/sfincs.py
index b3f965ce..d38a0039 100644
--- a/hydromt_sfincs/sfincs.py
+++ b/hydromt_sfincs/sfincs.py
@@ -18,9 +18,9 @@
 from hydromt.models.model_grid import GridModel
 from hydromt.raster import RasterDataArray
 from hydromt.vector import GeoDataArray, GeoDataset
+from hydromt.workflows.forcing import da_to_timedelta
 from pyproj import CRS
-from shapely.geometry import box
-from shapely.geometry import Polygon
+from shapely.geometry import box, LineString, MultiLineString, Polygon
 
 from . import DATADIR, plots, utils, workflows
 from .regulargrid import RegularGrid
@@ -40,6 +40,7 @@ class SfincsModel(GridModel):
         "observation_points": "obs",
         "weirs": "weir",
         "thin_dams": "thd",
+        "drainage_structures": "drn",
     }  # parsed to dict of geopandas.GeoDataFrame
     _FORCING_1D = {
         # timeseries (can be multiple), locations tuple
@@ -47,6 +48,7 @@ class SfincsModel(GridModel):
         "waves": (["bzi"], "bnd"),
         "discharge": (["dis"], "src"),
         "precip": (["precip"], None),
+        "wind": (["wnd"], None),
         "wavespectra": (["bhs", "btp", "bwd", "bds"], "bwv"),
         "wavemaker": (["whi", "wti", "wst"], "wvp"),  # TODO check names and test
     }
@@ -54,9 +56,12 @@ class SfincsModel(GridModel):
         # 2D forcing sfincs name, rename tuple
         "waterlevel": ("netbndbzsbzi", {"zs": "bzs", "zi": "bzi"}),
         "discharge": ("netsrcdis", {"discharge": "dis"}),
-        "precip": ("netampr", {"Precipitation": "precip"}),
-        "press": ("netamp", {"barometric_pressure": "press"}),
-        "wind": ("netamuamv", {"eastward_wind": "wind_u", "northward_wind": "wind_v"}),
+        "precip_2d": ("netampr", {"Precipitation": "precip_2d"}),
+        "press_2d": ("netamp", {"barometric_pressure": "press_2d"}),
+        "wind_2d": (
+            "netamuamv",
+            {"eastward_wind": "wind_u", "northward_wind": "wind_v"},
+        ),
     }
     _FORCING_SPW = {"spiderweb": "spw"}  # TODO add read and write functions
     _MAPS = ["msk", "dep", "scs", "manning", "qinf", "smax", "seff", "kr"]
@@ -78,6 +83,11 @@ class SfincsModel(GridModel):
         "bzi": {"standard_name": "wave height", "unit": "m"},
         "dis": {"standard_name": "discharge", "unit": "m3.s-1"},
         "precip": {"standard_name": "precipitation", "unit": "mm.hr-1"},
+        "precip_2d": {"standard_name": "precipitation", "unit": "mm.hr-1"},
+        "press_2d": {"standard_name": "barometric pressure", "unit": "Pa"},
+        "wind_u": {"standard_name": "eastward wind", "unit": "m/s"},
+        "wind_v": {"standard_name": "northward wind", "unit": "m/s"},
+        "wnd": {"standard_name": "wind", "unit": "m/s"},
     }
 
     def __init__(
@@ -943,7 +953,13 @@ def setup_river_outflow(
             self.set_geoms(gdf_riv, name="rivers_outflow")
 
     # Function to create constant spatially varying infiltration
-    def setup_constant_infiltration(self, qinf, reproj_method="average"):
+    def setup_constant_infiltration(
+        self,
+        qinf=None,
+        lulc=None,
+        reclass_table=None,
+        reproj_method="average",
+    ):
         """Setup spatially varying constant infiltration rate (qinffile).
 
         Adds model layers:
@@ -954,18 +970,48 @@ def setup_constant_infiltration(self, qinf, reproj_method="average"):
         ----------
         qinf : str, Path, or RasterDataset
             Spatially varying infiltration rates [mm/hr]
+        lulc: str, Path, or RasterDataset
+            Landuse/landcover data set
+        reclass_table: str, Path, or pd.DataFrame
+            Reclassification table to convert landuse/landcover to infiltration rates [mm/hr]
         reproj_method : str, optional
             Resampling method for reprojecting the infiltration data to the model grid.
             By default 'average'. For more information see, :py:meth:`hydromt.raster.RasterDataArray.reproject_like`
         """
 
         # get infiltration data
-        da_inf = self.data_catalog.get_rasterdataset(
-            qinf, bbox=self.mask.raster.transform_bounds(4326), buffer=10
-        )
-        da_inf = da_inf.raster.mask_nodata()  # set nodata to nan
+        if qinf is not None:
+            da_inf = self.data_catalog.get_rasterdataset(
+                qinf,
+                bbox=self.mask.raster.transform_bounds(4326),
+                buffer=10,
+            )
+        elif lulc is not None:
+            # landuse/landcover should always be combined with mapping
+            if reclass_table is None:
+                raise IOError(
+                    f"Infiltration mapping file should be provided for {lulc}"
+                )
+            da_lulc = self.data_catalog.get_rasterdataset(
+                lulc,
+                bbox=self.mask.raster.transform_bounds(4326),
+                buffer=10,
+                variables=["lulc"],
+            )
+            df_map = self.data_catalog.get_dataframe(
+                reclass_table,
+                variables=["qinf"],
+                index_col=0,  # driver kwargs
+            )
+            # reclassify
+            da_inf = da_lulc.raster.reclassify(df_map)["qinf"]
+        else:
+            raise ValueError(
+                "Either qinf or lulc must be provided when setting up constant infiltration."
+            )
 
         # reproject infiltration data to model grid
+        da_inf = da_inf.raster.mask_nodata()  # set nodata to nan
         da_inf = da_inf.raster.reproject_like(self.mask, method=reproj_method)
 
         # check on nan values
@@ -1339,6 +1385,88 @@ def setup_structures(
         self.set_geoms(gdf, stype)
         self.set_config(f"{stype}file", f"sfincs.{stype}")
 
+    def setup_drainage_structures(
+        self,
+        structures: Union[str, Path, gpd.GeoDataFrame],
+        stype: str = "pump",
+        discharge: float = 0.0,
+        merge: bool = True,
+        **kwargs,
+    ):
+        """Setup drainage structures.
+
+        Adds model layer:
+        * **drn** geom: drainage pump or culvert
+
+        Parameters
+        ----------
+        structures : str, Path
+            Path, data source name, or geopandas object to structure line geometry file.
+            The line should consist of only 2 points (else first and last points are used), ordered from up to downstream.
+            The "type" (1 for pump and 2 for culvert), "par1" ("discharge" also accepted) variables are optional.
+            If "type" or "par1" are not provided, they are based on stype or discharge arguments.
+        stype : {'pump', 'culvert'}, optional
+            Structure type, by default "pump". stype is converted to integer "type" to match with SFINCS expectations.
+        discharge : float, optional
+            Discharge of the structure, by default 0.0. For culverts, this is the maximum discharge,
+            since actual discharge depends on waterlevel gradient
+        merge : bool, optional
+            If True, merge with existing drainage structures, by default True.
+        """
+
+        stype = stype.lower()
+        svalues = {"pump": 1, "culvert": 2}
+        if stype not in svalues:
+            raise ValueError('stype must be one of "pump", "culvert"')
+        svalue = svalues[stype]
+
+        # read, clip and reproject
+        gdf_structures = self.data_catalog.get_geodataframe(
+            structures, geom=self.region, **kwargs
+        ).to_crs(self.crs)
+
+        # check if type (int) is present in gdf, else overwrite from args
+        # TODO also add check if type is interger?
+        if "type" not in gdf_structures:
+            gdf_structures["type"] = svalue
+        # if discharge is provided, rename to par1
+        if "discharge" in gdf_structures:
+            gdf_structures = gdf_structures.rename(columns={"discharge": "par1"})
+
+        # add par1, par2, par3, par4, par5 if not present
+        # NOTE only par1 is used in the model
+        if "par1" not in gdf_structures:
+            gdf_structures["par1"] = discharge
+        if "par2" not in gdf_structures:
+            gdf_structures["par2"] = 0
+        if "par3" not in gdf_structures:
+            gdf_structures["par3"] = 0
+        if "par4" not in gdf_structures:
+            gdf_structures["par4"] = 0
+        if "par5" not in gdf_structures:
+            gdf_structures["par5"] = 0
+
+        # multi to single lines
+        lines = gdf_structures.explode(column="geometry").reset_index(drop=True)
+        # get start [0] and end [1] points
+        endpoints = lines.boundary.explode().unstack()
+        # merge start and end points into a single linestring
+        gdf_structures["geometry"] = endpoints.apply(
+            lambda x: LineString(x.values.tolist()), axis=1
+        )
+
+        # combine with existing structures if present
+        if merge and "drn" in self.geoms:
+            gdf0 = self._geoms.pop("drn")
+            gdf_structures = gpd.GeoDataFrame(
+                pd.concat([gdf_structures, gdf0], ignore_index=True)
+            )
+            self.logger.info(f"Adding {stype} structures to existing structures.")
+
+        # set structures
+        self.set_geoms(gdf_structures, "drn")
+        self.set_config("drnfile", f"sfincs.drn")
+
     ### FORCING
     def set_forcing_1d(
         self,
@@ -1484,6 +1612,7 @@ def setup_waterlevel_forcing(
             Path, data source name, or pandas data object for tabular timeseries.
         locations: str, Path, gpd.GeoDataFrame, optional
             Path, data source name, or geopandas object for bnd point locations.
+            It should contain a 'index' column matching the column names in `timeseries`.
         offset: str, Path, xr.Dataset, float, optional
             Path, data source name, constant value or xarray raster data for gridded offset
             between vertical reference of elevation and waterlevel data,
@@ -1527,14 +1656,14 @@ def setup_waterlevel_forcing(
                 index_col=0,
             )
             df_ts.columns = df_ts.columns.map(int)  # parse column names to integers
-        else:
-            raise ValueError("Either geodataset or timeseries must be provided")
 
-        # optionally read location data (if not already read from geodataset)
+        # read location data (if not already read from geodataset)
         if gdf_locs is None and locations is not None:
             gdf_locs = self.data_catalog.get_geodataframe(
                 locations, geom=region, buffer=buffer, crs=self.crs
             ).to_crs(self.crs)
+            if "index" in gdf_locs.columns:
+                gdf_locs = gdf_locs.set_index("index")
         elif gdf_locs is None and "bzs" in self.forcing:
             gdf_locs = self.forcing["bzs"].vector.to_gdf()
         elif gdf_locs is None:
@@ -1559,7 +1688,7 @@ def setup_waterlevel_forcing(
             )
 
         # set/ update forcing
-        self.set_forcing_1d(df_ts, gdf_locs, name="bzs", merge=merge)
+        self.set_forcing_1d(df_ts=df_ts, gdf_locs=gdf_locs, name="bzs", merge=merge)
 
     def setup_waterlevel_bnd_from_mask(
         self,
@@ -1606,7 +1735,12 @@ def setup_waterlevel_bnd_from_mask(
         self.set_forcing_1d(gdf_locs=gdf, name="bzs", merge=merge)
 
     def setup_discharge_forcing(
-        self, geodataset=None, timeseries=None, locations=None, merge=True
+        self,
+        geodataset=None,
+        timeseries=None,
+        locations=None,
+        merge=True,
+        buffer: float = None,
     ):
         """Setup discharge forcing.
 
@@ -1629,8 +1763,12 @@ def setup_discharge_forcing(
             Path, data source name, or pandas data object for tabular timeseries.
         locations: str, Path, gpd.GeoDataFrame, optional
             Path, data source name, or geopandas object for bnd point locations.
+            It should contain a 'index' column matching the column names in `timeseries`.
         merge : bool, optional
             If True, merge with existing forcing data, by default True.
+        buffer: float, optional
+            Buffer [m] around model boundary within the model region
+            select discharge gauges, by default None.
 
         See Also
         --------
@@ -1638,12 +1776,16 @@ def setup_discharge_forcing(
         """
         gdf_locs, df_ts = None, None
         tstart, tstop = self.get_model_time()  # model time
+        # buffer
+        region = self.region
+        if buffer is not None:  # TODO this assumes the model crs is projected
+            region = region.boundary.buffer(buffer).clip(self.region)
         # read waterlevel data from geodataset or geodataframe
         if geodataset is not None:
             # read and clip data in time & space
             da = self.data_catalog.get_geodataset(
                 geodataset,
-                bbox=self.mask.raster.transform_bounds(4326),
+                geom=region,
                 variables=["discharge"],
                 time_tuple=(tstart, tstop),
                 crs=self.crs,
@@ -1659,21 +1801,21 @@ def setup_discharge_forcing(
                 index_col=0,
             )
             df_ts.columns = df_ts.columns.map(int)  # parse column names to integers
-        else:
-            raise ValueError("Either geodataset or timeseries must be provided")
 
-        # optionally read location data (if not already read from geodataset)
+        # read location data (if not already read from geodataset)
         if gdf_locs is None and locations is not None:
             gdf_locs = self.data_catalog.get_geodataframe(
-                locations, geom=self.region, crs=self.crs
+                locations, geom=region, crs=self.crs
             ).to_crs(self.crs)
+            if "index" in gdf_locs.columns:
+                gdf_locs = gdf_locs.set_index("index")
         elif gdf_locs is None and "dis" in self.forcing:
             gdf_locs = self.forcing["dis"].vector.to_gdf()
         elif gdf_locs is None:
             raise ValueError("No discharge boundary (src) points provided.")
 
         # set/ update forcing
-        self.set_forcing_1d(df_ts, gdf_locs, name="dis", merge=merge)
+        self.set_forcing_1d(df_ts=df_ts, gdf_locs=gdf_locs, name="dis", merge=merge)
 
     def setup_discharge_forcing_from_grid(
         self,
@@ -1779,7 +1921,7 @@ def setup_discharge_forcing_from_grid(
         self.set_geoms(ds_snapped.vector.to_gdf(), "src_snapped")
 
     def setup_precip_forcing_from_grid(
-        self, precip=None, dst_res=None, aggregate=False, **kwargs
+        self, precip, dst_res=None, aggregate=False, **kwargs
     ):
         """Setup precipitation forcing from a gridded spatially varying data source.
 
@@ -1838,20 +1980,22 @@ def setup_precip_forcing_from_grid(
                 dst_crs=self.crs, dst_res=dst_res, **kwargs
             ).fillna(0)
 
-            # resample in time
-            precip_out = hydromt.workflows.resample_time(
-                precip_out,
-                freq=pd.to_timedelta("1H"),
-                conserve_mass=True,
-                upsampling="bfill",
-                downsampling="sum",
-                logger=self.logger,
-            ).rename("precip")
+            # only resample in time if freq < 1H, else keep input values
+            if da_to_timedelta(precip_out) < pd.to_timedelta("1H"):
+                precip_out = hydromt.workflows.resample_time(
+                    precip_out,
+                    freq=pd.to_timedelta("1H"),
+                    conserve_mass=True,
+                    upsampling="bfill",
+                    downsampling="sum",
+                    logger=self.logger,
+                )
+            precip_out = precip_out.rename("precip_2d")
 
             # add to forcing
-            self.set_forcing(precip_out, name="precip")
+            self.set_forcing(precip_out, name="precip_2d")
 
-    def setup_precip_forcing(self, timeseries):
+    def setup_precip_forcing(self, timeseries=None, magnitude=None):
         """Setup spatially uniform precipitation forcing (precip).
 
         Adds model layers:
@@ -1865,15 +2009,26 @@ def setup_precip_forcing(self, timeseries):
             and location IDs in the first row,
             see :py:meth:`hydromt.open_timeseries_from_table`, for details.
             Note: tabulated timeseries files cannot yet be set through the data_catalog yml file.
+        magnitude: float
+            Precipitation magnitude [mm/hr] to use if no timeseries is provided.
         """
         tstart, tstop = self.get_model_time()
-        df_ts = self.data_catalog.get_dataframe(
-            timeseries,
-            time_tuple=(tstart, tstop),
-            # kwargs below only applied if timeseries not in data catalog
-            parse_dates=True,
-            index_col=0,
-        )
+        if timeseries is not None:
+            df_ts = self.data_catalog.get_dataframe(
+                timeseries,
+                time_tuple=(tstart, tstop),
+                # kwargs below only applied if timeseries not in data catalog
+                parse_dates=True,
+                index_col=0,
+            )
+        elif magnitude is not None:
+            times = pd.date_range(*self.get_model_time(), freq="10T")
+            df_ts = pd.DataFrame(
+                index=times, data=np.full((len(times), 1), magnitude, dtype=float)
+            )
+        else:
+            raise ValueError("Either timeseries or magnitude must be provided")
+
         if isinstance(df_ts, pd.DataFrame):
             df_ts = df_ts.squeeze()
         if not isinstance(df_ts, pd.Series):
@@ -1882,6 +2037,173 @@ def setup_precip_forcing(self, timeseries):
         df_ts.index.name = "time"
         self.set_forcing(df_ts.to_xarray(), name="precip")
 
+    def setup_pressure_forcing_from_grid(
+        self, press, dst_res=None, fill_value=101325, **kwargs
+    ):
+        """Setup pressure forcing from a gridded spatially varying data source.
+
+        Adds one model layer:
+
+        * **netampfile** forcing: distributed barometric pressure [Pa]
+
+        Parameters
+        ----------
+        press, str, Path, xr.Dataset, xr.DataArray
+            Path to pressure rasterdataset netcdf file or xarray dataset.
+
+            * Required variables: ['press' (Pa)]
+            * Required coordinates: ['time', 'y', 'x']
+
+        dst_res: float
+            output resolution (m), by default None and computed from source data.
+
+        fill_value: float
+            value to use when no data is available.
+            Standard atmospheric pressure (101325 Pa) is used if no value is given.
+        """
+        # get data for model domain and config time range
+        press = self.data_catalog.get_rasterdataset(
+            press,
+            geom=self.region,
+            buffer=2,
+            time_tuple=self.get_model_time(),
+            variables=["press"],
+        )
+
+        # reproject to model utm crs
+        # NOTE: currently SFINCS errors (stack overflow) on large files,
+        # downscaling to model grid is not recommended
+        kwargs0 = dict(align=dst_res is not None, method="nearest_index")
+        kwargs0.update(kwargs)
+        meth = kwargs0["method"]
+        self.logger.debug(f"Resample pressure using {meth}.")
+        press_out = press.raster.reproject(
+            dst_crs=self.crs, dst_res=dst_res, **kwargs
+        ).fillna(fill_value)
+
+        # only resample in time if freq < 1H, else keep input values
+        if da_to_timedelta(press_out) < pd.to_timedelta("1H"):
+            press_out = hydromt.workflows.resample_time(
+                press_out,
+                freq=pd.to_timedelta("1H"),
+                conserve_mass=False,
+                upsampling="interpolate",
+                downsampling="interpolate",
+                logger=self.logger,
+            )
+
+        press_out = press_out.rename("press_2d")
+
+        # add to forcing
+        self.set_forcing(press_out, name="press_2d")
+
+    def setup_wind_forcing_from_grid(self, wind, dst_res=None, **kwargs):
+        """Setup pressure forcing from a gridded spatially varying data source.
+
+        Adds one model layer:
+
+        * **netamuamv** forcing: distributed wind [m/s]
+
+        Parameters
+        ----------
+        wind, str, Path, xr.Dataset
+            Path to wind rasterdataset (including eastward and northward components) netcdf file or xarray dataset.
+
+            * Required variables: ['wind_u' (m/s), 'wind_v' (m/s)]
+            * Required coordinates: ['time', 'y', 'x']
+
+        dst_res: float
+            output resolution (m), by default None and computed from source data.
+        """
+        # get data for model domain and config time range
+        wind = self.data_catalog.get_rasterdataset(
+            wind,
+            geom=self.region,
+            buffer=2,
+            time_tuple=self.get_model_time(),
+            variables=["wind_u", "wind_v"],
+        )
+
+        # reproject to model utm crs
+        # NOTE: currently SFINCS errors (stack overflow) on large files,
+        # downscaling to model grid is not recommended
+        kwargs0 = dict(align=dst_res is not None, method="nearest_index")
+        kwargs0.update(kwargs)
+        meth = kwargs0["method"]
+        self.logger.debug(f"Resample wind using {meth}.")
+
+        wind = wind.raster.reproject(
+            dst_crs=self.crs, dst_res=dst_res, **kwargs
+        ).fillna(0)
+
+        # only resample in time if freq < 1H, else keep input values
+        if da_to_timedelta(wind) < pd.to_timedelta("1H"):
+            wind_out = xr.Dataset()
+            # resample in time
+            for var in wind.data_vars:
+                wind_out[var] = hydromt.workflows.resample_time(
+                    wind[var],
+                    freq=pd.to_timedelta("1H"),
+                    conserve_mass=False,
+                    upsampling="interpolate",
+                    downsampling="interpolate",
+                    logger=self.logger,
+                )
+        else:
+            wind_out = wind
+
+        # add to forcing
+        self.set_forcing(wind_out, name="wind_2d")
+
+    def setup_wind_forcing(self, timeseries=None, magnitude=None, direction=None):
+        """Setup spatially uniform wind forcing (wind).
+
+        Adds model layers:
+
+        * **windfile** forcing: uniform wind magnitude [m/s] and direction [deg]
+
+        Parameters
+        ----------
+        timeseries, str, Path
+            Path to tabulated timeseries csv file with time index in first column,
+            magnitude in second column and direction in third column
+            see :py:meth:`hydromt.open_timeseries_from_table`, for details.
+            Note: tabulated timeseries files cannot yet be set through the data_catalog yml file.
+        magnitude: float
+            Magnitude of the wind [m/s]
+        direction: float
+            Direction where the wind is coming from [deg], e.g. 0 is north, 90 is east, etc.
+        """
+        tstart, tstop = self.get_model_time()
+        if timeseries is not None:
+            df_ts = self.data_catalog.get_dataframe(
+                timeseries,
+                time_tuple=(tstart, tstop),
+                # kwargs below only applied if timeseries not in data catalog
+                parse_dates=True,
+                index_col=0,
+            )
+        elif magnitude is not None and direction is not None:
+            df_ts = pd.DataFrame(
+                index=pd.date_range(*self.get_model_time(), periods=2),
+                data=np.array([[magnitude, direction], [magnitude, direction]]),
+                columns=["mag", "dir"],
+            )
+        else:
+            raise ValueError(
+                "Either timeseries or magnitude and direction must be provided"
+            )
+
+        df_ts.name = "wnd"
+        df_ts.index.name = "time"
+        df_ts.columns.name = "index"
+        da = xr.DataArray(
+            df_ts.values,
+            dims=("time", "index"),
+            coords={"time": df_ts.index, "index": ["mag", "dir"]},
+        )
+        self.set_forcing(da, name="wnd")
+
     def setup_tiles(
         self,
         path: Union[str, Path] = None,
@@ -1992,7 +2314,7 @@ def setup_tiles(
     def plot_forcing(self, fn_out=None, **kwargs):
         """Plot model timeseries forcing.
 
-        For distributed forcing a spatial avarage is plotted.
+        For distributed forcing a spatial avarage, minimum or maximum is plotted.
 
         Parameters
         ----------
@@ -2297,7 +2619,7 @@ def write_subgrid(self):
     def read_geoms(self):
         """Read geometry files and save to `geoms` attribute.
         Known geometry files mentioned in the sfincs.inp configuration file are read,
-        including: bnd/src/obs xy files and thd/weir structure files.
+        including: bnd/src/obs xy(n) files, thd/weir structure files and drn drainage structure files.
 
         If other geojson files are present in a "gis" subfolder folder, those are read as well.
         """
@@ -2316,6 +2638,8 @@ def read_geoms(self):
                     gdf = utils.linestring2gdf(struct, crs=self.crs)
                 elif gname == "obs":
                     gdf = utils.read_xyn(fn, crs=self.crs)
+                elif gname == "drn":
+                    gdf = utils.read_drn(fn, crs=self.crs)
                 else:
                     gdf = utils.read_xy(fn, crs=self.crs)
                 self.set_geoms(gdf, name=gname)
@@ -2359,6 +2683,8 @@ def write_geoms(self, data_vars: Union[List, str] = None):
                         utils.write_geoms(fn, struct, stype=gname)
                     elif gname == "obs":
                         utils.write_xyn(fn, gdf, crs=self.crs)
+                    elif gname == "drn":
+                        utils.write_drn(fn, gdf)
                     else:
                         utils.write_xy(fn, gdf, fmt="%8.2f")
 
@@ -2412,17 +2738,22 @@ def read_forcing(self, data_vars: List = None):
                         self.logger.warning(f"{xy_name}file not found at {xy_fn}")
                 else:
                     gdf = utils.read_xy(xy_fn, crs=self.crs)
-                    # read attribute data from gis files
+                    # read attribute data from gis files; merge based on index
                     gis_fn = join(self.root, "gis", f"{xy_name}.geojson")
                     if isfile(gis_fn):
                         gdf1 = gpd.read_file(gis_fn)
                         if "index" in gdf1.columns:
                             gdf1 = gdf1.set_index("index")
-                            gdf.index = gdf1.index.values
-                            ds = ds.assign_coords(index=gdf1.index.values)
-                        if np.any(gdf1.columns != "geometry"):
-                            gdf = gpd.sjoin(gdf, gdf1, how="left")[gdf1.columns]
-                    # set locations as coordinates dataset
+                        if not np.all(np.isin(gdf.index, gdf1.index)):
+                            self.logger.warning(
+                                f"Index in {xy_name}file does not match {gis_fn}"
+                            )
+                        else:
+                            for col in gdf1.columns:
+                                if col not in gdf.columns:
+                                    gdf[col] = gdf1.loc[gdf.index, col]
+                    # set locations and attributes as coordinates of dataset
+                    ds = ds.assign_coords(index=gdf.index.values)
                     ds = GeoDataset.from_gdf(gdf, ds, index_dim="index")
             # save in self.forcing
             if len(ds) > 1:
@@ -2802,42 +3133,35 @@ def write_vector(
 
     ## model configuration
 
-    def read_config(self, config_fn: str = "sfincs.inp", epsg: int = None) -> None:
+    def read_config(self, config_fn: str = None, epsg: int = None) -> None:
         """Parse config from SFINCS input file.
-        If in write-only mode the config is initialized with default settings.
+        If in write-only mode the config is initialized with default settings
+        unless a path to a template config file is provided.
 
         Parameters
         ----------
         config_fn: str
-            Filename of config file, by default "sfincs.inp".
-            If in a different folder than the model root, the root is updated.
+            Filename of config file, by default None.
         epsg: int
-            EPSG code of the model CRS. Only used if missing in the SFINCS input file, by default None.
+            EPSG code of the model CRS. Only used if missing in the SFINCS input file,
+            by default None.
         """
         inp = SfincsInput()  # initialize with defaults
-        if self._read:  # in read-only or append mode, try reading config_fn
-            if not isfile(config_fn) and not isabs(config_fn) and self._root:
-                # path relative to self.root
+        if config_fn is not None or self._read:
+            if config_fn is None:  # read from default location
+                config_fn = self._config_fn
+            if not isabs(config_fn) and self._root:  # read from model root
                 config_fn = abspath(join(self.root, config_fn))
-            elif isfile(config_fn) and abspath(dirname(config_fn)) != self._root:
-                # new root
-                mode = (
-                    "r+"
-                    if self._write and self._read
-                    else ("w" if self._write else "r")
-                )
-                root = abspath(dirname(config_fn))
-                self.logger.warning(f"updating the model root to: {root}")
-                self.set_root(root=root, mode=mode)
-            else:
+            if not isfile(config_fn):
                 raise IOError(f"SFINCS input file not found {config_fn}")
-            # read config_fn
+            # read inp file
             inp.read(inp_fn=config_fn)
         # overwrite / initialize config attribute
         self._config = inp.to_dict()
         if epsg is not None and "epsg" not in self.config:
-            self.config.update(epsg=epsg)
-        self.update_grid_from_config()  # update grid properties based on sfincs.inp
+            self.set_config("epsg", int(epsg))
+        # update grid properties based on sfincs.inp
+        self.update_grid_from_config()
 
     def write_config(self, config_fn: str = "sfincs.inp"):
         """Write config to <root/config_fn>"""
diff --git a/hydromt_sfincs/utils.py b/hydromt_sfincs/utils.py
index 1f2491ed..86a87a19 100644
--- a/hydromt_sfincs/utils.py
+++ b/hydromt_sfincs/utils.py
@@ -40,6 +40,8 @@
     "write_xyn",
     "read_geoms",
     "write_geoms",
+    "read_drn",
+    "write_drn",
     "gdf2linestring",
     "gdf2polygon",
     "linestring2gdf",
@@ -315,7 +317,7 @@ def write_timeseries(
     data[:, 0] = (df.index - tref).total_seconds()
     # calculate required width for time column; hard coded single decimal precision
     # format for other columns is based on fmt`argument
-    w = int(np.floor(np.log10(data[-1, 0]))) + 3
+    w = int(np.floor(np.log10(abs(data[-1, 0])))) + 3
     fmt_lst = [f"%{w}.1f"] + [fmt for _ in range(df.columns.size)]
     fmt_out = " ".join(fmt_lst)
     with open(fn, "w") as f:
@@ -623,6 +625,98 @@ def read_geoms(fn: Union[str, Path]) -> List[Dict]:
     return feats
 
 
+def write_drn(fn: Union[str, Path], gdf_drainage: gpd.GeoDataFrame, fmt="%.4f") -> None:
+    """Write structure files from list of dictionaries.
+
+    Parameters
+    ----------
+    fn : str, Path
+        Path to structure file.
+    drainage : gpd.GeoDataFrame
+        Dataframe with drainage structure parameters and geometry.
+    fmt : str
+        Format for float values.
+    """
+
+    # expected columns for drainage structures
+    col_names = [
+        "xsnk",
+        "ysnk",
+        "xsrc",
+        "ysrc",
+        "type",
+        "par1",
+        "par2",
+        "par3",
+        "par4",
+        "par5",
+    ]
+
+    gdf = copy.deepcopy(gdf_drainage)
+    # get geometry linestring and convert to xsnk, ysnk, xsrc, ysrc
+    endpoints = gdf.boundary.explode().unstack()
+    gdf["xsnk"] = endpoints[0].x
+    gdf["ysnk"] = endpoints[0].y
+    gdf["xsrc"] = endpoints[1].x
+    gdf["ysrc"] = endpoints[1].y
+    gdf.drop(["geometry"], axis=1, inplace=True)
+
+    # reorder columns based on col_names
+    gdf = gdf[col_names]
+
+    # write to file
+    gdf.to_csv(fn, sep=" ", index=False, header=False, float_format=fmt)
+
+
+def read_drn(fn: Union[str, Path], crs: int = None) -> gpd.GeoDataFrame:
+    """Read drainage structure files to geodataframe.
+
+    Parameters
+    ----------
+    fn : str, Path
+        Path to drainge structure file.
+    crs : int
+        EPSG code for coordinate reference system.
+
+    Returns
+    -------
+    gpd.GeoDataFrame
+        Dataframe with drainage structure parameters and geometry.
+    """
+
+    # expected columns for drainage structures
+    col_names = [
+        "xsnk",
+        "ysnk",
+        "xsrc",
+        "ysrc",
+        "type",
+        "par1",
+        "par2",
+        "par3",
+        "par4",
+        "par5",
+    ]
+
+    # read structure file
+    df = pd.read_csv(fn, sep="\\s+", names=col_names)
+
+    # get geometry linestring
+    geom = [
+        LineString([(xsnk, ysnk), (xsrc, ysrc)])
+        for xsnk, ysnk, xsrc, ysrc in zip(
+            df["xsnk"], df["ysnk"], df["xsrc"], df["ysrc"]
+        )
+    ]
+    df.drop(["xsnk", "ysnk", "xsrc", "ysrc"], axis=1, inplace=True)
+
+    # convert to geodataframe
+    gdf = gpd.GeoDataFrame(df, geometry=geom)
+    if crs is not None:
+        gdf.set_crs(crs, inplace=True)
+    return gdf
+
+
 ## OUTPUT: sfincs_map.nc, sfincs_his.nc ##
 
 
diff --git a/hydromt_sfincs/workflows/curvenumber.py b/hydromt_sfincs/workflows/curvenumber.py
index e523bf85..910ede7c 100644
--- a/hydromt_sfincs/workflows/curvenumber.py
+++ b/hydromt_sfincs/workflows/curvenumber.py
@@ -62,10 +62,7 @@ def scs_recovery_determination(da_landuse, da_HSG, da_Ksat, df_map, da_mask_bloc
     # very high 	100 - Inf
     da_kr = da_Ksat.raster.reproject_like(da_kr, method="average").load()
     da_kr = np.minimum(da_kr, 100)  # not higher than 100
-    da_kr = da_kr * 0.141732
-    # from micrometers per second to inch/hr    (constant)
-    da_kr = np.sqrt(da_kr) / 75
-    # recovery in percentage of Smax per hour   (Eq. 4.36)
+    da_kr = da_kr * 3.6  # from micrometers per second to mm/hr    (constant)
 
     # Ensure no NaNs
     da_smax = da_smax.fillna(0)
diff --git a/hydromt_sfincs/workflows/discharge.py b/hydromt_sfincs/workflows/discharge.py
index b601da1b..3ad1730e 100644
--- a/hydromt_sfincs/workflows/discharge.py
+++ b/hydromt_sfincs/workflows/discharge.py
@@ -54,35 +54,35 @@ def snap_discharge(
         logger.debug(
             f"Snapping {discharge_name} points to best matching uparea cell within wdw (size={wdw})."
         )
+        # find cells in window with smallest difference in uparea
         upa0 = xr.DataArray(gdf[uparea_name], dims=("index"))
         upa_dff = np.abs(
             ds_wdw[uparea_name].where(ds_wdw[uparea_name] > 0).load() - upa0
         )
+        i_wdw = upa_dff.fillna(np.inf).argmin("wdw")
+        # find valid cells based on error criteria
         upa_check = np.logical_or((upa_dff / upa0) <= rel_error, upa_dff <= abs_error)
         valid = np.logical_and(valid, upa_check)
-        # combine valid local cells with best matching windows cells if local cell invalid
-        # i_loc = int((1 + 2 * wdw) ** 2 / 2)  # center cell
-        # i_wdw = upa_dff.argmin("wdw").where(~valid.isel(wdw=i_loc), i_loc).load()
-        # find best matching uparea cell in window
-        i_wdw = upa_dff.argmin("wdw").load()
     else:
         logger.debug(
             f"No {uparea_name} variable found in ds or gdf; "
             f"sampling {discharge_name} points from nearest grid cell."
         )
-        # add distance (measured in cells)
+        # calculate distance to center cell (measured in cells)
         ar_wdw = np.abs(np.arange(-wdw, wdw + 1))
         dist = np.hypot(**np.meshgrid(ar_wdw, ar_wdw)).ravel()
         ds_wdw["dist"] = xr.Variable(
             ("index", "wdw"), np.tile(dist, (ds_wdw["index"].size, 1))
         )
+        # find nearest valid cell in window
         i_wdw = ds_wdw["dist"].where(valid, np.inf).argmin("wdw").load()
+    # filter valid cells
     idx_valid = np.where(valid.isel(wdw=i_wdw).values)[0]
     if idx_valid.size < gdf.index.size:
         logger.warning(
             f"{idx_valid.size}/{gdf.index.size} {discharge_name} points successfully snapped."
         )
     i_wdw = i_wdw.isel(index=idx_valid)
+    # return discharge at valid cells
     ds_out = ds_wdw.isel(wdw=i_wdw.load(), index=idx_valid)
-
-    return ds_out  # .reset_coords()[discharge_name]
+    return ds_out
diff --git a/tests/data/local_data.yml b/tests/data/local_data.yml
index b41119d3..a6bacee6 100644
--- a/tests/data/local_data.yml
+++ b/tests/data/local_data.yml
@@ -16,3 +16,8 @@ weir:
   data_type: GeoDataFrame
   driver: vector
   path: local_data/weir.geojson
+
+drainage:
+  data_type: GeoDataFrame
+  driver: vector
+  path: local_data/drainage.geojson
diff --git a/tests/data/local_data/drainage.geojson b/tests/data/local_data/drainage.geojson
new file mode 100644
index 00000000..b0566430
--- /dev/null
+++ b/tests/data/local_data/drainage.geojson
@@ -0,0 +1,9 @@
+{
+"type": "FeatureCollection",
+"name": "drainage",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } },
+"features": [
+{ "type": "Feature", "properties": { "id": 1, "discharge": 10.0, "type": 1 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 322021.879048462666105, 5045547.683568798005581 ], [ 322378.691038294811733, 5045133.781660593114793 ] ] ] } },
+{ "type": "Feature", "properties": { "id": 2, "discharge": 40.0, "type": 2 }, "geometry": { "type": "MultiLineString", "coordinates": [ [ [ 319160.246890009148046, 5045262.233976932242513 ], [ 319167.383129805792123, 5044976.784385066479445 ] ] ] } }
+]
+}
diff --git a/tests/data/sfincs_test.yml b/tests/data/sfincs_test.yml
index a32011dd..13f7c98f 100644
--- a/tests/data/sfincs_test.yml
+++ b/tests/data/sfincs_test.yml
@@ -48,6 +48,9 @@ setup_structures:
   structures: weir
   stype: thd
 
+setup_drainage_structures:
+  structures: drainage
+
 setup_river_inflow:
   hydrography: merit_hydro
   river_upa: 10
diff --git a/tests/data/sfincs_test/gis/drn.geojson b/tests/data/sfincs_test/gis/drn.geojson
new file mode 100644
index 00000000..afebb1fc
--- /dev/null
+++ b/tests/data/sfincs_test/gis/drn.geojson
@@ -0,0 +1,8 @@
+{
+"type": "FeatureCollection",
+"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::32633" } },
+"features": [
+{ "type": "Feature", "properties": { "id": 1, "type": 1, "par1": 10.0, "par2": 0, "par3": 0, "par4": 0, "par5": 0 }, "geometry": { "type": "LineString", "coordinates": [ [ 322021.879048462666105, 5045547.683568798005581 ], [ 322378.691038294811733, 5045133.781660593114793 ] ] } },
+{ "type": "Feature", "properties": { "id": 2, "type": 2, "par1": 40.0, "par2": 0, "par3": 0, "par4": 0, "par5": 0 }, "geometry": { "type": "LineString", "coordinates": [ [ 319160.246890009148046, 5045262.233976932242513 ], [ 319167.383129805792123, 5044976.784385066479445 ] ] } }
+]
+}
diff --git a/tests/data/sfincs_test/precip.nc b/tests/data/sfincs_test/precip_2d.nc
similarity index 100%
rename from tests/data/sfincs_test/precip.nc
rename to tests/data/sfincs_test/precip_2d.nc
diff --git a/tests/data/sfincs_test/sfincs.drn b/tests/data/sfincs_test/sfincs.drn
new file mode 100644
index 00000000..d06a208d
--- /dev/null
+++ b/tests/data/sfincs_test/sfincs.drn
@@ -0,0 +1,2 @@
+322021.8790 5045547.6836 322378.6910 5045133.7817 1 10.0000 0 0 0 0
+319160.2469 5045262.2340 319167.3831 5044976.7844 2 40.0000 0 0 0 0
diff --git a/tests/data/sfincs_test/sfincs.inp b/tests/data/sfincs_test/sfincs.inp
index e00e1c20..c8428c88 100644
--- a/tests/data/sfincs_test/sfincs.inp
+++ b/tests/data/sfincs_test/sfincs.inp
@@ -52,4 +52,5 @@ outputformat         = net
 cdnrb                = 3
 cdwnd                = 0.0 28.0 50.0
 cdval                = 0.001 0.0025 0.0015
-netamprfile          = precip.nc
+netamprfile          = precip_2d.nc
+drnfile              = sfincs.drn
diff --git a/tests/data/sfincs_test/sfincs_his.nc b/tests/data/sfincs_test/sfincs_his.nc
index 6aabea5b..34d4e889 100644
Binary files a/tests/data/sfincs_test/sfincs_his.nc and b/tests/data/sfincs_test/sfincs_his.nc differ
diff --git a/tests/data/sfincs_test/sfincs_map.nc b/tests/data/sfincs_test/sfincs_map.nc
index 928362e0..edbe093a 100644
Binary files a/tests/data/sfincs_test/sfincs_map.nc and b/tests/data/sfincs_test/sfincs_map.nc differ
diff --git a/tests/test_1model_class.py b/tests/test_1model_class.py
index 87af3ae2..920e60d2 100644
--- a/tests/test_1model_class.py
+++ b/tests/test_1model_class.py
@@ -106,7 +106,7 @@ def test_infiltration(tmpdir):
     assert np.isclose(
         mod1.grid["seff"].where(mod.mask > 0).sum(), 32.929287 * effective
     )
-    assert np.isclose(mod1.grid["kr"].where(mod.mask > 0).sum(), 1.7879527)
+    assert np.isclose(mod1.grid["kr"].where(mod.mask > 0).sum(), 330.588)
 
 
 def test_structs(tmpdir):
@@ -138,6 +138,28 @@ def test_structs(tmpdir):
     assert isfile(join(mod.root, "sfincs.weir"))
 
 
+def test_drainage_structures(tmpdir):
+    root = TESTMODELDIR
+    mod = SfincsModel(root=root, mode="r+")
+    # read
+    mod.set_config("drnfile", "sfincs.drn")
+    mod.read_grid()
+    mod.read_geoms()
+    assert "drn" in mod.geoms
+    nr_drainage_structures = len(mod.geoms["drn"].index)
+    # write drn file only
+    tmp_root = str(tmpdir.join("drainage_struct_test"))
+    mod.set_root(tmp_root, mode="w")
+    mod.write_geoms(data_vars=["drn"])
+    assert isfile(join(mod.root, "sfincs.drn"))
+    assert not isfile(join(mod.root, "sfincs.obs"))
+    fn_drn_gis = join(mod.root, "gis", "drn.geojson")
+    assert isfile(fn_drn_gis)
+    # add more drainage structures
+    mod.setup_drainage_structures(fn_drn_gis, merge=True)
+    assert len(mod.geoms["drn"].index) == nr_drainage_structures * 2
+
+
 def test_results():
     root = TESTMODELDIR
     mod = SfincsModel(root=root, mode="r")