Add stratigraphy from GeoTOP to CPT data

This example shows how stratigraphic layer boundaries from GeoTOP can easily be added to CPT data. This way, CPT parameters can easily be aggregated to get averages for geological units. For this example we are going to 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 CptCollection.

from pathlib import Path

import seaborn as sns
from matplotlib import pyplot as plt

import geost
from geost.analysis.combine import add_voxelmodel_variable
from geost.bro import GeoTop

cpts = geost.read_cpt_table(r'c:\Users\knaake\OneDrive - Stichting Deltares\Documents\data\cpt_data_usp.parquet')
cpts.header.gdf.explore(style_kwds=dict(color="red", weight=6))

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

Adding information from a voxelmodel

Any information from voxelmodels can be added. For this example we show how to add the stratigraphy from GeoTOP to the CPTs. First we read GeoTOP directly from the OpenDaP server for the USP area.

geotop = GeoTop.from_opendap(data_vars=['strat'], bbox=cpts.header.gdf.total_bounds)
print(geotop)

GeoTop
Data variables:
    strat    (y, x, z) float32 240kB ...
Dimensions: {'y': 12, 'x': 16, 'z': 313}
Resolution (y, x, z): (100.0, 100.0, 0.5)

As you can see, this prints a GeoTop instance along with the dimensions and resolution of GeoTOP. We can use this GeoTop instance to add the variable “strat” to the CPT data. This adds a column “strat” to the data object of the CptCollection. Below we add the variable and print the resulting column:

result = add_voxelmodel_variable(cpts, geotop, "strat")
print(result.data['strat'])

0        1000.0
1        1000.0
2        2010.0
3        2010.0
4        3030.0
          ...  
74725       NaN
74726       NaN
74727       NaN
74728       NaN
74729       NaN
Name: strat, Length: 74632, dtype: float32

Note that some of the resulting values are “NaN” (Not a Number) which occurs when a CPT falls outside of the 3D model extent. In this case, some CPTs are deeper than the maximum depth of GeoTOP. Now we could easily aggregate any CPT parameter according to the new “strat” variable. For example, make a boxplot of the average cone resistance (“qc”):

fig, ax = plt.subplots()
sns.boxplot(ax=ax, data=result.data.df, x="strat", y="cone_resistance", showfliers=False)
ax.tick_params(axis='x', rotation=45)

Each number on the x-axis represents a geological unit. These can be grouped in any way off course a user would like for further analyses. These however will not be shown in this tutorial.