{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Retrieving and analysing BRO soil cores \n",
    "In this example, we will download soil cores (BRO BHR-P objects) and calculate the total\n",
    "thickness of sand, clay, silt and peat layers for each or the boreholes.\n",
    "\n",
    "First we import GeoST and retrieve soil cores from the area of Landgoed Oostbroek\n",
    "(located between Utrecht Science Park and Zeist). GeoST directly connects to the BRO\n",
    "REST-API service. This allows the user to download BRO objects from a spatial \n",
    "query and apply GeoST functionality to the loaded objects. In this case we will use \n",
    "a bounding box:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\onselen\\Development\\geost\\geost\\validation\\validate.py:54: ValidationWarning: \n",
      "Validation failed for schema 'Point header'.\n",
      "Details:\n",
      "non-nullable series 'surface' contains null values:\n",
      "0   NaN\n",
      "Name: surface, dtype: float64\n",
      "\n",
      "  warnings.warn(\n",
      "C:\\Users\\onselen\\Development\\geost\\geost\\validation\\validate.py:54: ValidationWarning: \n",
      "Validation failed for schema 'Layer data non-inclined'.\n",
      "Details:\n",
      "non-nullable series 'surface' contains null values:\n",
      "0   NaN\n",
      "1   NaN\n",
      "2   NaN\n",
      "3   NaN\n",
      "Name: surface, dtype: float64\n",
      "\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BoreholeCollection:\n",
      "# header = 87\n",
      "\n",
      "The header looks like this:\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\onselen\\Development\\geost\\geost\\validation\\validate.py:54: ValidationWarning: \n",
      "Validation failed for schema 'Point header'.\n",
      "Details:\n",
      "non-nullable series 'surface' contains null values:\n",
      "43   NaN\n",
      "58   NaN\n",
      "Name: surface, dtype: float64\n",
      "\n",
      "  warnings.warn(\n",
      "C:\\Users\\onselen\\Development\\geost\\geost\\validation\\validate.py:54: ValidationWarning: \n",
      "Validation failed for schema 'Layer data non-inclined'.\n",
      "Details:\n",
      "non-nullable series 'surface' contains null values:\n",
      "210   NaN\n",
      "211   NaN\n",
      "212   NaN\n",
      "213   NaN\n",
      "290   NaN\n",
      "291   NaN\n",
      "292   NaN\n",
      "293   NaN\n",
      "Name: surface, dtype: float64\n",
      "\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "nr                            BHR000000228773\n",
       "crs               urn:ogc:def:crs:EPSG::28992\n",
       "surface                                   NaN\n",
       "vertical_datum                            NAP\n",
       "begin_depth                               0.0\n",
       "end                                       NaN\n",
       "ghg                                       0.0\n",
       "glg                                       NaN\n",
       "landuse                      graslandBlijvend\n",
       "x                                    141639.0\n",
       "y                                    455437.0\n",
       "geometry                POINT (141639 455437)\n",
       "is_valid                                False\n",
       "Name: 58, dtype: object"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import geost\n",
    "\n",
    "bhrp_cores = geost.bro_api_read(\n",
    "    \"BHR-P\", bbox=(141_000, 455_200, 142_500, 456_000)\n",
    ")  # xmin, ymin, xmax, ymax\n",
    "\n",
    "print(bhrp_cores)\n",
    "print(\"\\nThe header looks like this:\")\n",
    "bhrp_cores.header.gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we received two validation warnings because two of the boreholes (BHR000000228773 and BHR000000190070)\n",
    "do not have elevation data as they are located in a water body. Objects that do not validate are dropped by default in GeoST. See more on validation and validation settings in the [User guide](../user_guide/validation.ipynb).\n",
    "\n",
    "BRO objects come in the Rijksdriehoek coordinate reference system by default (EPSG: 28992),\n",
    "but the example project requires us to work in UTM 31N coordinates (EPSG: 32631). Hence we \n",
    "convert the CRS and check the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bhrp_cores.change_horizontal_reference(32631)  # This is an in-place operation\n",
    "\n",
    "print(bhrp_cores.horizontal_reference)\n",
    "bhrp_cores.header.gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Furthermore, we'd only like to work with boreholes that have a minimum length of 1.5 m:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bhrp_cores = bhrp_cores.select_by_length(\n",
    "    1.5\n",
    ")  # Not an in-place operation, so assign to new or overwrite existing variable!\n",
    "\n",
    "print(bhrp_cores)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like we ditched only two boreholes that were too short.\n",
    "\n",
    "Next, we need to get the main lithology (Zand, Leem, Klei or Veen = Sand, Silt, Clay or Peat) \n",
    "of each described layer in our collection of boreholes. We can use the \"standard_name\" column\n",
    "in the data, which is a lithological description including an indication of admixture:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bhrp_cores.data[\"standardSoilName\"].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we are only interested in the main lithology we need to simplify the names. There \n",
    "is no method for this provided by GeoST, so we will have to create our own. Luckily, GeoST\n",
    "is built on Pandas, which means that you can easily apply your own logic if\n",
    "built-in GeoST methods don't get you where you want. In this case, we will use a string method of a Pandas Series to use a regex (i.e. regular expression) which extracts the main lithology. This will reduce e.g. \"matigSiltigeKlei\" to \"Klei\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bhrp_cores.data[\"lithology\"] = bhrp_cores.data[\"standardSoilName\"].str.extract(\n",
    "    r\"(Klei|Zand|Leem|Veen)$\",\n",
    "    expand=False,  # Find one of these words at the end of the standardSoilName\n",
    ")\n",
    "\n",
    "# Check and show that it worked:\n",
    "bhrp_cores.data[\"lithology\"].unique()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can proceed to compute the total thickness of our main lithologies by applying\n",
    "the built-in method \"get_cumulative_thickness\" on our just created \"main_lithology\" column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# By including the \"include_in_header=True\" argument, the result will be added in-place\n",
    "# to the header table of the BoreholeCollection. If False, the function will only return\n",
    "# a Pandas dataframe with the results.\n",
    "bhrp_cores.get_cumulative_thickness(\n",
    "    \"lithology\", [\"Klei\", \"Zand\", \"Leem\", \"Veen\"], include_in_header=True\n",
    ")\n",
    "\n",
    "# Check if the new columns are present and correct in the header table\n",
    "bhrp_cores.header.gdf.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the purpose of this example we will quickly plot the results on a map below. Hover your mouse\n",
    "over the points to display header information, including the layer thicknesses that we just added.\n",
    "From this point you would normally export the data or continue with further data analyses."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bhrp_cores.header.gdf.explore(style_kwds=dict(color=\"red\", weight=6))"
   ]
  }
 ],
 "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.13.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}