From 97bbc0aca8171f1507e82d898b386a71fdf4d7e7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 16 Oct 2023 14:29:49 +0200
Subject: [PATCH 01/12] Add nhi.py and test

---
 nlmod/read/nhi.py     | 56 +++++++++++++++++++++++++++++++++++++++++++
 tests/test_021_nhi.py | 20 ++++++++++++++++
 2 files changed, 76 insertions(+)
 create mode 100644 nlmod/read/nhi.py
 create mode 100644 tests/test_021_nhi.py

diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
new file mode 100644
index 00000000..f7488f18
--- /dev/null
+++ b/nlmod/read/nhi.py
@@ -0,0 +1,56 @@
+import os
+import requests
+import logging
+import numpy as np
+import rioxarray
+from ..dims.resample import structured_da_to_ds
+
+logger = logging.getLogger(__name__)
+
+def download_file(url, pathname, overwrite=False):
+    filename = url.split("/")[-1]
+    fname = os.path.join(pathname, filename)
+    if overwrite or not os.path.isfile(fname):
+        logger.info(f"Downloading {filename}")
+        r = requests.get(url, allow_redirects=True)
+        with open(fname, "wb") as file:
+            file.write(r.content)
+    return fname
+            
+def add_buisdrainage(ds, pathname=None, cond_var="buisdrain_cond", depth_var="buisdrain_depth", cond_method="average", depth_method="mode"):
+    if pathname is None:
+        pathname = ds.cachedir
+    # download files
+    url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m"
+    url = f"{url_bas}/buisdrain_c_ras25/buisdrain_c_ras25.nc"
+    fname_c = download_file(url, pathname)
+    url = f"{url_bas}/buisdrain_d_ras25/buisdrain_d_ras25.nc"
+    fname_d = download_file(url, pathname)
+    
+    # use method = "average" on conductance, to account for locations without pipe
+    # drainage, where the conductance is 0
+    buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0]
+    cond = 1 / buisdrain_c
+    # set conductance to 0 where resistance is infinite or 0
+    cond = cond.where(~(np.isinf(cond) | np.isnan(cond)), 0.0)
+    cond = cond.rio.write_crs(buisdrain_c.rio.crs)
+    # resample to model grid
+    ds[cond_var] = structured_da_to_ds(cond, ds, method=cond_method)
+    
+    # use mthod = "mode" on depth, to retreive the depth that occurs most in each cell
+    mask_and_scale = False
+    buisdrain_d = rioxarray.open_rasterio(fname_d, mask_and_scale=mask_and_scale)[0]
+    if mask_and_scale:
+        nodata = np.nan
+    else:
+        nodata = buisdrain_d._FillValue
+    # set buisdrain_d to nodata where it is 0
+    mask = buisdrain_d != 0
+    buisdrain_d = buisdrain_d.where(mask, nodata).rio.write_crs(buisdrain_d.rio.crs)
+    # resample to model grid
+    ds[depth_var] = structured_da_to_ds(buisdrain_d, ds, method=depth_method, nodata=nodata)
+    if not mask_and_scale:
+        # set nodata values to NaN
+        ds[depth_var] = ds[depth_var].where(ds[depth_var] != buisdrain_d._FillValue)
+    
+    return ds
\ No newline at end of file
diff --git a/tests/test_021_nhi.py b/tests/test_021_nhi.py
new file mode 100644
index 00000000..8991817a
--- /dev/null
+++ b/tests/test_021_nhi.py
@@ -0,0 +1,20 @@
+import os
+import numpy as np
+import tempfile
+import nlmod
+
+tmpdir = tempfile.gettempdir()
+
+
+def test_buidrainage():
+    model_ws = os.path.join(tmpdir, "buidrain")
+    ds = nlmod.get_ds([110_000, 130_000, 435_000, 445_000], model_ws=model_ws)
+    ds = nlmod.read.nhi.add_buisdrainage(ds)
+
+    # assert that all locations with a specified depth also have a positive conductance
+    mask = ~ds["buisdrain_depth"].isnull()
+    assert np.all(ds["buisdrain_cond"].data[mask] > 0)
+
+    # assert that all locations with a positive conductance also have a specified depth
+    mask = ds["buisdrain_cond"] > 0
+    assert np.all(~np.isnan(ds["buisdrain_depth"].data[mask]))

From 8d9f0d13b83a07eb668604a495b363899197fc9e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 16 Oct 2023 17:37:53 +0200
Subject: [PATCH 02/12] Improve nhi.py

---
 nlmod/dims/resample.py     |  4 +--
 nlmod/gwf/surface_water.py |  4 ++-
 nlmod/read/__init__.py     |  1 +
 nlmod/read/nhi.py          | 58 +++++++++++++++++++++++++++++---------
 4 files changed, 50 insertions(+), 17 deletions(-)

diff --git a/nlmod/dims/resample.py b/nlmod/dims/resample.py
index 1768b9fe..50a17cfe 100644
--- a/nlmod/dims/resample.py
+++ b/nlmod/dims/resample.py
@@ -446,7 +446,7 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
     Parameters
     ----------
     da : xarray.DataArray
-        THe data-array to be resampled, with dimensions x and y.
+        The data-array to be resampled, with dimensions x and y.
     ds : xarray.Dataset
         The model dataset.
     method : string or rasterio.enums.Resampling, optional
@@ -457,7 +457,7 @@ def structured_da_to_ds(da, ds, method="average", nodata=np.NaN):
         method is 'linear' or 'nearest' da.interp() is used. Otherwise
         da.rio.reproject_match() is used. The default is "average".
     nodata : float, optional
-        THe nodata value in input and output. THe default is np.NaN.
+        The nodata value in input and output. The default is np.NaN.
 
     Returns
     -------
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index e9101233..91c3d9f0 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -1024,8 +1024,10 @@ def add_season_timeseries(
     package,
     summer_months=(4, 5, 6, 7, 8, 9),
     filename="season.ts",
-    seasons=("winter", "summer"),
+    seasons=None,
 ):
+    if seasons is None:
+        seasons = ["winter", "summer"]
     tmin = pd.to_datetime(ds.time.start)
     if tmin.month in summer_months:
         ts_data = [(0.0, 0.0, 1.0)]
diff --git a/nlmod/read/__init__.py b/nlmod/read/__init__.py
index 63dbf4f4..efbcfac4 100644
--- a/nlmod/read/__init__.py
+++ b/nlmod/read/__init__.py
@@ -13,6 +13,7 @@
     rws,
     waterboard,
     webservices,
+    nhi,
 )
 from .geotop import get_geotop
 from .regis import get_regis
diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
index f7488f18..ec7e2db5 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -7,8 +7,10 @@
 
 logger = logging.getLogger(__name__)
 
-def download_file(url, pathname, overwrite=False):
-    filename = url.split("/")[-1]
+
+def download_file(url, pathname, filename=None, overwrite=False):
+    if filename is None:
+        filename = url.split("/")[-1]
     fname = os.path.join(pathname, filename)
     if overwrite or not os.path.isfile(fname):
         logger.info(f"Downloading {filename}")
@@ -16,27 +18,50 @@ def download_file(url, pathname, overwrite=False):
         with open(fname, "wb") as file:
             file.write(r.content)
     return fname
-            
-def add_buisdrainage(ds, pathname=None, cond_var="buisdrain_cond", depth_var="buisdrain_depth", cond_method="average", depth_method="mode"):
-    if pathname is None:
-        pathname = ds.cachedir
-    # download files
+
+
+def download_buisdrainage(pathname, overwrite=False):
     url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m"
+
+    # download resistance
     url = f"{url_bas}/buisdrain_c_ras25/buisdrain_c_ras25.nc"
-    fname_c = download_file(url, pathname)
+    fname_c = download_file(url, pathname, overwrite=overwrite)
+
+    # download drain depth
     url = f"{url_bas}/buisdrain_d_ras25/buisdrain_d_ras25.nc"
-    fname_d = download_file(url, pathname)
-    
+    fname_d = download_file(url, pathname, overwrite=overwrite)
+
+    return fname_c, fname_d
+
+
+def add_buisdrainage(
+    ds,
+    pathname=None,
+    cond_var="buisdrain_cond",
+    depth_var="buisdrain_depth",
+    cond_method="average",
+    depth_method="mode",
+):
+    if pathname is None:
+        pathname = ds.cachedir
+    # download files if needed
+    fname_c, fname_d = download_buisdrainage(pathname)
+
+    # make sure crs is set of ds
+    if ds.rio.crs is None:
+        ds = ds.rio.write_crs(28992)
+
     # use method = "average" on conductance, to account for locations without pipe
     # drainage, where the conductance is 0
     buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0]
+    # calculate a conductance (per m2) from a resistance
     cond = 1 / buisdrain_c
     # set conductance to 0 where resistance is infinite or 0
     cond = cond.where(~(np.isinf(cond) | np.isnan(cond)), 0.0)
     cond = cond.rio.write_crs(buisdrain_c.rio.crs)
     # resample to model grid
     ds[cond_var] = structured_da_to_ds(cond, ds, method=cond_method)
-    
+
     # use mthod = "mode" on depth, to retreive the depth that occurs most in each cell
     mask_and_scale = False
     buisdrain_d = rioxarray.open_rasterio(fname_d, mask_and_scale=mask_and_scale)[0]
@@ -48,9 +73,14 @@ def add_buisdrainage(ds, pathname=None, cond_var="buisdrain_cond", depth_var="bu
     mask = buisdrain_d != 0
     buisdrain_d = buisdrain_d.where(mask, nodata).rio.write_crs(buisdrain_d.rio.crs)
     # resample to model grid
-    ds[depth_var] = structured_da_to_ds(buisdrain_d, ds, method=depth_method, nodata=nodata)
+    ds[depth_var] = structured_da_to_ds(
+        buisdrain_d, ds, method=depth_method, nodata=nodata
+    )
     if not mask_and_scale:
         # set nodata values to NaN
         ds[depth_var] = ds[depth_var].where(ds[depth_var] != buisdrain_d._FillValue)
-    
-    return ds
\ No newline at end of file
+
+    # from cm to m
+    ds[depth_var] = ds[depth_var] / 100.0
+
+    return ds

From 1b0053f8c81a0a8b3f787171c6dfdd8a21eb3660 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Tue, 17 Oct 2023 15:30:11 +0200
Subject: [PATCH 03/12] Small changes (unrelated to buisdrainage)

---
 nlmod/dims/layers.py |  9 +++++++++
 nlmod/dims/time.py   |  7 ++++++-
 nlmod/read/knmi.py   | 14 +++++++++++---
 3 files changed, 26 insertions(+), 4 deletions(-)

diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py
index 0580ffdc..af7497ed 100644
--- a/nlmod/dims/layers.py
+++ b/nlmod/dims/layers.py
@@ -684,6 +684,8 @@ def set_model_top(ds, top, min_thickness=0.0):
     """
     if "gridtype" not in ds.attrs:
         raise (KeyError("Make sure the Dataset is built by nlmod"))
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_model_top does not support top with a layer dimension"))
     if isinstance(top, (float, int)):
         top = xr.full_like(ds["top"], top)
     if not top.shape == ds["top"].shape:
@@ -707,6 +709,8 @@ def set_model_top(ds, top, min_thickness=0.0):
 def set_layer_top(ds, layer, top):
     """Set the top of a layer."""
     assert layer in ds.layer
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_layer_top does not support top with a layer dimension"))
     lay = np.where(ds.layer == layer)[0][0]
     if lay == 0:
         # change the top of the model
@@ -730,6 +734,8 @@ def set_layer_top(ds, layer, top):
 def set_layer_botm(ds, layer, botm):
     """Set the bottom of a layer."""
     assert layer in ds.layer
+    if "layer" in ds["top"].dims:
+        raise (ValueError("set_layer_botm does not support top with a layer dimension"))
     lay = np.where(ds.layer == layer)[0][0]
     # if lay > 0 and np.any(botm > ds["botm"][lay - 1]):
     #    raise (Exception("set_layer_botm cannot change botm of higher layers yet"))
@@ -779,6 +785,9 @@ def remove_thin_layers(ds, min_thickness=0.1):
 
     THe thickness of the removed cells is added to the first active layer below
     """
+    if "layer" in ds["top"].dims:
+        msg = "remove_thin_layers does not support top with a layer dimension"
+        raise (ValueError(msg))
     thickness = calculate_thickness(ds)
     for lay_org in range(len(ds.layer)):
         # determine where the layer is too thin
diff --git a/nlmod/dims/time.py b/nlmod/dims/time.py
index 697b5bbe..d74ebcfd 100644
--- a/nlmod/dims/time.py
+++ b/nlmod/dims/time.py
@@ -500,7 +500,9 @@ def dataframe_to_flopy_timeseries(
     filename=None,
     time_series_namerecord=None,
     interpolation_methodrecord="stepwise",
+    append=False,
 ):
+    assert not df.isna().any(axis=None)
     if ds is not None:
         # set index to days after the start of the simulation
         df = df.copy()
@@ -516,7 +518,10 @@ def dataframe_to_flopy_timeseries(
 
     if isinstance(interpolation_methodrecord, str):
         interpolation_methodrecord = [interpolation_methodrecord] * len(df.columns)
-    package.ts.initialize(
+
+    # initialize or append a new package
+    method = package.ts.append_package if append else package.ts.initialize
+    method(
         filename=filename,
         timeseries=timeseries,
         time_series_namerecord=time_series_namerecord,
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index 92f530e0..a31d5444 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -131,7 +131,7 @@ def _add_ts_to_ds(timeseries, loc_sel, variable, ds):
     """Add a timeseries to a variable at location loc_sel in model DataSet."""
     end = pd.Timestamp(ds.time.data[-1])
     if timeseries.index[-1] < end:
-        raise ValueError(f"no recharge available at {timeseries.name} for date {end}")
+        raise ValueError(f"no data available at {timeseries.name} for date {end}")
 
     # fill recharge data array
     model_recharge = pd.Series(index=ds.time, dtype=float)
@@ -277,11 +277,19 @@ def get_knmi_at_locations(ds, start="2010", end=None, most_common_station=False)
 
     # get knmi data stations closest to any grid cell
     oc_knmi_prec = hpd.ObsCollection.from_knmi(
-        stns=stns_rd, starts=[start], ends=[end], meteo_vars=["RD"]
+        stns=stns_rd,
+        starts=[start],
+        ends=[end],
+        meteo_vars=["RD"],
+        fill_missing_obs=True,
     )
 
     oc_knmi_evap = hpd.ObsCollection.from_knmi(
-        stns=stns_ev24, starts=[start], ends=[end], meteo_vars=["EV24"]
+        stns=stns_ev24,
+        starts=[start],
+        ends=[end],
+        meteo_vars=["EV24"],
+        fill_missing_obs=True,
     )
 
     return locations, oc_knmi_prec, oc_knmi_evap

From b98c2daba430c0d3acd5285b766bf3b56e3c2d04 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Wed, 18 Oct 2023 12:30:03 +0200
Subject: [PATCH 04/12] Add some options to modpath

---
 nlmod/modpath/modpath.py | 16 ++++++++++++----
 1 file changed, 12 insertions(+), 4 deletions(-)

diff --git a/nlmod/modpath/modpath.py b/nlmod/modpath/modpath.py
index 4b3a4588..20db899e 100644
--- a/nlmod/modpath/modpath.py
+++ b/nlmod/modpath/modpath.py
@@ -430,7 +430,15 @@ def pg_from_pd(nodes, localx=0.5, localy=0.5, localz=0.5):
 
 
 def sim(
-    mpf, particlegroups, direction="backward", gwf=None, ref_time=None, stoptime=None
+    mpf,
+    particlegroups,
+    direction="backward",
+    gwf=None,
+    ref_time=None,
+    stoptime=None,
+    simulationtype="combined",
+    weaksinkoption="pass_through",
+    weaksourceoption="pass_through",
 ):
     """Create a modpath backward simulation from a particle group.
 
@@ -472,10 +480,10 @@ def sim(
 
     mpsim = flopy.modpath.Modpath7Sim(
         mpf,
-        simulationtype="combined",
+        simulationtype=simulationtype,
         trackingdirection=direction,
-        weaksinkoption="pass_through",
-        weaksourceoption="pass_through",
+        weaksinkoption=weaksinkoption,
+        weaksourceoption=weaksourceoption,
         referencetime=ref_time,
         stoptimeoption=stoptimeoption,
         stoptime=stoptime,

From 338b628e574ed0bdca17c96449800753ad74b1db Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 10:17:10 +0200
Subject: [PATCH 05/12] Fix Issue #261

---
 nlmod/gwf/gwf.py | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/nlmod/gwf/gwf.py b/nlmod/gwf/gwf.py
index 2ce266db..7bcde3b8 100644
--- a/nlmod/gwf/gwf.py
+++ b/nlmod/gwf/gwf.py
@@ -215,10 +215,6 @@ def _disv(ds, model, length_units="METERS", pname="disv", **kwargs):
         xorigin = ds.attrs["xorigin"]
         yorigin = ds.attrs["yorigin"]
         angrot = ds.attrs["angrot"]
-    elif "extent" in ds.attrs.keys():
-        xorigin = ds.attrs["extent"][0]
-        yorigin = ds.attrs["extent"][2]
-        angrot = 0.0
     else:
         xorigin = 0.0
         yorigin = 0.0

From c9a2c4109c5d7b8991cae71295ae386692ec013a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 10:40:01 +0200
Subject: [PATCH 06/12] Fix codacy stuff and improve some testst

---
 nlmod/read/nhi.py               |  8 ++++----
 tests/test_005_external_data.py | 15 +++++++++++++--
 tests/test_021_nhi.py           |  2 ++
 3 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
index ec7e2db5..42a2b8d8 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -8,13 +8,13 @@
 logger = logging.getLogger(__name__)
 
 
-def download_file(url, pathname, filename=None, overwrite=False):
+def download_file(url, pathname, filename=None, overwrite=False, timeout=120):
     if filename is None:
         filename = url.split("/")[-1]
     fname = os.path.join(pathname, filename)
     if overwrite or not os.path.isfile(fname):
         logger.info(f"Downloading {filename}")
-        r = requests.get(url, allow_redirects=True)
+        r = requests.get(url, allow_redirects=True, timeout=timeout)
         with open(fname, "wb") as file:
             file.write(r.content)
     return fname
@@ -68,7 +68,7 @@ def add_buisdrainage(
     if mask_and_scale:
         nodata = np.nan
     else:
-        nodata = buisdrain_d._FillValue
+        nodata = buisdrain_d.attrs["_FillValue"]
     # set buisdrain_d to nodata where it is 0
     mask = buisdrain_d != 0
     buisdrain_d = buisdrain_d.where(mask, nodata).rio.write_crs(buisdrain_d.rio.crs)
@@ -78,7 +78,7 @@ def add_buisdrainage(
     )
     if not mask_and_scale:
         # set nodata values to NaN
-        ds[depth_var] = ds[depth_var].where(ds[depth_var] != buisdrain_d._FillValue)
+        ds[depth_var] = ds[depth_var].where(ds[depth_var] != nodata)
 
     # from cm to m
     ds[depth_var] = ds[depth_var] / 100.0
diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index 8019b418..97307f6e 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -1,6 +1,7 @@
 import test_001_model
-
+import pandas as pd
 import nlmod
+import pytest
 
 
 def test_get_recharge():
@@ -13,7 +14,9 @@ def test_get_recharge():
 
 def test_get_reacharge_most_common():
     # model with sea
-    ds = test_001_model.get_ds_from_cache("basic_sea_model")
+    # ds = test_001_model.get_ds_from_cache("basic_sea_model")
+    ds = nlmod.get_ds([100000, 110000, 420000, 430000])
+    nlmod.time.set_ds_time(ds, start="2000", time=pd.date_range("2010", "2025"))
 
     # add knmi recharge to the model dataset
     ds.update(nlmod.read.knmi.get_recharge(ds, most_common_station=True))
@@ -31,6 +34,14 @@ def test_get_recharge_steady_state():
     ds.update(nlmod.read.knmi.get_recharge(ds))
 
 
+def test_get_recharge_not_available():
+    ds = nlmod.get_ds([100000, 110000, 420000, 430000])
+    time = pd.date_range("2010", pd.Timestamp.now())
+    ds = nlmod.time.set_ds_time(ds, start="2000", time=time)
+    with pytest.raises(ValueError):
+        ds.update(nlmod.read.knmi.get_recharge(ds))
+
+
 def test_ahn_within_extent():
     extent = [95000.0, 105000.0, 494000.0, 500000.0]
     da = nlmod.read.ahn.get_ahn_from_wcs(extent)
diff --git a/tests/test_021_nhi.py b/tests/test_021_nhi.py
index 8991817a..af768339 100644
--- a/tests/test_021_nhi.py
+++ b/tests/test_021_nhi.py
@@ -2,10 +2,12 @@
 import numpy as np
 import tempfile
 import nlmod
+import pytest
 
 tmpdir = tempfile.gettempdir()
 
 
+@pytest.mark.slow
 def test_buidrainage():
     model_ws = os.path.join(tmpdir, "buidrain")
     ds = nlmod.get_ds([110_000, 130_000, 435_000, 445_000], model_ws=model_ws)

From 5dc73139716c5ac4b95c3035ded63aba42f274dd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 10:50:58 +0200
Subject: [PATCH 07/12] Fix the unit of the conductance of buisdrainage

---
 nlmod/read/nhi.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
index 42a2b8d8..b080b962 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -61,6 +61,8 @@ def add_buisdrainage(
     cond = cond.rio.write_crs(buisdrain_c.rio.crs)
     # resample to model grid
     ds[cond_var] = structured_da_to_ds(cond, ds, method=cond_method)
+    # multiply by area to get a conductance
+    ds[cond_var] = ds[cond_var] * ds["area"]
 
     # use mthod = "mode" on depth, to retreive the depth that occurs most in each cell
     mask_and_scale = False

From 8d9ab6f5c376340e2b7a827e94f39d57e6f351d5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 11:46:13 +0200
Subject: [PATCH 08/12] Fix test

---
 tests/test_005_external_data.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index 97307f6e..7bef7adc 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -14,9 +14,8 @@ def test_get_recharge():
 
 def test_get_reacharge_most_common():
     # model with sea
-    # ds = test_001_model.get_ds_from_cache("basic_sea_model")
     ds = nlmod.get_ds([100000, 110000, 420000, 430000])
-    nlmod.time.set_ds_time(ds, start="2000", time=pd.date_range("2010", "2025"))
+    ds = nlmod.time.set_ds_time(ds, start="2021", time=pd.date_range("2022", "2023"))
 
     # add knmi recharge to the model dataset
     ds.update(nlmod.read.knmi.get_recharge(ds, most_common_station=True))

From 3c7609f0c0404b6a4de0601f6578da85f2d3e0c2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Dav=C3=ADd=20Brakenhoff?= <d.brakenhoff@artesia-water.nl>
Date: Mon, 23 Oct 2023 14:32:48 +0200
Subject: [PATCH 09/12] some minor fixes - added one docstring, fixed some
 typos, sort imports

---
 nlmod/dims/layers.py            |  2 +-
 nlmod/gwf/surface_water.py      | 25 +++++++++++++++++--------
 nlmod/read/knmi.py              |  4 +++-
 nlmod/read/nhi.py               | 16 ++++++++++------
 pyproject.toml                  |  2 +-
 tests/test_005_external_data.py |  7 ++++---
 6 files changed, 36 insertions(+), 20 deletions(-)

diff --git a/nlmod/dims/layers.py b/nlmod/dims/layers.py
index b630257d..c3a6e4b5 100644
--- a/nlmod/dims/layers.py
+++ b/nlmod/dims/layers.py
@@ -788,7 +788,7 @@ def set_minimum_layer_thickness(ds, layer, min_thickness, change="botm"):
 def remove_thin_layers(ds, min_thickness=0.1):
     """Remove layers from cells with a thickness less than min_thickness
 
-    THe thickness of the removed cells is added to the first active layer below
+    The thickness of the removed cells is added to the first active layer below
     """
     if "layer" in ds["top"].dims:
         msg = "remove_thin_layers does not support top with a layer dimension"
diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index f0410331..f79ca188 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -1026,6 +1026,21 @@ def add_season_timeseries(
     filename="season.ts",
     seasons=None,
 ):
+    """Add time series indicating which season is active (e.g. summer/winter).
+
+    Parameters
+    ----------
+    ds : xarray.Dataset
+        xarray dataset used for time discretization
+    package : flopy.mf6 package
+        Modflow 6 package to add time series to
+    summer_months : tuple, optional
+        summer months, by default (4, 5, 6, 7, 8, 9)
+    filename : str, optional
+        name of time series file, by default "season.ts"
+    seasons : tuple of str, optional
+        name of the seasons timeseries , by default None
+    """
     if seasons is None:
         seasons = ["winter", "summer"]
     tmin = pd.to_datetime(ds.time.start)
@@ -1056,14 +1071,8 @@ def add_season_timeseries(
 
 
 def rivdata_from_xylist(gwf, xylist, layer, stage, cond, rbot):
-    # TODO: temporary fix until flopy is patched
-    if gwf.modelgrid.grid_type == "structured":
-        gi = flopy.utils.GridIntersect(gwf.modelgrid, rtree=False)
-        cellids = gi.intersect(xylist, shapetype="linestring")["cellids"]
-    else:
-        gi = flopy.utils.GridIntersect(gwf.modelgrid)
-        cellids = gi.intersects(xylist, shapetype="linestring")["cellids"]
-
+    gi = flopy.utils.GridIntersect(gwf.modelgrid, method="vertex")
+    cellids = gi.intersect(xylist, shapetype="linestring")["cellids"]
     riv_data = []
     for cid in cellids:
         if len(cid) == 2:
diff --git a/nlmod/read/knmi.py b/nlmod/read/knmi.py
index a31d5444..9499817f 100644
--- a/nlmod/read/knmi.py
+++ b/nlmod/read/knmi.py
@@ -131,7 +131,9 @@ def _add_ts_to_ds(timeseries, loc_sel, variable, ds):
     """Add a timeseries to a variable at location loc_sel in model DataSet."""
     end = pd.Timestamp(ds.time.data[-1])
     if timeseries.index[-1] < end:
-        raise ValueError(f"no data available at {timeseries.name} for date {end}")
+        raise ValueError(
+            f"no data available for time series'{timeseries.name}' on date {end}"
+        )
 
     # fill recharge data array
     model_recharge = pd.Series(index=ds.time, dtype=float)
diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
index b080b962..e2c812e8 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -1,8 +1,10 @@
-import os
-import requests
 import logging
+import os
+
 import numpy as np
+import requests
 import rioxarray
+
 from ..dims.resample import structured_da_to_ds
 
 logger = logging.getLogger(__name__)
@@ -47,12 +49,13 @@ def add_buisdrainage(
     # download files if needed
     fname_c, fname_d = download_buisdrainage(pathname)
 
-    # make sure crs is set of ds
+    # make sure crs is set on ds
     if ds.rio.crs is None:
         ds = ds.rio.write_crs(28992)
 
-    # use method = "average" on conductance, to account for locations without pipe
-    # drainage, where the conductance is 0
+    # use cond_methd for conductance
+    # (default is "average" to account for locations without pipe drainage, where the 
+    # conductance is 0)
     buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0]
     # calculate a conductance (per m2) from a resistance
     cond = 1 / buisdrain_c
@@ -64,7 +67,8 @@ def add_buisdrainage(
     # multiply by area to get a conductance
     ds[cond_var] = ds[cond_var] * ds["area"]
 
-    # use mthod = "mode" on depth, to retreive the depth that occurs most in each cell
+    # use depth_method to retrieve the depth
+    # (default is "mode" for depth that occurs most in each cell)
     mask_and_scale = False
     buisdrain_d = rioxarray.open_rasterio(fname_d, mask_and_scale=mask_and_scale)[0]
     if mask_and_scale:
diff --git a/pyproject.toml b/pyproject.toml
index e3de3d4a..d63b685c 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -6,7 +6,7 @@ build-backend = "setuptools.build_meta"
 [project]
 name = "nlmod"
 dynamic = ["version"]
-description = "nlmod is an open-source Python package for building Modflow 6 groundwater models from online data sources in The Netherlands"
+description = "Python package to build, run and visualize MODFLOW 6 groundwater models in the Netherlands."
 license = { file = "LICENSE" }
 readme = "README.md"
 authors = [
diff --git a/tests/test_005_external_data.py b/tests/test_005_external_data.py
index 7bef7adc..da72fe7a 100644
--- a/tests/test_005_external_data.py
+++ b/tests/test_005_external_data.py
@@ -1,7 +1,8 @@
-import test_001_model
 import pandas as pd
-import nlmod
 import pytest
+import test_001_model
+
+import nlmod
 
 
 def test_get_recharge():
@@ -12,7 +13,7 @@ def test_get_recharge():
     ds.update(nlmod.read.knmi.get_recharge(ds))
 
 
-def test_get_reacharge_most_common():
+def test_get_recharge_most_common():
     # model with sea
     ds = nlmod.get_ds([100000, 110000, 420000, 430000])
     ds = nlmod.time.set_ds_time(ds, start="2021", time=pd.date_range("2022", "2023"))

From c60c2d657442e1fbfa39227012baafb7de817e61 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 15:01:58 +0200
Subject: [PATCH 10/12] Add docstrings to nhi

---
 nlmod/read/nhi.py | 87 +++++++++++++++++++++++++++++++++++++++++++++--
 1 file changed, 85 insertions(+), 2 deletions(-)

diff --git a/nlmod/read/nhi.py b/nlmod/read/nhi.py
index e2c812e8..858e3a16 100644
--- a/nlmod/read/nhi.py
+++ b/nlmod/read/nhi.py
@@ -10,7 +10,31 @@
 logger = logging.getLogger(__name__)
 
 
-def download_file(url, pathname, filename=None, overwrite=False, timeout=120):
+def download_file(url, pathname, filename=None, overwrite=False, timeout=120.0):
+    """
+    Download a file from the NHI website.
+
+    Parameters
+    ----------
+    url : str
+        The url of the file to download.
+    pathname : str
+        The pathname to which the file is downloaded.
+    filename : str, optional
+        The name of the file to contain the downloadded data. When filename is None, it
+        is derived from url. The default is None.
+    overwrite : bool, optional
+        Overwrite the file if it allready exists. The default is False.
+    timeout : float, optional
+        How many seconds to wait for the server to send data before giving up. The
+        default is 120.
+
+    Returns
+    -------
+    fname : str
+        The full path of the downloaded file.
+
+    """
     if filename is None:
         filename = url.split("/")[-1]
     fname = os.path.join(pathname, filename)
@@ -23,6 +47,24 @@ def download_file(url, pathname, filename=None, overwrite=False, timeout=120):
 
 
 def download_buisdrainage(pathname, overwrite=False):
+    """
+    Download resistance and depth of buisdrainage from the NHI website
+
+    Parameters
+    ----------
+    pathname : str
+        The pathname to which the files are downloaded.
+    overwrite : bool, optional
+        Overwrite the files if they allready exists. The default is False.
+
+    Returns
+    -------
+    fname_c : str
+        The full path of the downloaded file containing the resistance of buisdrainage.
+    fname_d : str
+        The full path of the downloaded file containing the depth of buisdrainage.
+
+    """
     url_bas = "https://thredds.data.nhi.nu/thredds/fileServer/opendap/models/nhi3_2/25m"
 
     # download resistance
@@ -44,6 +86,47 @@ def add_buisdrainage(
     cond_method="average",
     depth_method="mode",
 ):
+    """
+    Add data about the buisdrainage to the model Dataset.
+
+    This data consists of the conductance of buisdrainage (m2/d) and the depth of
+    buisdrainage (m to surface level). With the default settings for `cond_method` and
+    `depth_method`, the conductance is calculated from the area weighted average of the
+    resistance in a cell, while the depth is set to the value which appears most often
+    in a cell. Cells without data (or 0 in the case of the resistance) are ignored in
+    these calculations.
+
+    The original data of the resistance and the depth of buisdrainage at a 25x25 m scale
+    is downloaded to `pathname` when not found.
+
+    Parameters
+    ----------
+    ds : xr.Dataset
+        The model Dataset.
+    pathname : str, optional
+        The pathname containing the downloaded files or the pathname to which the files
+        are downloaded. When pathname is None, it is set the the cachedir. The default
+        is None.
+    cond_var : str, optional
+        The name of the variable in ds to contain the data about the conductance of
+        buisdrainage. The default is "buisdrain_cond".
+    depth_var : str, optional
+        The name of the variable in ds to contain the data about the depth of
+        buisdrainage. The default is "buisdrain_depth".
+    cond_method : str, optional
+        The method to transform the conductance of buisdrainage to the model Dataset.
+        The default is "average".
+    depth_method : str, optional
+        The method to transform the depth of buisdrainage to the model Dataset. The
+        default is "mode".
+
+    Returns
+    -------
+    ds : xr.Dataset
+        The model dataset with added variables with the names `cond_var` and
+        `depth_var`.
+
+    """
     if pathname is None:
         pathname = ds.cachedir
     # download files if needed
@@ -54,7 +137,7 @@ def add_buisdrainage(
         ds = ds.rio.write_crs(28992)
 
     # use cond_methd for conductance
-    # (default is "average" to account for locations without pipe drainage, where the 
+    # (default is "average" to account for locations without pipe drainage, where the
     # conductance is 0)
     buisdrain_c = rioxarray.open_rasterio(fname_c, mask_and_scale=True)[0]
     # calculate a conductance (per m2) from a resistance

From 8fc0101193c74da34044ad01e82591d7b5edb60a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 20:57:31 +0200
Subject: [PATCH 11/12] Replace seasons by winter_name and summer_name in
 surface_water.py

---
 nlmod/gwf/surface_water.py | 25 ++++++++++++++++---------
 1 file changed, 16 insertions(+), 9 deletions(-)

diff --git a/nlmod/gwf/surface_water.py b/nlmod/gwf/surface_water.py
index f79ca188..aed11a13 100644
--- a/nlmod/gwf/surface_water.py
+++ b/nlmod/gwf/surface_water.py
@@ -1015,7 +1015,13 @@ def gdf_to_seasonal_pkg(
         **kwargs,
     )
     # add timeseries for the seasons 'winter' and 'summer'
-    add_season_timeseries(ds, package, summer_months=summer_months, seasons=seasons)
+    add_season_timeseries(
+        ds,
+        package,
+        summer_months=summer_months,
+        winter_name="winter",
+        summer_name="summer",
+    )
     return package
 
 
@@ -1024,7 +1030,8 @@ def add_season_timeseries(
     package,
     summer_months=(4, 5, 6, 7, 8, 9),
     filename="season.ts",
-    seasons=None,
+    winter_name="winter",
+    summer_name="summer",
 ):
     """Add time series indicating which season is active (e.g. summer/winter).
 
@@ -1035,14 +1042,14 @@ def add_season_timeseries(
     package : flopy.mf6 package
         Modflow 6 package to add time series to
     summer_months : tuple, optional
-        summer months, by default (4, 5, 6, 7, 8, 9)
+        summer months. The default is (4, 5, 6, 7, 8, 9), so from april to september.
     filename : str, optional
-        name of time series file, by default "season.ts"
-    seasons : tuple of str, optional
-        name of the seasons timeseries , by default None
+        name of time series file. The default is "season.ts".
+    winter_name : str, optional
+        The name of the time-series with ones in winter. The default is "winter".
+    summer_name : str, optional
+        The name of the time-series with ones in summer. The default is "summer".
     """
-    if seasons is None:
-        seasons = ["winter", "summer"]
     tmin = pd.to_datetime(ds.time.start)
     if tmin.month in summer_months:
         ts_data = [(0.0, 0.0, 1.0)]
@@ -1065,7 +1072,7 @@ def add_season_timeseries(
     package.ts.initialize(
         filename=filename,
         timeseries=ts_data,
-        time_series_namerecord=seasons,
+        time_series_namerecord=[winter_name, summer_name],
         interpolation_methodrecord=["stepwise", "stepwise"],
     )
 

From 1f7d70ab996653816dc61ce823f260abc52565f3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ruben=20Calj=C3=A9?= <r.calje@artesia-water.nl>
Date: Mon, 23 Oct 2023 21:49:40 +0200
Subject: [PATCH 12/12] Remove warnings and and output from unsaturated
 notebooks

---
 docs/examples/17_unsaturated_zone_flow.ipynb | 16 ++++++++++++++--
 1 file changed, 14 insertions(+), 2 deletions(-)

diff --git a/docs/examples/17_unsaturated_zone_flow.ipynb b/docs/examples/17_unsaturated_zone_flow.ipynb
index 9dbd49eb..68c675ad 100644
--- a/docs/examples/17_unsaturated_zone_flow.ipynb
+++ b/docs/examples/17_unsaturated_zone_flow.ipynb
@@ -38,6 +38,18 @@
     "nlmod.show_versions()"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "34f28281-fa68-4618-b063-8b9604306974",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Ignore warning about 1d layers, in numpy 1.26.1 and flopy 3.4.3\n",
+    "import warnings\n",
+    "warnings.filterwarnings(\"ignore\", message=\"Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future.\")"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "9813a7da",
@@ -241,7 +253,7 @@
    "outputs": [],
    "source": [
     "# run model\n",
-    "nlmod.sim.write_and_run(sim, ds)\n",
+    "nlmod.sim.write_and_run(sim, ds, silent=True)\n",
     "\n",
     "# get heads\n",
     "head_rch = nlmod.gwf.get_heads_da(ds)\n",
@@ -311,7 +323,7 @@
    "outputs": [],
    "source": [
     "# run model\n",
-    "nlmod.sim.write_and_run(sim, ds)\n",
+    "nlmod.sim.write_and_run(sim, ds, silent=True)\n",
     "\n",
     "# get heads\n",
     "head_rch = nlmod.gwf.get_heads_da(ds)\n",