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 Collection.
# import seaborn as sns
import pandas as pd
from matplotlib import pyplot as plt
import geost
from geost.analysis.combine import add_voxelmodel_variable
from geost.bro import GeoTop
from geost.bro.bro_geotop import StratGeotop
cpts = geost.data.cpts_usp()
cpts.header.explore(style_kwds=dict(color="red", weight=6))
Downloading file 'cpts_usp.parquet' from 'https://github.com/Deltares-research/geost/raw/main/data/cpts_usp.parquet' to '/home/runner/.cache/geost'.
/home/runner/work/geost/geost/geost/utils/conversion.py:181: MixedDepthWarning: Data contains a mix of surveys with depth respect to surface and with depth starting at 0. GeoST methods including depth expect all surveys to have a. depth starting at 0 and increasing downward. Please adjust your data accordingly.
warnings.warn(
/home/runner/work/geost/geost/geost/validation/validate.py:77: ValidationWarning:
============================================================
⚠️ VALIDATION ISSUE (1/1)
============================================================
Column : 'depth'
Message : Column 'depth' must be positive downwards (i.e. >= 0), but some rows violate this condition.
# surveys : 1
# rows : 1
============================================================
warnings.warn(
/home/runner/work/geost/geost/geost/base.py:473: AlignmentWarning: Header covers more/other objects than present in the data table. consider running the method 'reset_header' to update the header.
warnings.warn(
✅ Invalid surveys were dropped from the DataFrame because geost.config.validation.DROP_INVALID=True
📖 See the user guide section on validation for advanced handling of validation issues: https://deltares-research.github.io/geost/user_guide/validation.html
NOTE: Header has been reset to align with data because AUTO_ALIGN is enabled in the GeoST configuration.
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.total_bounds)
print(geotop)
GeoTop
Data variables:
strat (y, x, z) float32 260kB dask.array<chunksize=(13, 16, 76), meta=np.ndarray>
Dimensions: {'y': 13, '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 Collection. First, we must create a “depth” column in the CPT data because this is expected to be present in the data. We can use the column “penetration_length” for this. Below we add the variable and print the resulting column:
cpts.data["depth"] = cpts.data["depth"]
cpts = add_voxelmodel_variable(cpts, geotop, "strat")
print(cpts.data["strat"])
0 2010.0
1 2010.0
2 2010.0
3 2010.0
4 2010.0
...
71846 3100.0
71847 3100.0
71848 3100.0
71849 3100.0
71850 3100.0
Name: strat, Length: 71851, dtype: float32
Note that some of the resulting values are “NaN” 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()
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]')
Each number on the x-axis represents a geological unit. However, these numbers are off course not really intuitive and knowing the corresponding geological unit is hard. GeoST provides the StratGeotop class which makes it easy to select desired stratigraphic units with or relabel the unit numbers in the plot above into more meaningfull names.
First, we select all the units that have been merged with the CPT data:
units = StratGeotop.select_values(cpts.data["strat"].unique())
units
[<HoloceneUnits.EC: 2010>,
<ChannelBeltUnits.BEC: 6100>,
<ChannelBeltUnits.CEC: 6200>,
<OlderUnits.BXWISIKO: 3030>,
<OlderUnits.BX: 3100>,
<OlderUnits.KRBXDE: 4010>,
<OlderUnits.DR: 5000>,
<OlderUnits.UR: 5060>,
<OlderUnits.ST: 5070>,
<AntropogenicUnits.AAOP: 1000>]
This returns a list of enum types for each unit. An enum contains a “name” with corresponding “value”. Now, we can make a dictionary with both and replace the numbers with the names in the “strat” column. Let’s print the result and plot the figure again:
replace_dict = {unit.value: unit.name for unit in units}
cpts.data.replace({"strat": replace_dict}, 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 EC
1 EC
2 EC
3 EC
4 EC
..
71846 BX
71847 BX
71848 BX
71849 BX
71850 BX
Name: strat, Length: 71851, dtype: object
Text(0, 0.5, 'Cone resistance [MPa]')
The x-axis in the plot now has some more meaningful abbreviations for geological units which can be found in the stratigraphic nomenclature. Most units belong to the “Boven-Noordzee Groep” (prefix “NU”). For example, the unit “EC” in the plot with the corresponding code can be found here.
The StratGeotop class can be used for all sorts of selections and groupings a user would like for further analyses.