{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Retrieving and analysing BRO soil cores \n", "\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", "
\n", "As of February 2024, BHR-P objects are supported for demonstration purposes. BHR-GT,\n", "BHR-G and BHR-CPT are planned. The BHR-P object parser to GeoST BoreholeCollections is\n", "currently quite slow and will be refactored in a future update.\n", "
\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", "database through an API. 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": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BoreholeCollection:\n", "# header = 85\n", "\n", "The header looks like this:\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nrxysurfaceendgeometry
0BHR0000000039231420694555721.980.48POINT (142069 455572)
1BHR0000000041931421724554441.850.35POINT (142172 455444)
2BHR0000000087791420524556582.330.83POINT (142052 455658)
3BHR0000000090511422444554892.260.76POINT (142244 455489)
4BHR0000000217851419674555932.330.83POINT (141967 455593)
\n", "
" ], "text/plain": [ " nr x y surface end geometry\n", "0 BHR000000003923 142069 455572 1.98 0.48 POINT (142069 455572)\n", "1 BHR000000004193 142172 455444 1.85 0.35 POINT (142172 455444)\n", "2 BHR000000008779 142052 455658 2.33 0.83 POINT (142052 455658)\n", "3 BHR000000009051 142244 455489 2.26 0.76 POINT (142244 455489)\n", "4 BHR000000021785 141967 455593 2.33 0.83 POINT (141967 455593)" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import geost\n", "\n", "bhrp_cores = geost.get_bro_objects_from_bbox(\"BHR-P\", 141000, 142500, 455200, 456000)\n", "\n", "print(bhrp_cores)\n", "print('\\nThe header looks like this:')\n", "bhrp_cores.header.gdf.head()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "EPSG:32631\n" ] }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nrxysurfaceendgeometry
0BHR000000003923650624.3004745.773135e+061.980.48POINT (650624.3 5773135.146)
1BHR000000004193650731.4486135.773011e+061.850.35POINT (650731.449 5773010.605)
2BHR000000008779650604.4834675.773221e+062.330.83POINT (650604.483 5773220.538)
3BHR000000009051650801.9282425.773058e+062.260.76POINT (650801.928 5773057.945)
4BHR000000021785650521.6687715.773153e+062.330.83POINT (650521.669 5773152.781)
\n", "
" ], "text/plain": [ " nr x y surface end \\\n", "0 BHR000000003923 650624.300474 5.773135e+06 1.98 0.48 \n", "1 BHR000000004193 650731.448613 5.773011e+06 1.85 0.35 \n", "2 BHR000000008779 650604.483467 5.773221e+06 2.33 0.83 \n", "3 BHR000000009051 650801.928242 5.773058e+06 2.26 0.76 \n", "4 BHR000000021785 650521.668771 5.773153e+06 2.33 0.83 \n", "\n", " geometry \n", "0 POINT (650624.3 5773135.146) \n", "1 POINT (650731.449 5773010.605) \n", "2 POINT (650604.483 5773220.538) \n", "3 POINT (650801.928 5773057.945) \n", "4 POINT (650521.669 5773152.781) " ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "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()\n" ] }, { "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": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "BoreholeCollection:\n", "# header = 83\n" ] } ], "source": [ "bhrp_cores = bhrp_cores.select_by_length(1.5) # 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": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "\n", "[ 'sterkSiltigeKlei', 'matigSiltigeKlei', 'matigZandigeKlei',\n", " 'kleiigZand', 'uiterstSiltigeKlei', 'sterkZandigeKlei',\n", " 'zwakSiltigZand', 'sterkZandigeLeem', 'zwakZandigeLeem',\n", " 'mineraalarmVeen', 'matigSiltigZand', 'zwakKleiigVeen',\n", " 'sterkSiltigZand']\n", "Length: 13, dtype: string" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bhrp_cores.data[\"standard_name\"].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 create a function\n", "that makes use of regular expression and apply it to the dataframe with layer data (bhrp_cores.data.df).\n", "The function will reduce e.g. \"matigSiltigeKlei\" to \"Klei\"." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array(['Klei', 'Zand', 'Leem', 'Veen'], dtype=object)" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import re\n", "\n", "# Create a regex function that returns \"Klei\", \"Zand\", \"Leem\", or \"Veen\" if it is found\n", "# at the end of a string.\n", "get_main_lithology = lambda x: re.search(r\"(Klei|Zand|Leem|Veen)$\", x).group()\n", "\n", "# Apply the function to the column \"standard_name\" and create a new column\n", "# \"main_lithology\" to store the result.\n", "bhrp_cores.data.df[\"main_lithology\"] = bhrp_cores.data.df[\"standard_name\"].apply(get_main_lithology)\n", "\n", "# Check and show that it worked:\n", "bhrp_cores.data.df[\"main_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_layer_thickness\" on our just created \"main_lithology\" column." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
nrxysurfaceendgeometryKlei_thicknessLeem_thicknessVeen_thicknessZand_thickness
0BHR000000003923650624.3004745.773135e+061.980.48POINT (650624.3 5773135.146)1.500.00.00.00
1BHR000000004193650731.4486135.773011e+061.850.35POINT (650731.449 5773010.605)1.200.00.00.30
2BHR000000008779650604.4834675.773221e+062.330.83POINT (650604.483 5773220.538)1.500.00.00.00
3BHR000000009051650801.9282425.773058e+062.260.76POINT (650801.928 5773057.945)1.500.00.00.00
4BHR000000021785650521.6687715.773153e+062.330.83POINT (650521.669 5773152.781)0.850.00.00.65
\n", "
" ], "text/plain": [ " nr x y surface end \\\n", "0 BHR000000003923 650624.300474 5.773135e+06 1.98 0.48 \n", "1 BHR000000004193 650731.448613 5.773011e+06 1.85 0.35 \n", "2 BHR000000008779 650604.483467 5.773221e+06 2.33 0.83 \n", "3 BHR000000009051 650801.928242 5.773058e+06 2.26 0.76 \n", "4 BHR000000021785 650521.668771 5.773153e+06 2.33 0.83 \n", "\n", " geometry Klei_thickness Leem_thickness \\\n", "0 POINT (650624.3 5773135.146) 1.50 0.0 \n", "1 POINT (650731.449 5773010.605) 1.20 0.0 \n", "2 POINT (650604.483 5773220.538) 1.50 0.0 \n", "3 POINT (650801.928 5773057.945) 1.50 0.0 \n", "4 POINT (650521.669 5773152.781) 0.85 0.0 \n", "\n", " Veen_thickness Zand_thickness \n", "0 0.0 0.00 \n", "1 0.0 0.30 \n", "2 0.0 0.00 \n", "3 0.0 0.00 \n", "4 0.0 0.65 " ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "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_layer_thickness(\n", " \"main_lithology\",\n", " [\"Klei\", \"Zand\", \"Leem\", \"Veen\"],\n", " 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": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
Make this Notebook Trusted to load map: File -> Trust Notebook
" ], "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "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.12.7" } }, "nbformat": 4, "nbformat_minor": 2 }