Fix crs of buidrain data

.github/workflows/ci.yml

@@ -17,7 +17,7 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: [3.11]
+        python-version: [3.11, 3.12]
 
     steps:
       - uses: actions/checkout@v7

.readthedocs.yml

@@ -9,7 +9,7 @@ version: 2
 build:
   os: ubuntu-lts-latest
   tools:
-    python: "3.11"
+    python: "3.12"
 
 # Build documentation in the docs/ directory with Sphinx
 sphinx:
@@ -26,4 +26,3 @@ python:
       extra_requirements:
         - nbtest
         - rtd
-

docs/advanced_stress_packages/01_lake.ipynb

@@ -426,7 +426,7 @@
    "outputs": [],
    "source": [
     "lak_stage = pd.read_csv(os.path.join(ds.model_ws, \"lak_STAGE.csv\"), index_col=0)\n",
-    "lak_stage.index = pd.to_datetime(ds.time.start) + pd.to_timedelta(lak_stage.index, \"d\")\n",
+    "lak_stage.index = pd.to_datetime(ds.time.start) + pd.to_timedelta(lak_stage.index, \"D\")\n",
     "lak_stage.columns = [x.capitalize() for x in lak_stage.columns]\n",
     "\n",
     "f, ax = plt.subplots(figsize=(10, 6), layout=\"constrained\")\n",
@@ -533,7 +533,7 @@
     "    for key in lake_bgt.keys():\n",
     "        data[key] = [x[x[\"node\"] == lakeno + 1][\"q\"].sum() for x in lake_bgt[key]]\n",
     "\n",
-    "    index = pd.to_datetime(ds.time.start) + pd.to_timedelta(cbf_lak.get_times(), \"d\")\n",
+    "    index = pd.to_datetime(ds.time.start) + pd.to_timedelta(cbf_lak.get_times(), \"D\")\n",
     "    df = pd.DataFrame(data, index=index)\n",
     "\n",
     "    f, ax = plt.subplots(figsize=(10, 6), layout=\"constrained\")\n",
@@ -548,7 +548,7 @@
     "    outflow = df_block.where(df_block < 0, 0.0)\n",
     "    outflow = outflow.loc[:, ~(outflow == 0).all()]\n",
     "\n",
-    "    ymax = max(inflow.sum(1).max(), outflow.sum(1).max())\n",
+    "    ymax = max(inflow.sum(axis=1).max(), outflow.sum(axis=1).max())\n",
     "\n",
     "    inflow.plot.area(color=colors, ax=ax, linewidth=0)\n",
     "    handles_in, labels_in = ax.get_legend_handles_labels()\n",

docs/examples/09_schoonhoven.ipynb

@@ -854,7 +854,7 @@
     "# add grid\n",
     "dcs.plot_grid()\n",
     "# add labels with layer names\n",
-    "dcs.plot_layers(alpha=0.0, min_label_area=1000)"
+    "dcs.plot_layers(alpha=0.0, min_label_area=1000);"
    ]
   }
  ],

docs/examples/16_groundwater_transport.ipynb

@@ -72,7 +72,7 @@
     "start_time = \"2010-1-1\"\n",
     "starting_head = 1.0\n",
     "\n",
-    "municipalities = nlmod.read.administrative.download_municipalities(\n",
+    "municipalities = nlmod.read.administrative.download_municipalities_gdf(\n",
     "    extent=extent_hbossche\n",
     ")"
    ]
@@ -90,7 +90,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "fig, ax = nlmod.plot.get_map(extent_hbossche, background=\"OpenStreetMap.Mapnik\")"
+    "fig, ax = nlmod.plot.get_map(extent_hbossche, background=True)"
    ]
   },
   {
@@ -254,7 +254,9 @@
    "outputs": [],
    "source": [
     "# aggregate chloride to our layer model using weighted mean\n",
-    "cli_da = nlmod.layers.aggregate_by_weighted_mean_to_ds(ds, cli, \"3d-chloride\")\n",
+    "cli_da = nlmod.layers.aggregate_by_weighted_mean_to_ds(\n",
+    "    ds, cli, \"3d-chloride\", source_botm_name=\"bottom\"\n",
+    ")\n",
     "\n",
     "# interpolate NaNs nearest\n",
     "for ilay in range(cli_da.shape[0]):\n",
@@ -496,7 +498,7 @@
    "outputs": [],
    "source": [
     "# plot using flopy\n",
-    "fig, ax = nlmod.plot.get_map(extent_hbossche, background=\"OpenStreetMap.Mapnik\")\n",
+    "fig, ax = nlmod.plot.get_map(extent_hbossche, background=True)\n",
     "pmv = fp.plot.PlotMapView(model=gwf, layer=0, ax=ax)\n",
     "# pc = pmv.plot_array(c.isel(time=0), cmap=\"Spectral_r\")\n",
     "pmv.plot_bc(\"GHB\", plotAll=True, alpha=0.1, label=\"GHB\")\n",

docs/examples/nested_models.ipynb

@@ -31,6 +31,16 @@
     "import nlmod"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nlmod.util.get_color_logger(\"INFO\")\n",
+    "nlmod.show_versions()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},

nlmod/dims/base.py

@@ -688,7 +688,7 @@ def get_ds(
     coords = {"x": x, "y": y, "layer": layer}
     if angrot != 0.0:
         affine = grid.get_affine_mod_to_world(attrs)
-        xc, yc = affine * np.meshgrid(x, y)
+        xc, yc = affine @ np.meshgrid(x, y)
         coords["xc"] = (("y", "x"), xc)
         coords["yc"] = (("y", "x"), yc)
 

nlmod/dims/grid.py

@@ -832,7 +832,7 @@ def ds_to_gridprops(ds_in, gridprops, method="nearest", icvert_nodata=-1):
             )
     if is_rotated(ds_out):
         affine = get_affine_mod_to_world(ds_out)
-        ds_out["xc"], ds_out["yc"] = affine * (ds_out.x, ds_out.y)
+        ds_out["xc"], ds_out["yc"] = affine @ (ds_out.x, ds_out.y)
 
     if "area" in gridprops:
         if "area" in ds_out:
@@ -2553,7 +2553,7 @@ def get_affine_mod_to_world(ds):
     xorigin = attrs["xorigin"]
     yorigin = attrs["yorigin"]
     angrot = attrs["angrot"]
-    return Affine.translation(xorigin, yorigin) * Affine.rotation(angrot)
+    return Affine.translation(xorigin, yorigin) @ Affine.rotation(angrot)
 
 
 def get_affine_world_to_mod(ds):
@@ -2562,7 +2562,7 @@ def get_affine_world_to_mod(ds):
     xorigin = attrs["xorigin"]
     yorigin = attrs["yorigin"]
     angrot = attrs["angrot"]
-    return Affine.rotation(-angrot) * Affine.translation(-xorigin, -yorigin)
+    return Affine.rotation(-angrot) @ Affine.translation(-xorigin, -yorigin)
 
 
 def get_affine(ds, sx=None, sy=None):
@@ -2588,10 +2588,10 @@ def get_affine(ds, sx=None, sy=None):
         yoff = yorigin + dy * np.cos(angrot * np.pi / 180)
         return (
             Affine.translation(xoff, yoff)
-            * Affine.scale(sx, sy)
-            * Affine.rotation(angrot)
+            @ Affine.scale(sx, sy)
+            @ Affine.rotation(angrot)
         )
     else:
         xoff = attrs["extent"][0]
         yoff = attrs["extent"][3]
-        return Affine.translation(xoff, yoff) * Affine.scale(sx, sy)
+        return Affine.translation(xoff, yoff) @ Affine.scale(sx, sy)

nlmod/observations.py

@@ -6,10 +6,11 @@
 import xarray as xr
 from geopandas import GeoDataFrame, points_from_xy
 from scipy.spatial import Delaunay
-from tqdm import tqdm
 
 from nlmod.dims import grid
 
+from .util import tqdm
+
 logger = logging.getLogger(__name__)
 
 

nlmod/read/nhi.py

@@ -86,6 +86,7 @@ def add_buisdrainage(
     depth_var="buisdrain_depth",
     cond_method="average",
     depth_method="mode",
+    fix_crs=None,
 ):
     """Add data about the buisdrainage to the model Dataset.
 
@@ -114,6 +115,11 @@ def add_buisdrainage(
     depth_method : str, optional
         The method to transform the depth of buisdrainage to the model Dataset. The
         default is "mode".
+    fix_crs : bool, optional
+        Fix the coordinate shift (~100 m to the northeast) in the original NHI
+        buisdrainage dataset caused by a missing datum transformation (Bessel 1841 to
+        WGS84) in the source dataset. If None, a warning is issued and fix_crs
+        defaults to False to preserve backwards compatibility. The default is None.
 
     Returns
     -------
@@ -130,11 +136,12 @@ def add_buisdrainage(
     )
     ds_out = discretize_buisdrainage(
         ds,
-        pathname,
-        cond_var,
-        depth_var,
-        cond_method,
-        depth_method,
+        pathname=pathname,
+        cond_var=cond_var,
+        depth_var=depth_var,
+        cond_method=cond_method,
+        depth_method=depth_method,
+        fix_crs=fix_crs,
     )
     ds.update(ds_out)
     return ds
@@ -147,6 +154,7 @@ def discretize_buisdrainage(
     depth_var="buisdrain_depth",
     cond_method="average",
     depth_method="mode",
+    fix_crs=None,
 ):
     """Add data about the buisdrainage to the model Dataset.
 
@@ -180,13 +188,27 @@ def discretize_buisdrainage(
     depth_method : str, optional
         The method to transform the depth of buisdrainage to the model Dataset. The
         default is "mode".
+    fix_crs : bool, optional
+        Fix the coordinate shift (~100 m to the northeast) in the original NHI
+        buisdrainage dataset caused by a missing datum transformation (Bessel 1841 to
+        WGS84) in the source dataset. If None, a warning is issued and fix_crs
+        defaults to False to preserve backwards compatibility. The default is None.
 
     Returns
     -------
     ds : xr.Dataset
         The model dataset with added variables with the names `cond_var` and
         `depth_var`.
     """
+    if fix_crs is None:
+        logger.warning(
+            "`fix_crs` was not specified and defaults to False. Set `fix_crs=True` to "
+            "fix the ~100 m coordinate shift (to the northeast) in the original NHI "
+            "buisdrainage dataset caused by a missing datum transformation in the "
+            "source dataset."
+        )
+        fix_crs = False
+
     if pathname is None:
         pathname = ds.cachedir
     # download files if needed
@@ -198,10 +220,17 @@ def discretize_buisdrainage(
 
     ds_out = get_ds_empty(ds, keep_coords=("y", "x"))
 
+    src_crs = (
+        "+proj=longlat +ellps=bessel "
+        "+towgs84=565.417,50.3319,465.552,-0.398957,0.343988,-1.8774,4.0725 +no_defs"
+    )
+
     # 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]
+    if fix_crs:
+        buisdrain_c = buisdrain_c.rio.write_crs(src_crs)
     # calculate a conductance (per m2) from a resistance
     cond = 1 / buisdrain_c
     # set conductance to 0 where resistance is infinite or 0
@@ -216,6 +245,8 @@ def discretize_buisdrainage(
     # (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 fix_crs:
+        buisdrain_d = buisdrain_d.rio.write_crs(src_crs)
     if mask_and_scale:
         nodata = np.nan
     else:

nlmod/read/rws.py

@@ -8,6 +8,7 @@
 import numpy as np
 import xarray as xr
 from rioxarray.merge import merge_arrays
+from rasterio.env import Env
 
 from .. import NLMOD_DATADIR, cache, dims, util
 from ..util import tqdm
@@ -340,9 +341,9 @@ def get_gdr_configuration() -> dict:
     config["bodemhoogte_1m"] = {
         "url": (
             "https://geo.rijkswaterstaat.nl/arcgis/rest/services/GDR/"
-            "bodemhoogte_index/FeatureServer"
+            "bodemhoogte_index/MapServer"
         ),
-        "layer": 1,
+        "layer": 1,  # bodemhoogte_1mtr
     }
     # NOTE: the 20m resolution is no longer available from the GDR service via a
     # geodataframe containing the url.
@@ -529,34 +530,38 @@ def download_bathymetry(
     xmin, xmax, ymin, ymax = extent
     dataarrays = []
 
-    for _, row in tqdm(
-        gdf.iterrows(), desc="Downloading bathymetry", total=gdf.index.size
-    ):
-        url = row["geotiff"]
-        ds = xr.open_dataset(url, engine="rasterio")
-        ds = ds.assign_coords({"y": ds["y"].round(0), "x": ds["x"].round(0)})
-        da = (
-            ds["band_data"]
-            .sel(band=1, x=slice(xmin, xmax), y=slice(ymax, ymin))
-            .drop_vars("band")
-        )
-        if chunks:
-            da = da.chunk(chunks)
-        dataarrays.append(da)
-
-    if len(dataarrays) > 1:
-        da = merge_arrays(
-            dataarrays,
-            bounds=[xmin, ymin, xmax, ymax],
-            res=res,
-            method=method,
-        )
-    else:
-        da = dataarrays[0]
-        if res is not None:
-            da = da.rio.reproject(
-                da.rio.crs,
+    rasterio_env = {
+        "GDAL_DISABLE_READDIR_ON_OPEN": "YES",
+    }
+    with Env(**rasterio_env):
+        for _, row in tqdm(
+            gdf.iterrows(), desc="Downloading bathymetry", total=gdf.index.size
+        ):
+            url = row["geotiff"]
+            ds = xr.open_dataset(url, engine="rasterio")
+            ds = ds.assign_coords({"y": ds["y"].round(0), "x": ds["x"].round(0)})
+            da = (
+                ds["band_data"]
+                .sel(band=1, x=slice(xmin, xmax), y=slice(ymax, ymin))
+                .drop_vars("band")
+            )
+            if chunks:
+                da = da.chunk(chunks)
+            dataarrays.append(da)
+
+        if len(dataarrays) > 1:
+            da = merge_arrays(
+                dataarrays,
+                bounds=[xmin, ymin, xmax, ymax],
                 res=res,
-                resampling=method,
+                method=method,
             )
+        else:
+            da = dataarrays[0]
+            if res is not None:
+                da = da.rio.reproject(
+                    da.rio.crs,
+                    res=res,
+                    resampling=method,
+                )
     return da

pyproject.toml

@@ -41,7 +41,6 @@ classifiers = [
     "Intended Audience :: Other Audience",
     "License :: OSI Approved :: MIT License",
     "Programming Language :: Python :: 3 :: Only",
-    "Programming Language :: Python :: 3.10",
     "Programming Language :: Python :: 3.11",
     "Programming Language :: Python :: 3.12",
     "Programming Language :: Python :: 3.13",
@@ -84,6 +83,8 @@ rtd = [
     "sphinx_rtd_theme",
     "nbconvert",
     "openpyxl",
+    "numpy<2.5.0", # until there is a new version of netCDF4 (now 1.7.4) this solves "Setting the shape on a NumPy array has been deprecated"
+    "affine<3.0.1", # until there are newer versions of rasterio (now 1.5.1) and rioxarray (now 0.23.0), this solves "Use @ matmul instead of * mul operator for matrix multiplication"
 ]
 lgn = ["xarray-spatial", "exactextract"]