Download bofek

.gitignore

@@ -145,7 +145,7 @@ cython_debug/
 nlmod/bin/*
 !nlmod/bin/
 !nlmod/bin/mp7_2_002_provisional
-flowchartnlmod.pptx
+nlmod/data/bofek/*
 
 tests/data/*
 !tests/data/**/

nlmod/gwf/recharge.py

@@ -1,4 +1,5 @@
 import logging
+import numbers
 
 import flopy
 import numpy as np
@@ -370,7 +371,7 @@ def ds_to_uzf(
 
     # generate packagedata
     surfdep = _get_value_from_ds_datavar(ds, "surfdep", surfdep, return_da=True)
-    vks = _get_value_from_ds_datavar(ds, "vk", vks, return_da=True)
+    vks = _get_value_from_ds_datavar(ds, "vks", vks, return_da=True)
     thtr = _get_value_from_ds_datavar(ds, "thtr", thtr, return_da=True)
     thts = _get_value_from_ds_datavar(ds, "thts", thts, return_da=True)
     thti = _get_value_from_ds_datavar(ds, "thti", thti, return_da=True)
@@ -496,7 +497,9 @@ def ds_to_uzf(
         # account for surfdep, as this decreases the height of the top of the upper cell
         # otherwise modflow may return an error
         thickness = calculate_thickness(ds)
-        thickness = [thickness[x] - landflag[x] * surfdep / 2 for x in cellids_obs]
+        if isinstance(surfdep, numbers.Number):
+            surfdep = xr.ones_like(thickness) * surfdep
+        thickness = [thickness[x] - landflag[x] * surfdep[x] / 2 for x in cellids_obs]
 
         # create observations list
         obsdepths = []

nlmod/read/__init__.py

@@ -3,6 +3,7 @@
     administrative,
     ahn,
     bgt,
+    bofek,
     bro,
     brp,
     geotop,

nlmod/read/ahn.py

@@ -32,8 +32,8 @@ def get_ahn(ds=None, identifier="AHN4_5M_M", method="average", extent=None, **kw
 
     Parameters
     ----------
-    ds : xr.Dataset
-        dataset with the model information.
+    ds : xr.Dataset, optional
+        dataset with the model information. If None the extent is used.
     identifier : str, optional
         Possible values for the different AHN-versions are (casing is important):
             AHN1: 'AHN1_5M'
@@ -49,7 +49,7 @@ def get_ahn(ds=None, identifier="AHN4_5M_M", method="average", extent=None, **kw
         nlmod.resample.structured_da_to_ds for possible values. The default is
         'average'.
     extent : list, tuple or np.array, optional
-        extent. The default is None.
+        extent xmin, xmax, ymin, ymax. Only used if ds is None. The default is None.
 
     Returns
     -------

nlmod/read/bofek.py

@@ -0,0 +1,83 @@
+import io
+import logging
+import shutil
+import zipfile
+from pathlib import Path
+
+import geopandas as gpd
+import requests
+
+from nlmod import NLMOD_DATADIR, cache, dims, util
+
+logger = logging.getLogger(__name__)
+
+
+@cache.cache_pickle
+def get_gdf_bofek(ds=None, extent=None):
+    """Get geodataframe of bofek 2020 wihtin the extent of the model.
+
+    It does so by downloading a zip file (> 100 MB) and extracting the relevant
+    geodatabase. Therefore the function can be slow, ~35 seconds depending on your
+    internet connection.
+
+    Parameters
+    ----------
+    ds : xr.DataSet, optional
+        dataset containing relevant model information. The default is None.
+    extent : list, tuple or np.array, optional
+        extent xmin, xmax, ymin, ymax. Only used if ds is None. The default is None.
+
+    Returns
+    -------
+    gdf_bofek : GeoDataframe
+        Bofek2020 geodataframe with a column 'BOFEK2020' containing the bofek cluster
+        codes
+    """
+    import py7zr
+
+    if extent is None and ds is not None:
+        extent = dims.get_extent(ds)
+
+    # set paths
+    tmpdir = Path(NLMOD_DATADIR)
+    fname_7z = tmpdir / 'BOFEK2020_GIS.7z'
+    fname_bofek = tmpdir / 'GIS' / 'BOFEK2020_bestanden' / 'BOFEK2020.gdb'
+    fname_bofek_geojson = tmpdir / 'bofek' / 'BOFEK2020.geojson'
+    bofek_zip_url = 'https://www.wur.nl/nl/show/bofek-2020-gis-1.htm'
+
+    if not fname_bofek_geojson.exists():
+        # download zip
+        logger.info('Downloading BOFEK2020 GIS data (~35 seconds)')
+        r = requests.get(bofek_zip_url, timeout=3600)
+        z = zipfile.ZipFile(io.BytesIO(r.content))
+
+        # extract 7z
+        logger.debug('Extracting zipped BOFEK2020 GIS data')
+        z.extractall(tmpdir)
+        with py7zr.SevenZipFile(fname_7z, mode='r') as z:
+            z.extract(targets=['GIS/BOFEK2020_bestanden/BOFEK2020.gdb'],
+                      path=tmpdir,
+                      recursive=True)
+
+        # clean up
+        logger.debug('Remove zip files')
+        Path(fname_7z).unlink()
+
+        # read geodatabase
+        logger.debug('convert geodatabase to geojson')
+        gdf_bofek = gpd.read_file(fname_bofek)
+
+        # save to geojson
+        gdf_bofek.to_file(fname_bofek_geojson, driver='GeoJSON')
+
+        # remove geodatabase
+        shutil.rmtree(fname_bofek)
+
+    # read geojson
+    logger.debug(f'read bofek2020 geojson from {fname_bofek_geojson}')
+    gdf_bofek = gpd.read_file(fname_bofek_geojson)
+
+    # slice to extent
+    gdf_bofek = util.gdf_within_extent(gdf_bofek, extent)
+
+    return gdf_bofek

nlmod/util.py

@@ -719,6 +719,8 @@ def gdf_within_extent(gdf, extent):
             raise TypeError(f"Only accepts single geometry type not {geom_types}")
     elif geom_types[0] == "Polygon":
         gdf = gpd.overlay(gdf, gdf_extent)
+    elif geom_types[0] == "MultiPolygon":
+        gdf = gpd.overlay(gdf, gdf_extent)
     elif geom_types[0] == "LineString":
         gdf = gpd.sjoin(gdf, gdf_extent)
     elif geom_types[0] == "Point":

pyproject.toml

@@ -57,7 +57,15 @@ repository = "https://github.com/gwmod/nlmod"
 documentation = "https://nlmod.readthedocs.io/en/latest/"
 
 [project.optional-dependencies]
-full = ["nlmod[knmi]", "gdown", "geocube", "rasterstats", "contextily", "scikit-image"]
+full = [
+    "nlmod[knmi]",
+    "gdown",
+    "geocube",
+    "rasterstats",
+    "contextily",
+    "scikit-image",
+    "py7zr"
+]
 knmi = ["h5netcdf", "nlmod[grib]"]
 grib = ["cfgrib", "ecmwflibs"]
 test = ["pytest>=7", "pytest-cov", "pytest-dependency"]

tests/test_005_external_data.py

@@ -131,3 +131,14 @@ def test_get_surface_water_ghb():
 def test_get_brp():
     extent = [116500, 120000, 439000, 442000]
     nlmod.read.brp.get_percelen(extent)
+
+
+# disable because slow (~35 seconds depending on internet connection)
+# def test_get_bofek():
+#     # model with sea
+#     ds = test_001_model.get_ds_from_cache("basic_sea_model")
+
+#     # add knmi recharge to the model dataset
+#     gdf_bofek = nlmod.read.bofek.get_gdf_bofek(ds)
+
+#     assert not gdf_bofek.empty, "Bofek geodataframe is empty"