Download bofek

.gitignore

@@ -145,7 +145,7 @@ cython_debug/
 nlmod/bin/*
 !nlmod/bin/
 !nlmod/bin/mp7_2_002_provisional
-nlmod/data/GIS
+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/bofek.py

@@ -1,16 +1,23 @@
-import requests
-import zipfile
 import io
-import geopandas as gpd
+import logging
+import shutil
+import zipfile
 from pathlib import Path
+import requests
+import geopandas as gpd
 from .. import NLMOD_DATADIR, cache, dims, util
 
 
+logger = logging.getLogger(__name__)
+
+
 @cache.cache_pickle
-def get_gdf_bofek(ds=None, extent=None, timeout=120):
-    """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.
+def get_gdf_bofek(ds=None, extent=None, timeout=3600):
+    """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
     ----------
@@ -19,32 +26,37 @@ def get_gdf_bofek(ds=None, extent=None, timeout=120):
     extent : list, tuple or np.array, optional
         extent xmin, xmax, ymin, ymax. Only used if ds is None. The default is None.
     timeout : int, optional
-        timeout time of request in seconds. Default is 120.
+        timeout time of request in seconds. Default is 3600.
 
     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.exists():
+    if not fname_bofek_geojson.exists():
         # download zip
+        logger.info("Downloading BOFEK2020 GIS data (~35 seconds)")
         r = requests.get(bofek_zip_url, timeout=timeout)
+
+        logger.debug("Extracting zipped BOFEK2020 GIS data")
         with zipfile.ZipFile(io.BytesIO(r.content)) as z:
             # extract 7z
             z.extractall(tmpdir)
+
         with py7zr.SevenZipFile(fname_7z, mode="r") as z:
             z.extract(
                 targets=["GIS/BOFEK2020_bestanden/BOFEK2020.gdb"],
@@ -53,10 +65,22 @@ def get_gdf_bofek(ds=None, extent=None, timeout=120):
             )
 
         # clean up
+        logger.debug("Remove zip files")
         Path(fname_7z).unlink()
 
-    # read geodatabase
-    gdf_bofek = gpd.read_file(fname_bofek)
+        # 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)
 
     if extent is not None:
         # slice to extent