Download bofek (#392)

* Download bofek2020

* add downloaded bofek data to gitignore

* codacy

* black formatting of bofek.py

* add download progress bar
uncomment test, but skip using pytest because of slowness

---------

Co-authored-by: Martin Vonk <vonk.mart@gmail.com>
Co-authored-by: Ruben Caljé <r.calje@artesia-water.nl>
Co-authored-by: Davíd Brakenhoff <d.brakenhoff@artesia-water.nl>

Author: 4 people

.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/data/bofek/.gitkeep

@@ -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/gwf/recharge.py

@@ -1,16 +1,24 @@
-import requests
+import logging
+import shutil
 import zipfile
-import io
-import geopandas as gpd
 from pathlib import Path
-from .. import NLMOD_DATADIR, cache, dims, util
+
+import geopandas as gpd
+import requests
+from tqdm.auto import tqdm
+
+from nlmod 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 +27,49 @@ 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
-        r = requests.get(bofek_zip_url, timeout=timeout)
-        with zipfile.ZipFile(io.BytesIO(r.content)) as z:
+        logger.info("Downloading BOFEK2020 GIS data (~35 seconds)")
+        r = requests.get(bofek_zip_url, timeout=timeout, stream=True)
+
+        # show download progress
+        total_size = int(r.headers.get("content-length", 0))
+        block_size = 1024
+        with tqdm(
+            total=total_size, unit="B", unit_scale=True, desc="Downloading BOFEK"
+        ) as progress_bar:
+            with open(tmpdir / "bofek.zip", "wb") as file:
+                for data in r.iter_content(block_size):
+                    progress_bar.update(len(data))
+                    file.write(data)
+
+        # unpack zips
+        logger.debug("Extracting zipped BOFEK2020 GIS data")
+        with zipfile.ZipFile(tmpdir / "bofek.zip") 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 +78,23 @@ def get_gdf_bofek(ds=None, extent=None, timeout=120):
             )
 
         # clean up
+        logger.debug("Remove zip files")
+        Path(tmpdir / "bofek.zip").unlink()
         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

nlmod/read/bofek.py

@@ -3,14 +3,12 @@
 import nlmod
 
 
-# @pytest.mark.skip(reason="too slow")
 def test_get_regis(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
     regis_ds = nlmod.read.regis.get_regis(extent)
 
     assert regis_ds.sizes["layer"] == 20
 
 
-# @pytest.mark.skip(reason="too slow")
 def test_get_regis_botm_layer_BEk1(
     extent=[98700.0, 99000.0, 489500.0, 489700.0],
     botm_layer="MSc",
@@ -44,7 +42,6 @@ def test_get_geotop(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
     nlmod.plot.geotop_lithok_on_map(geotop_ds, z=-20.2, ax=ax)
 
 
-# @pytest.mark.skip(reason="too slow")
 def test_get_regis_geotop(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
     regis_geotop_ds = nlmod.read.regis.get_combined_layer_models(
         extent, use_regis=True, use_geotop=True
@@ -53,7 +50,6 @@ def test_get_regis_geotop(extent=[98600.0, 99000.0, 489400.0, 489700.0]):
     assert regis_geotop_ds.sizes["layer"] == 24
 
 
-# @pytest.mark.skip(reason="too slow")
 def test_get_regis_geotop_keep_all_layers(
     extent=[98600.0, 99000.0, 489400.0, 489700.0],
 ):

tests/test_002_regis_geotop.py

@@ -134,11 +134,12 @@ def test_get_brp():
 
 
 # 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")
+@pytest.mark.skip(reason="slow")
+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)
+    # 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"
+    assert not gdf_bofek.empty, "Bofek geodataframe is empty"