Add stratigraphy from GeoTOP to CPT data#

This example shows how stratigraphic information from GeoTOP can easily be added to CPT data. This allows CPT parameters to be aggregated to calculate averages for geological units. For this example, we will use a selection of CPTs in the area of the Utrecht Science Park (USP).

We will first import the relevant modules and plot the locations of the CPTs in the Collection.

import xyzservices.providers as xyz  # For accessing background map providers
from matplotlib import pyplot as plt

import geost

cpts = geost.data.cpts_usp()
cpts.header.explore(tiles=xyz.nlmaps.grijs, style_kwds={"weight": 4})

Downloading file 'cpts_usp.parquet' from 'https://github.com/Deltares-research/geost/raw/main/data/cpts_usp.parquet' to '/home/runner/.cache/geost'.
Make this Notebook Trusted to load map: File -> Trust Notebook

Adding information from a voxel model#

Any information from voxel models can be added to the CPT data. In this example, we show how to add stratigraphic information from GeoTOP to the CPTs. First, we read GeoTOP directly from the OpenDAP server for the USP area.

geotop = geost.read_geotop_from_opendap(
    data_vars=["strat"], bbox=cpts.header.total_bounds
)
print(geotop)

<xarray.Dataset> Size: 262kB
Dimensions:      (y: 13, x: 16, z: 313)
Coordinates:
  * y            (y) float64 104B 4.548e+05 4.548e+05 ... 4.558e+05 4.56e+05
  * x            (x) float64 128B 1.396e+05 1.396e+05 ... 1.41e+05 1.41e+05
  * z            (z) float32 1kB -49.75 -49.25 -48.75 ... 105.2 105.8 106.2
    spatial_ref  int64 8B 0
Data variables:
    strat        (y, x, z) float32 260kB dask.array<chunksize=(13, 16, 313), meta=np.ndarray>
Attributes:
    title:         geotop v01r6s1 (100.0m * 100.0m * 0.5m
    references:    http://www.dinoloket.nl/detaillering-van-de-bovenste-lagen...
    comment:       GeoTOP 1.6.1 (lithoklasse)
    disclaimer:    http://www.dinoloket.nl
    terms_of_use:  http://www.dinoloket.nl
    institution:   TNO / Geologische Dienst Nederland
    Conventions:   CF-1.4
    source:        Generated by nl.tno.dino.geo3dmodel.generator.GeoTOPVoxelM...
    history:       Generated on Thu May 01 22:54:36 CEST 2025

As you can see, this prints an xarray.DataArray containing the "strat" variable of GeoTOP. We can use this object to add the stratigraphic information to the CPT data. This adds a "strat" column to the data table of the Collection.

First, we need to create a "depth" column in the CPT data, as this is expected to be present. We can use the "penetration_length" column for this. Below, we add the variable and print the resulting column:

cpts.data["depth"] = cpts.data["depth"]
cpts = cpts.add_model_data(geotop, data_vars="strat")
print(cpts.data["strat"])

0        2010.0
1        2010.0
2        2010.0
3        2010.0
4        2010.0
          ...  
73664    3100.0
73665    3100.0
73666    3100.0
73667    3100.0
73668    3100.0
Name: strat, Length: 73331, dtype: float32

Note that some of the resulting values are "NaN", which occurs when a CPT falls outside the 3D model extent. In this case, some CPTs extend deeper than the maximum depth of GeoTOP.

We can now easily aggregate any CPT parameter according to the new "strat" variable. For example, we can create a boxplot of the average cone resistance ("qc"):

fig, ax = plt.subplots()
cpts.data.boxplot(
    ax=ax, column="coneresistance", by="strat", grid=False, showfliers=False
)
ax.tick_params(axis="x", rotation=45)
ax.set_title("")
ax.set_ylabel("Cone resistance [MPa]")

Text(0, 0.5, 'Cone resistance [MPa]')
../_images/1d03b969cf407b3418875afc4999fe84191ced8ece59b05c9e16f9cb5758a008.png

Each number on the x-axis represents a geological unit. However, these numbers are not really intuitive, and it can be difficult to determine which geological unit they represent. GeoST provides an option to quickly load the stratigraphic metadata for GeoTOP supplied by the TNO Geological Survey. This allows us to work with GeoTOP in a more intuitive way and makes it easy to select specific stratigraphic units or relabel the unit numbers in the plot above with more meaningful names.

First, we select all the units that have been added to the CPT data to see which units they are:

strat_units = geost.bro.geotop_strat_units()
units_in_data = strat_units.select_voxel_nr(cpts.data["strat"].unique())
units_in_data

Downloading file 'geotop_v01r6s1_metadata_strat.parquet' from 'https://github.com/Deltares-research/geost/raw/main/data/geotop_v01r6s1_metadata_strat.parquet' to '/home/runner/.cache/geost'.
UnitType.STRAT
STR_UNIT_CD DESCRIPTION SEQ_NR RED_DEC GREEN_DEC BLUE_DEC
VOXEL_NR
2010 NUEC2 Formatie van Echteld, gedeelte onder Formatie ... 28.0 170 255 245
3030 NUBXWI-SI-KO Formatie van Boxtel, laagpakketten van Wierden... 39.0 255 255 80
3100 NUBX Formatie van Boxtel 47.0 255 235 0
5000 NUDR Formatie van Drente 62.0 255 127 80
5060 NUUR2 Formatie van Urk, gedeelte onder Formatie van ... 68.0 189 183 107
6200 NUECgc Formatie van Echteld (geulafzettingen generati... 22.0 170 196 255
4010 NUKR-BXDE Formatie van Kreftenheye en Formatie van Boxte... 49.0 176 48 96
1000 NUAAOP Antropogene afzettingen, opgebrachte grond 1.0 200 200 200
6100 NUECgb Formatie van Echteld (geulafzettingen generati... 11.0 102 205 171
5070 NUST Formatie van Sterksel 69.0 205 92 92

This returns a GeotopUnits instance showing the stratigraphic units that occur in the CPT data. We can now use the unit attribute of the GeotopUnits instance to replace the numbers in the added "strat" column and plot the figure again:

cpts.data.replace({"strat": units_in_data.unit}, inplace=True)
print(cpts.data["strat"])

fig, ax = plt.subplots()
cpts.data.boxplot(
    ax=ax, column="coneresistance", by="strat", grid=False, showfliers=False
)
ax.tick_params(axis="x", rotation=45)
ax.set_title("")
ax.set_ylabel("Cone resistance [MPa]")

0        NUEC2
1        NUEC2
2        NUEC2
3        NUEC2
4        NUEC2
         ...  
73664     NUBX
73665     NUBX
73666     NUBX
73667     NUBX
73668     NUBX
Name: strat, Length: 73331, dtype: object

Text(0, 0.5, 'Cone resistance [MPa]')
../_images/11dc0d655335d459c983c4a57449f8e60b9cb42aa211c981d58474d78642acde.png

The x-axis in the plot now shows more meaningful abbreviations for the geological units, which can be found in the stratigraphic nomenclature. Most units belong to the “Boven-Noordzee Groep” and have the prefix “NU”. For example, the unit “EC” shown in the plot, together with its corresponding code, can be found here.

The GeotopUnits class can be used for various selections and groupings, allowing users to perform further analyses according to their specific needs.