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 REST-API service. 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.bro_api_read(
    "BHR-P", bbox=(141_000, 455_200, 142_500, 456_000)
)  # xmin, ymin, xmax, ymax

print(bhrp_cores)
print("\nThe header looks like this:")
bhrp_cores.header.gdf.head()

BoreholeCollection:
# header = 85

The header looks like this:

/home/runner/work/geost/geost/geost/validation/validate.py:46: ValidationWarning: 
Validation dropped 2 row(s) for schema 'Point header'.
Dropped indices: [43, 58]

  warnings.warn(
/home/runner/work/geost/geost/geost/validation/validate.py:46: ValidationWarning: 
Validation dropped 8 row(s) for schema 'Layer data non-inclined'.
Dropped indices: [210, 211, 212, 213, 290, 291, 292, 293]

  warnings.warn(

	nr	crs	surface	vertical_datum	end	ghg	glg	landuse	x	y	geometry
0	BHR000000003923	urn:ogc:def:crs:EPSG::28992	1.98	NAP	0.48	0.00	NaN	loofbos	142069.0	455572.0	POINT (142069 455572)
1	BHR000000004193	urn:ogc:def:crs:EPSG::28992	1.85	NAP	0.35	0.10	NaN	graslandBlijvend	142172.0	455444.0	POINT (142172 455444)
2	BHR000000008779	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.35	NaN	loofbos	142052.0	455658.0	POINT (142052 455658)
3	BHR000000009051	urn:ogc:def:crs:EPSG::28992	2.26	NAP	0.76	0.35	NaN	fruitteeltGroen	142244.0	455489.0	POINT (142244 455489)
4	BHR000000021785	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.50	NaN	loofbos	141967.0	455593.0	POINT (141967 455593)

Note that we received two validation warnings because two of the boreholes (BHR000000228773 and BHR000000190070) 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.

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	crs	surface	vertical_datum	end	ghg	glg	landuse	x	y	geometry
0	BHR000000003923	urn:ogc:def:crs:EPSG::28992	1.98	NAP	0.48	0.00	NaN	loofbos	650624.300474	5.773135e+06	POINT (650624.3 5773135.146)
1	BHR000000004193	urn:ogc:def:crs:EPSG::28992	1.85	NAP	0.35	0.10	NaN	graslandBlijvend	650731.448613	5.773011e+06	POINT (650731.449 5773010.605)
2	BHR000000008779	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.35	NaN	loofbos	650604.483467	5.773221e+06	POINT (650604.483 5773220.538)
3	BHR000000009051	urn:ogc:def:crs:EPSG::28992	2.26	NAP	0.76	0.35	NaN	fruitteeltGroen	650801.928242	5.773058e+06	POINT (650801.928 5773057.945)
4	BHR000000021785	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.50	NaN	loofbos	650521.668771	5.773153e+06	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["standardSoilName"].unique()

array(['sterkSiltigeKlei', 'matigSiltigeKlei', 'matigZandigeKlei',
       'kleiigZand', 'uiterstSiltigeKlei', 'sterkZandigeKlei',
       'zwakSiltigZand', 'sterkZandigeLeem', 'zwakZandigeLeem',
       'mineraalarmVeen', 'matigSiltigZand', 'zwakKleiigVeen',
       'sterkSiltigZand'], dtype=object)

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 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”.

bhrp_cores.data["lithology"] = bhrp_cores.data["standardSoilName"].str.extract(
    r"(Klei|Zand|Leem|Veen)$",
    expand=False,  # Find one of these words at the end of the standardSoilName
)

# Check and show that it worked:
bhrp_cores.data["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_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_thickness(
    "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	crs	surface	vertical_datum	end	ghg	glg	landuse	x	y	geometry	Klei_thickness	Zand_thickness
0	BHR000000003923	urn:ogc:def:crs:EPSG::28992	1.98	NAP	0.48	0.00	NaN	loofbos	650624.300474	5.773135e+06	POINT (650624.3 5773135.146)	1.50	0.00
1	BHR000000004193	urn:ogc:def:crs:EPSG::28992	1.85	NAP	0.35	0.10	NaN	graslandBlijvend	650731.448613	5.773011e+06	POINT (650731.449 5773010.605)	1.20	0.30
2	BHR000000008779	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.35	NaN	loofbos	650604.483467	5.773221e+06	POINT (650604.483 5773220.538)	1.50	0.00
3	BHR000000009051	urn:ogc:def:crs:EPSG::28992	2.26	NAP	0.76	0.35	NaN	fruitteeltGroen	650801.928242	5.773058e+06	POINT (650801.928 5773057.945)	1.50	0.00
4	BHR000000021785	urn:ogc:def:crs:EPSG::28992	2.33	NAP	0.83	0.50	NaN	loofbos	650521.668771	5.773153e+06	POINT (650521.669 5773152.781)	0.85	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))

Make this Notebook Trusted to load map: File -> Trust Notebook

Retrieving and analysing BRO soil cores#

This Page