Buisdrainage

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",

nlmod/dims/layers.py

@@ -689,6 +689,8 @@ def set_model_top(ds, top, min_thickness=0.0):
                 "Either add attribute or use nlmod functions to build the dataset."
             )
         )
+    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:
@@ -712,6 +714,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
@@ -735,6 +739,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"))
@@ -782,8 +788,11 @@ 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"
+        raise (ValueError(msg))
     thickness = calculate_thickness(ds)
     for lay_org in range(len(ds.layer)):
         # determine where the layer is too thin

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
     -------

nlmod/dims/time.py

@@ -503,7 +503,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()
@@ -519,7 +521,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,

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

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,8 +1030,26 @@ def add_season_timeseries(
     package,
     summer_months=(4, 5, 6, 7, 8, 9),
     filename="season.ts",
-    seasons=("winter", "summer"),
+    winter_name="winter",
+    summer_name="summer",
 ):
+    """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. The default is (4, 5, 6, 7, 8, 9), so from april to september.
+    filename : str, optional
+        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".
+    """
     tmin = pd.to_datetime(ds.time.start)
     if tmin.month in summer_months:
         ts_data = [(0.0, 0.0, 1.0)]
@@ -1048,20 +1072,14 @@ 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"],
     )
 
 
 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:

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,

nlmod/read/__init__.py

@@ -13,6 +13,7 @@
     rws,
     waterboard,
     webservices,
+    nhi,
 )
 from .geotop import get_geotop
 from .regis import get_regis

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 recharge 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)