{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# BRO GeoTop\n",
    "\n",
    "GeoST provides support to work with models that are provided by the [BRO](https://www.broloket.nl/ondergrondmodellen/kaart) such as GeoTOP. This section will explain the usage of [`GeoTop`](../api_reference/bro_geotop.rst), how to make easy selections and combinations with for example borehole data contained in a [`BoreholeCollection`](../api_reference/borehole_collection.rst).\n",
    "\n",
    "## Reading GeoTop\n",
    "\n",
    "GeoTOP data can be read from a local NetCDF file of downloaded data from the [BRO](https://www.broloket.nl/ondergrondmodellen/kaart) or directly accessed from the [OPeNDAP-server](https://dinodata.nl/opendap/GeoTOP/geotop.nc.html) of TNO Geological survey. See [`GeoTop.from_netcdf`](../api_reference/bro_geotop.rst) and [`GeoTop.from_opendap`](../api_reference/bro_geotop.rst) in the [API reference](../api_reference.md).\n",
    "\n",
    "## Working with GeoTop\n",
    "\n",
    "Let's load some example GeoTOP data for the Utrecht Science park:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import geost\n",
    "\n",
    "geotop = geost.data.geotop_usp()\n",
    "print(geotop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see this prints information on the type of class (i.e. [`GeoTop`](../api_reference/bro_geotop.rst)), the available data variables, xyz-dimensions and resolution. The [`GeoTop`](../api_reference/bro_geotop.rst) class inherits from GeoST [`VoxelModel`](../api_reference/voxelmodel.rst) and as such, inherits all the attributes and associated selection methods. For example, we can check the horizontal and vertical bounds:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(geotop.horizontal_bounds)\n",
    "print(geotop.vertical_bounds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This prints the \"xmin, ymin, xmax, ymax\" bounding box and the \"zmin, zmax\" of the data. [`GeoTop`](../api_reference/bro_geotop.rst) depends strongly on [xarray.Dataset](https://docs.xarray.dev/en/stable/generated/xarray.Dataset.html) and therefore provides familiar selection methods `.sel` and `.isel` to select specific coordinates or indices, or slices of coordinates or indices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sel = geotop.sel(x=[139_650, 139_750]) # Select specific x-coordinates\n",
    "sel = geotop.isel(x=[2, 5]) # Select specific x-indices\n",
    "\n",
    "sel = geotop.sel(x=slice(140_000, 140_500)) # Select a slice of x-coordinates\n",
    "sel = geotop.isel(x=slice(2, 5)) # Select a slice of x-indices\n",
    "\n",
    "sel = geotop.sel(x=[140_013.23, 140_333.14], method='nearest') # Select the nearest x-coordinates\n",
    "print(sel)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, a selection returns a new [`GeoTop`](../api_reference/bro_geotop.rst) instance for the desired coordinates. Note that the x-resolution has also been changed because the coordinates in the last selection are approximately 300 meters apart.\n",
    "\n",
    "### Selections\n",
    "\n",
    "Often it can be necessary to select GeoTop at specific locations, such as borehole locations, to compare the borehole data with GeoTop. We can easily select GeoTop at the locations of boreholes using the [`GeoTop.select_with_points`](../api_reference/bro_geotop.rst) method. Let's load a [`BoreholeCollection`](../api_reference/borehole_collection.rst) with borehole data for the Utrecht Science Park to see what happens: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "boreholes = geost.data.boreholes_usp() # BoreholeCollection\n",
    "\n",
    "geotop_at_locations = geotop.select_with_points(boreholes.header.gdf)\n",
    "print(geotop_at_locations)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an `xarray.Dataset` with a new \"idx\" dimension and note that variables like \"strat\" now have \"(idx, z)\" dimensions. The `.select_with_points` method selects the voxel columns each borehole is in, in the index order of the boreholes in the [`BoreholeCollection`](../api_reference/borehole_collection.rst). It does not return a `VoxelModel` instance because a `VoxelModel` is only valid if it contains \"x\", \"y\" and \"z\" dimensions.\n",
    "\n",
    "With the selection result, it is possible to compare a borehole with a specific voxel column. Let's for example compare the first borehole in the collection with the voxel column. First, let's print the \"header\" attribute of the borehole collection to find the \"nr\" of the first borehole:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(boreholes.header)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the first number (at index 0) is \"B31H0541\". Notice in the \"surface\" and \"end\" columns that the borehole is between +1.2 m NAP and -9.9 m NAP. We can select the first index from `geotop_at_locations`, select the depths between 1.2 and -9.9 m NAP and create a DataFrame to compare the data. In similar way, selection method for other use cases are available when working with [`GeoTop`](../api_reference/bro_geotop.rst). See the [API reference](../api_reference/bro_geotop.rst) for available methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Increase zmin and zmax by a little to make sure the complete depth of the borehole is covered\n",
    "voxel_column = geotop_at_locations.sel(idx=0, z=slice(-10.5, 1.5))\n",
    "voxel_column = voxel_column.to_dataframe().sort_index(ascending=False) # Make depth increase from top to bottom\n",
    "voxel_column # Check out the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Relabelling GeoTOP information\n",
    "\n",
    "Notice that in the `voxel_column` result, the columns \"strat\" (i.e. stratigraphy) and \"lithok\" (i.e. lithology) contain numbers, instead of names which makes it difficult to understand what the unit and lithology at specific depths are. This makes a comparison with the borehole data difficult.\n",
    "\n",
    "GeoST provides the [`StratGeotop`](../api_reference/geotop_selection.rst) class and [`Lithology`](../api_reference/geotop_selection.rst) class which make it easy to translate these numbers into more meaningful names. Let's import them and check which stratigraphic unit belongs to the number 2010 in the \"strat\" column and what the number 6 in the \"lithok\" means:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from geost.bro import Lithology, StratGeotop\n",
    "\n",
    "unit = StratGeotop.select_values(2010)\n",
    "print(unit)\n",
    "\n",
    "lith = Lithology.select_values(6)\n",
    "print(lith)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, the printed result says 2010 in the \"strat\" column belongs to \"HoloceneUnits.EC\". The first part relates to the main group the unit belongs to. [`StratGeotop`](../api_reference/geotop_selection.rst) contains four main groups of stratigraphic units as class attributes: \"holocene\", \"channel\", \"older\" and \"antropogenic\". The second part is the abbreviated name of the unit which can be found in the [stratigraphic nomenclature](https://www.dinoloket.nl/stratigrafische-nomenclator/boven-noordzee-groep). Most units belong to the \"Boven-Noordzee Groep\" (prefix \"NU\"). For example, the unit \"EC\" from the selection result corresponds to the \"Echteld formation\" and can be found [here](https://www.dinoloket.nl/stratigrafische-nomenclator/formatie-van-echteld).\n",
    "\n",
    "The printed result of `lith` says 6 in the \"lithok\" column refers to \"Lithology.medium_sand\" which is medium coarse sand.\n",
    "\n",
    "[`StratGeotop`](../api_reference/geotop_selection.rst) and [`Lithology`](../api_reference/geotop_selection.rst) can also easily be transformed into a dictionary so each value in the \"strat\" and \"lithok\" columns can be transformed into meaningful labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "strat_dict = StratGeotop.to_dict(key=\"value\")\n",
    "lith_dict = Lithology.to_dict(key=\"value\")\n",
    "\n",
    "voxel_column[\"strat\"] = voxel_column[\"strat\"].replace(strat_dict)\n",
    "voxel_column[\"lithok\"] = voxel_column[\"lithok\"].replace(lith_dict)\n",
    "voxel_column # See the result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, let's select borehole \"B31H0541\" from the [`BoreholeCollection`](../api_reference/borehole_collection.rst) (this was the borehole at the first index in `geotop_at_locations`) and then compare lithology and stratigraphy information in the borehole data with the voxel column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "borehole = boreholes.get(\"B31H0541\")\n",
    "# Check out relevant information in the borehole data\n",
    "borehole.data[['nr', 'surface', 'end', 'top', 'bottom', 'lith', 'strat_2003']]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Combining with Collections\n",
    "\n",
    "GeoST also makes it possible to add information from GeoTop directly to borehole data in a [`BoreholeCollection`](../api_reference/borehole_collection.rst) or CPT data in a [`CptCollection`](../api_reference/cpt_collection.rst). This is shown in the [examples](../examples/combine_geotop_with_cpts.ipynb) section."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "default",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}