Retrieving and analysing BRO soil cores#
In this example, we will download soil cores (BRO BHR-P objects) and calculate the total thickness of sand, clay, silt and peat layers for each or the boreholes.
First we import GeoST and retrieve soil cores from the area of Landgoed Oostbroek (located between Utrecht Science Park and Zeist). GeoST directly connects to the BRO database through an API. This allows the user to download BRO objects from a spatial query and apply GeoST functionality to the loaded objects. In this case we will use a bounding box:
import geost
bhrp_cores = geost.get_bro_objects_from_bbox("BHR-P", 141000, 142500, 455200, 456000)
print(bhrp_cores)
print('\nThe header looks like this:')
bhrp_cores.header.gdf.head()
BoreholeCollection:
# header = 85
The header looks like this:
| nr | x | y | surface | end | geometry | |
|---|---|---|---|---|---|---|
| 0 | BHR000000003923 | 142069 | 455572 | 1.98 | 0.48 | POINT (142069 455572) |
| 1 | BHR000000004193 | 142172 | 455444 | 1.85 | 0.35 | POINT (142172 455444) |
| 2 | BHR000000008779 | 142052 | 455658 | 2.33 | 0.83 | POINT (142052 455658) |
| 3 | BHR000000009051 | 142244 | 455489 | 2.26 | 0.76 | POINT (142244 455489) |
| 4 | BHR000000021785 | 141967 | 455593 | 2.33 | 0.83 | POINT (141967 455593) |
BRO objects come in the Rijksdriehoek coordinate reference system by default (EPSG: 28992), but the example project requires us to work in UTM 31N coordinates (EPSG: 32631). Hence we convert the CRS and check the result:
bhrp_cores.change_horizontal_reference(32631) # This is an in-place operation
print(bhrp_cores.horizontal_reference)
bhrp_cores.header.gdf.head()
EPSG:32631
| nr | x | y | surface | end | geometry | |
|---|---|---|---|---|---|---|
| 0 | BHR000000003923 | 650624.300474 | 5.773135e+06 | 1.98 | 0.48 | POINT (650624.3 5773135.146) |
| 1 | BHR000000004193 | 650731.448613 | 5.773011e+06 | 1.85 | 0.35 | POINT (650731.449 5773010.605) |
| 2 | BHR000000008779 | 650604.483467 | 5.773221e+06 | 2.33 | 0.83 | POINT (650604.483 5773220.538) |
| 3 | BHR000000009051 | 650801.928242 | 5.773058e+06 | 2.26 | 0.76 | POINT (650801.928 5773057.945) |
| 4 | BHR000000021785 | 650521.668771 | 5.773153e+06 | 2.33 | 0.83 | POINT (650521.669 5773152.781) |
Furthermore, we’d only like to work with boreholes that have a minimum length of 1.5 m
bhrp_cores = bhrp_cores.select_by_length(1.5) # Not an in-place operation, so assign to new or overwrite existing variable!
print(bhrp_cores)
BoreholeCollection:
# header = 83
It looks like we ditched only two boreholes that were too short.
Next, we need to get the main lithology (Zand, Leem, Klei or Veen = Sand, Silt, Clay or Peat) of each described layer in our collection of boreholes. We can use the “standard_name” column in the data, which is a lithological description including an indication of admixture:
bhrp_cores.data["standard_name"].unique()
<StringArray>
[ 'sterkSiltigeKlei', 'matigSiltigeKlei', 'matigZandigeKlei',
'kleiigZand', 'uiterstSiltigeKlei', 'sterkZandigeKlei',
'zwakSiltigZand', 'sterkZandigeLeem', 'zwakZandigeLeem',
'mineraalarmVeen', 'matigSiltigZand', 'zwakKleiigVeen',
'sterkSiltigZand']
Length: 13, dtype: string
Since we are only interested in the main lithology we need to simplify the names. There is no method for this provided by GeoST, so we will have to create our own. Luckily, GeoST is built on Pandas, which means that you can easily apply your own logic if built-in GeoST methods don’t get you where you want. In this case we will create a function that makes use of regular expression and apply it to the dataframe with layer data (bhrp_cores.data.df). The function will reduce e.g. “matigSiltigeKlei” to “Klei”.
import re
# Create a regex function that returns "Klei", "Zand", "Leem", or "Veen" if it is found
# at the end of a string.
get_main_lithology = lambda x: re.search(r"(Klei|Zand|Leem|Veen)$", x).group()
# Apply the function to the column "standard_name" and create a new column
# "main_lithology" to store the result.
bhrp_cores.data.df["main_lithology"] = bhrp_cores.data.df["standard_name"].apply(get_main_lithology)
# Check and show that it worked:
bhrp_cores.data.df["main_lithology"].unique()
array(['Klei', 'Zand', 'Leem', 'Veen'], dtype=object)
Now we can proceed to compute the total thickness of our main lithologies by applying the built-in method “get_cumulative_layer_thickness” on our just created “main_lithology” column.
# By including the "include_in_header=True" argument, the result will be added in-place
# to the header table of the BoreholeCollection. If False, the function will only return
# a Pandas dataframe with the results.
bhrp_cores.get_cumulative_layer_thickness(
"main_lithology",
["Klei", "Zand", "Leem", "Veen"],
include_in_header=True,
)
# Check if the new columns are present and correct in the header table
bhrp_cores.header.gdf.head()
| nr | x | y | surface | end | geometry | Klei_thickness | Leem_thickness | Veen_thickness | Zand_thickness | |
|---|---|---|---|---|---|---|---|---|---|---|
| 0 | BHR000000003923 | 650624.300474 | 5.773135e+06 | 1.98 | 0.48 | POINT (650624.3 5773135.146) | 1.50 | 0.0 | 0.0 | 0.00 |
| 1 | BHR000000004193 | 650731.448613 | 5.773011e+06 | 1.85 | 0.35 | POINT (650731.449 5773010.605) | 1.20 | 0.0 | 0.0 | 0.30 |
| 2 | BHR000000008779 | 650604.483467 | 5.773221e+06 | 2.33 | 0.83 | POINT (650604.483 5773220.538) | 1.50 | 0.0 | 0.0 | 0.00 |
| 3 | BHR000000009051 | 650801.928242 | 5.773058e+06 | 2.26 | 0.76 | POINT (650801.928 5773057.945) | 1.50 | 0.0 | 0.0 | 0.00 |
| 4 | BHR000000021785 | 650521.668771 | 5.773153e+06 | 2.33 | 0.83 | POINT (650521.669 5773152.781) | 0.85 | 0.0 | 0.0 | 0.65 |
For the purpose of this example we will quickly plot the results on a map below. Hover your mouse over the points to display header information, including the layer thicknesses that we just added. From this point you would normally export the data or continue with further data analyses.
bhrp_cores.header.gdf.explore(style_kwds=dict(color="red", weight=6))