{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1c5755f8",
   "metadata": {},
   "source": [
    "# Using information from GeoTOP\n",
    "Most geohydrological models in the Netherlands use the layer model REGIS as the basis for the geohydrological properties of the model. However, REGIS does not contain information for all layers, of which the holocene confining layer (HLc) is the most important one. In this notebook we will show how to use the voxel model GeoTOP to generate geohydrolocal properties for this layer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3e89ac62",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from shapely.geometry import LineString\n",
    "\n",
    "import nlmod"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0de3f47d",
   "metadata": {},
   "outputs": [],
   "source": [
    "nlmod.util.get_color_logger(\"INFO\")\n",
    "nlmod.show_versions()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "603caa14",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define an extent and a line from southwest to northeast through this extent\n",
    "extent = [100000.0, 105000.0, 499800.0, 500000.0]\n",
    "line = LineString([(extent[0], extent[2]), (extent[1], extent[3])])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38282cfc",
   "metadata": {},
   "source": [
    "We define a method called 'plot_kh_kv', which we can call in cells below to plot the horizontal and vertical conductivity in a cross-section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0df7d7a",
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_kh_kv(\n",
    "    ds,\n",
    "    layer=\"layer\",\n",
    "    variables=None,\n",
    "    zmin=-50.25,\n",
    "    min_label_area=None,\n",
    "    cmap=None,\n",
    "    norm=None,\n",
    "):\n",
    "    if variables is None:\n",
    "        variables = [\"kh\", \"kv\"]\n",
    "    if cmap is None:\n",
    "        cmap = plt.get_cmap(\"turbo_r\")\n",
    "    if norm is None:\n",
    "        boundaries = [0.0001, 0.001, 0.01, 0.1, 1, 5, 10, 20, 50, 100]\n",
    "        norm = matplotlib.colors.BoundaryNorm(boundaries, cmap.N)\n",
    "    for var in variables:\n",
    "        f, ax = plt.subplots(figsize=(10, 5))\n",
    "        cs = nlmod.plot.DatasetCrossSection(ds, line, layer=layer, ax=ax, zmin=zmin)\n",
    "        pc = cs.plot_array(ds[var], norm=norm, cmap=cmap)\n",
    "        if min_label_area is not None:\n",
    "            cs.plot_layers(alpha=0.0, min_label_area=min_label_area)\n",
    "            cs.plot_grid(vertical=False)\n",
    "        format = matplotlib.ticker.FuncFormatter(lambda y, _: \"{:g}\".format(y))\n",
    "        nlmod.plot.colorbar_inside(pc, bounds=[0.05, 0.05, 0.02, 0.9], format=format)\n",
    "        nlmod.plot.title_inside(var, ax=ax)\n",
    "        ax.set_xlabel(\"afstand langs doorsnede (m)\")\n",
    "        ax.set_ylabel(\"z (m NAP)\")\n",
    "        f.tight_layout(pad=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefeba72",
   "metadata": {},
   "source": [
    "## Download GeoTOP\n",
    "We download GeoTOP for a certain extent. We get an xarray.Dataset with voxels of 100 * 100 * 0.5 (dx * dy * dz) m, with variables 'strat' and 'lithok'. We also get the probaliblity of the occurence of each lithoclass in variables named 'kans_*' (since we set `probabilities=True`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb437747",
   "metadata": {},
   "outputs": [],
   "source": [
    "gt = nlmod.read.geotop.download_geotop(extent, probabilities=True)\n",
    "gt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ace59ae0",
   "metadata": {},
   "source": [
    "We plot the lithoclass (soil types) and the stratigraphy (geological units) in two cross-sections using the methods `nlmod.plot.geotop_lithok_in_cross_section` & `nlmod.plot.geotop_strat_in_cross_section`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b86b6b8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(10, 5))\n",
    "nlmod.plot.geotop_lithok_in_cross_section(line, gt)\n",
    "f.tight_layout(pad=0.0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a314f836",
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(10, 5))\n",
    "nlmod.plot.geotop_strat_in_cross_section(line, gt, ax=ax)\n",
    "f.tight_layout(pad=0.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5577c352",
   "metadata": {},
   "source": [
    "## Add hydrological properties (kh and kv)\n",
    "GeoTOP does not contain information about geohydroloical propties directly. We need to calculate this information, using the lithoclass, and optionally the stratigraphy (layer unit). We get this information from a DataFrame, whch needs to contain the columns 'lithok' and 'kh' (and optionally 'strat' and 'kv')."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9506f77d",
   "metadata": {},
   "source": [
    "### Based on lithok\n",
    "With `nlmod.read.geotop.get_lithok_props()` we get a default value for each of the 9 lithoclasses (lithoclass 4 is not used). These values are a rough estimate of the hydrologic conductivity. We recommend changing these values based on local conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4429541a",
   "metadata": {},
   "outputs": [],
   "source": [
    "df = nlmod.read.geotop.get_lithok_props()\n",
    "df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "932b8994",
   "metadata": {},
   "source": [
    "The method `nlmod.read.geotop.add_kh_and_kv` takes this DataFrame, applies it to the GeoTOP voxel-dataset `gt`, and adds the variables `kh` and `kv` to `gt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "11c61bec",
   "metadata": {},
   "outputs": [],
   "source": [
    "gt = nlmod.read.geotop.add_kh_and_kv(gt, df)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7f2e2f5c",
   "metadata": {},
   "source": [
    "When we plot kh and kv we see fine sands get a value of 5 m/d (green) and medium fine sands get a value of 10 m/d (light blue). We see the peat (0.001 m/d) and clay (0.01 m/d) layers as zones with lower conductivities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fa2b2a04",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(gt, layer=\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18450570",
   "metadata": {},
   "source": [
    "### Based on lithok and strat\n",
    "We can also load one of the other  DataFrames that are built into nlmod, using the method `nlmod.geotop.get_kh_kv_table()`. Using this method, a table for certain location can be loaded. Right now, the only allowed value for `kind` is 'Brabant', which is the default value, and loads the hydrological properties per lithoclass and stratigraphic unit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d2624408",
   "metadata": {},
   "outputs": [],
   "source": [
    "df_brab = nlmod.read.geotop.get_kh_kv_table(kind=\"Brabant\")\n",
    "df_brab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1b8d252b",
   "metadata": {},
   "source": [
    "We use this table to add a kh and kv value for each voxel, in variables named 'kh' and 'kv'."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8d8ca0fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "gt_brab = nlmod.read.geotop.add_kh_and_kv(gt.copy(), df_brab)\n",
    "gt_brab"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "015399df",
   "metadata": {},
   "source": [
    "We can plot these values along the same line we plotted the lithoclass-values in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0665b527",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(gt_brab, layer=\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d93c5189",
   "metadata": {},
   "source": [
    "## Aggregate voxels to GeoTOP-layers\n",
    "In the previous example we have added geodrological properties to the voxels of GeoTOP and plotted them. The layers of a groundwatermodel generally are thicker than the thickness of the voxels (0.5 m). Therefore we need to aggregate the voxel-data into the layers of the model. We show this process by using the stratigraphy-data of GeoTOP to form a layer model, using the method `nlmod.read.geotop.to_model_layers`.\n",
    "\n",
    "When values for `kh` and `kv` are present in `gt`, this method also calculates the geohydrological properties of the layer model with the method `nlmod.read.geotop.aggregate_to_ds`. The method calculates the combined horizontal transmissivity, and the combined vertical resistance of all (parts of) voxels in a layer, and calculates a `kh` and `kv` value from this transmissivity and resistance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a39f8e71",
   "metadata": {},
   "outputs": [],
   "source": [
    "gtl = nlmod.read.geotop.to_model_layers(gt)\n",
    "gtl"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9917260a",
   "metadata": {},
   "source": [
    "We can plot the kh and kv value for each of the layers with the same method we used for the voxel-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "500ffa1d",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(gtl, min_label_area=1000.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f274bf06",
   "metadata": {},
   "source": [
    "## Aggregate voxels to REGIS-layers\n",
    "We can use any layer model in `nlmod.read.geotop.aggregate_to_ds()`, also one from REGIS. Let's demonstrate this by downloading REGIS-data for the same extent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dd95c4ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "regis = nlmod.read.download_regis(extent)\n",
    "regis"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87d0e5d0",
   "metadata": {},
   "source": [
    "First we plot the original hydrological properties of REGIS. We see that kh is defined for the aquifers (top plot) and kv is defined for the aquitards (bottom plot). Neither kh or kv is defined for the top layer HLc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fede3a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(regis, min_label_area=1000.0, zmin=-100.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "72fe6c86",
   "metadata": {},
   "source": [
    "We will add varibales `kh_gt` and `kv_gt` that contain the kh- and kv-values calculated from GeoTOP. Layers that do not contain any voxel will get a NaN-value for kh and kv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a32ee6f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# make sure there are no NaNs in top and botm of layers\n",
    "regis = nlmod.read.geotop.aggregate_to_ds(gt, regis, kh=\"kh_gt\", kv=\"kv_gt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0a618866",
   "metadata": {},
   "source": [
    "When we plot the kh and kv value, we see all layers above -50 m NAP now contain values, calculated from the GeoTOP-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3a3f1eb",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(regis, min_label_area=1000.0, zmin=-100.0, variables=[\"kh_gt\", \"kv_gt\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "44a537a1",
   "metadata": {},
   "source": [
    "We can plot kh of one of the layers on a map, for both REGIS and GeoTOP. We generally see that conductivities in GeoTOP are a bit lower."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "681c0d58",
   "metadata": {},
   "outputs": [],
   "source": [
    "layer = \"KRz3\"\n",
    "var = \"kh\"\n",
    "norm = matplotlib.colors.Normalize(0.0, 40.0)\n",
    "\n",
    "f, axes = nlmod.plot.get_map(extent, nrows=2)\n",
    "pc = nlmod.plot.data_array(regis[var].loc[layer], ax=axes[0], norm=norm)\n",
    "nlmod.plot.colorbar_inside(pc, bounds=[0.02, 0.05, 0.02, 0.9], ax=axes[0])\n",
    "nlmod.plot.title_inside(\"REGIS\", ax=axes[0])\n",
    "pc = nlmod.plot.data_array(regis[f\"{var}_gt\"].loc[layer], ax=axes[1], norm=norm)\n",
    "nlmod.plot.title_inside(\"GeoTOP\", ax=axes[1])\n",
    "nlmod.plot.colorbar_inside(pc, bounds=[0.02, 0.05, 0.02, 0.9], ax=axes[1]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "be7292b5",
   "metadata": {},
   "source": [
    "## Using stochastic data from GeoTOP\n",
    "In the previous section we used the most likely values from the lithoclass data of GeoTOP. GeoTOP is constructed by generating 100 realisations of this data. Using these realisations a probablity is determined for the occurence in each pixel for each of the lithoclassses. We can also use these probabilities to determine the kh and kv-value of each voxel. We do this by settting the `stochastic` parameter in `nlmod.read.geotop.add_kh_and_kv` to True. The kh and kv values are now calculated by a weighted average of the lithoclass data in each voxel, where the weights are determined by the probablilities. By default an arithmetic weighted mean is used for kh and a harmonic weighted mean for kv, but these methods can be chosen by the user."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d1de603a",
   "metadata": {},
   "outputs": [],
   "source": [
    "gt = nlmod.read.geotop.add_kh_and_kv(gt, df, stochastic=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27411195",
   "metadata": {},
   "source": [
    "We can plot the kh- and kv-values again. Using the stochastic data generally results in smoother values for kh and kv."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "41bc6161",
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_kh_kv(gt, layer=\"z\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "864daed4",
   "metadata": {},
   "source": [
    "## Geulen\n",
    "\n",
    "Most of the stratigraphic units in Geotop are layered on top of each other. Meaning that these stratipgrahic units are always completely on top or completely below the other units. The so-called \"geulen\" are the exception to this, they can occur above, in-between and below other units. This makes it harderd to construct a layered groundwatermodel from the stratigraphic units.\n",
    "\n",
    "Currently there are three options to create a layer model based on the stratigraphic units in geotop:\n",
    "1. `add_to_layer_below`, a geul is added to the layer directly below the geul\n",
    "2. `add_to_layer_above`, a geul is added to the layer directly above the geul\n",
    "3. `split_layers`, the geul is added as seperate layers in the layer model. If a geul occurs both above and below another unit this unit is also split into two layers.\n",
    "\n",
    "Below an example to showcase the differences between these methods. For this we use a different extent with 2 geulen in the stratigraphic geotop data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "af2b02ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# This extent has two Geulen in it\n",
    "extent = (106000.0, 114000.0, 490500.0, 490700.0)\n",
    "line1 = LineString([(extent[0], extent[2]+50), (extent[1], extent[2]+50)])\n",
    "line2 = LineString([(extent[0], extent[2]+150), (extent[1], extent[2]+150)])\n",
    "\n",
    "# download geotop\n",
    "gt = nlmod.read.geotop.download_geotop(extent)\n",
    "\n",
    "# add hydrological properties\n",
    "df = nlmod.read.geotop.get_lithok_props()\n",
    "gt = nlmod.read.geotop.add_kh_and_kv(gt, df)\n",
    "\n",
    "# set the logger to debug to get all the information\n",
    "logger = nlmod.util.get_color_logger(\"DEBUG\")\n",
    "gtl1 = nlmod.read.geotop.to_model_layers(gt, method_geulen='add_to_layer_below')\n",
    "gtl2 = nlmod.read.geotop.to_model_layers(gt, method_geulen='add_to_layer_above')\n",
    "gtl3 = nlmod.read.geotop.to_model_layers(gt, method_geulen='split_layers')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7a5145e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# plot layer models\n",
    "logger = nlmod.util.get_color_logger(\"INFO\")\n",
    "stratdf = nlmod.read.geotop.get_strat_props().set_index('code')\n",
    "f, axes = plt.subplots(figsize=(10, 15), nrows=4, sharex=True)\n",
    "\n",
    "nlmod.plot.geotop_strat_in_cross_section(line1, gt, ax=axes[0], legend=False)\n",
    "axes[0].set_title('original stratigraphy Geotop')\n",
    "\n",
    "cs1 = nlmod.plot.DatasetCrossSection(gtl1, line1, layer='layer', ax=axes[1])\n",
    "out1 = cs1.plot_layers(colors=stratdf, edgecolor='k',  min_label_area=3000)\n",
    "axes[1].set_title('add_to_layer_below')\n",
    "\n",
    "cs2 = nlmod.plot.DatasetCrossSection(gtl2, line1, layer='layer', ax=axes[2])\n",
    "out2 = cs2.plot_layers(colors=stratdf, edgecolor='k',  min_label_area=3000)\n",
    "axes[2].set_title('add_to_layer_above')\n",
    "\n",
    "cs3 = nlmod.plot.DatasetCrossSection(gtl3, line1, layer='layer', ax=axes[3])\n",
    "out3 = cs3.plot_layers(colors=stratdf, edgecolor='k', min_label_area=3000)\n",
    "axes[3].set_title('split_layers');\n",
    "\n",
    "print('units in original gt:', sum(~np.isnan(np.unique(gt['strat']))))\n",
    "print('layers in gtl1:', len(gtl1.layer))\n",
    "print('layers in gtl2:', len(gtl2.layer))\n",
    "print('layers in gtl3:', len(gtl3.layer))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nlmod",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}