BRO GeoTop

GeoST provides support to work with models that are provided by the BRO such as GeoTOP. This section will explain the usage of GeoTop, how to make easy selections and combinations with for example borehole data contained in a BoreholeCollection.

Reading GeoTop

GeoTOP data can be read from a local NetCDF file of downloaded data from the BRO or directly accessed from the OPeNDAP-server of TNO Geological survey. See GeoTop.from_netcdf and GeoTop.from_opendap in the API reference.

Working with GeoTop

Let’s load some example GeoTOP data for the Utrecht Science park:

import geost

geotop = geost.data.geotop_usp()
print(geotop)

Downloading file 'geotop_usp.nc' from 'https://github.com/Deltares-research/geost/raw/main/data/geotop_usp.nc' to '/home/runner/.cache/geost'.

GeoTop
Data variables:
    strat    (y, x, z) float32 357kB ...
    lithok   (y, x, z) float32 357kB ...
    kans_1   (y, x, z) float32 357kB ...
    kans_2   (y, x, z) float32 357kB ...
    kans_3   (y, x, z) float32 357kB ...
    kans_4   (y, x, z) float32 357kB ...
    kans_5   (y, x, z) float32 357kB ...
    kans_6   (y, x, z) float32 357kB ...
    kans_7   (y, x, z) float32 357kB ...
    kans_8   (y, x, z) float32 357kB ...
    kans_9   (y, x, z) float32 357kB ...
    onz_lk   (y, x, z) float32 357kB ...
    onz_ls   (y, x, z) float32 357kB ...
Dimensions: {'x': 19, 'y': 15, 'z': 313}
Resolution (y, x, z): (100.0, 100.0, 0.5)

As you can see this prints information on the type of class (i.e. GeoTop), the available data variables, xyz-dimensions and resolution. The GeoTop class inherits from GeoST VoxelModel and as such, inherits all the attributes and associated selection methods. For example, we can check the horizontal and vertical bounds:

print(geotop.horizontal_bounds)
print(geotop.vertical_bounds)

(139500.0, 454700.0, 141400.0, 456200.0)
(-50.0, 106.5)

This prints the “xmin, ymin, xmax, ymax” bounding box and the “zmin, zmax” of the data. GeoTop depends strongly on xarray.Dataset and therefore provides familiar selection methods .sel and .isel to select specific coordinates or indices, or slices of coordinates or indices:

sel = geotop.sel(x=[139_650, 139_750]) # Select specific x-coordinates
sel = geotop.isel(x=[2, 5]) # Select specific x-indices

sel = geotop.sel(x=slice(140_000, 140_500)) # Select a slice of x-coordinates
sel = geotop.isel(x=slice(2, 5)) # Select a slice of x-indices

sel = geotop.sel(x=[140_013.23, 140_333.14], method='nearest') # Select the nearest x-coordinates
print(sel)

GeoTop
Data variables:
    strat    (y, x, z) float32 38kB ...
    lithok   (y, x, z) float32 38kB ...
    kans_1   (y, x, z) float32 38kB ...
    kans_2   (y, x, z) float32 38kB ...
    kans_3   (y, x, z) float32 38kB ...
    kans_4   (y, x, z) float32 38kB ...
    kans_5   (y, x, z) float32 38kB ...
    kans_6   (y, x, z) float32 38kB ...
    kans_7   (y, x, z) float32 38kB ...
    kans_8   (y, x, z) float32 38kB ...
    kans_9   (y, x, z) float32 38kB ...
    onz_lk   (y, x, z) float32 38kB ...
    onz_ls   (y, x, z) float32 38kB ...
Dimensions: {'x': 2, 'y': 15, 'z': 313}
Resolution (y, x, z): (100.0, 300.0, 0.5)

As you can see, a selection returns a new GeoTop instance for the desired coordinates. Note that the x-resolution has also been changed because the coordinates in the last selection are approximately 300 meters apart.

Selections

Often it can be necessary to select GeoTop at specific locations, such as borehole locations, to compare the borehole data with GeoTop. We can easily select GeoTop at the locations of boreholes using the GeoTop.select_with_points method. Let’s load a BoreholeCollection with borehole data for the Utrecht Science Park to see what happens:

boreholes = geost.data.boreholes_usp() # BoreholeCollection

geotop_at_locations = geotop.select_with_points(boreholes.header.gdf)
print(geotop_at_locations)

<xarray.Dataset> Size: 1MB
Dimensions:      (idx: 67, z: 313)
Coordinates:
    x            (idx) float64 536B 1.396e+05 1.396e+05 ... 1.41e+05 1.412e+05
    y            (idx) float64 536B 4.56e+05 4.55e+05 ... 4.56e+05 4.56e+05
  * z            (z) float32 1kB -49.75 -49.25 -48.75 ... 105.2 105.8 106.2
    spatial_ref  int64 8B 0
  * idx          (idx) int64 536B 0 1 2 3 4 5 6 7 8 ... 59 60 61 62 63 64 65 66
Data variables: (12/13)
    strat        (idx, z) float32 84kB ...
    lithok       (idx, z) float32 84kB ...
    kans_1       (idx, z) float32 84kB ...
    kans_2       (idx, z) float32 84kB ...
    kans_3       (idx, z) float32 84kB ...
    kans_4       (idx, z) float32 84kB ...
    ...           ...
    kans_6       (idx, z) float32 84kB ...
    kans_7       (idx, z) float32 84kB ...
    kans_8       (idx, z) float32 84kB ...
    kans_9       (idx, z) float32 84kB ...
    onz_lk       (idx, z) float32 84kB ...
    onz_ls       (idx, z) float32 84kB ...
Attributes:
    title:         geotop v01r6 (100.0m * 100.0m * 0.5m
    references:    http://www.dinoloket.nl/detaillering-van-de-bovenste-lagen...
    comment:       GeoTOP 1.6 (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 Oct 05 17:01:25 CEST 2023

The result is an xarray.Dataset with a new “idx” dimension and note that variables like “strat” now have “(idx, z)” dimensions. The .select_with_points method selects the voxel columns each borehole is in, in the index order of the boreholes in the BoreholeCollection. It does not return a VoxelModel instance because a VoxelModel is only valid if it contains “x”, “y” and “z” dimensions.

With the selection result, it is possible to compare a borehole with a specific voxel column. Let’s for example compare the first borehole in the collection with the voxel column. First, let’s print the “header” attribute of the borehole collection to find the “nr” of the first borehole:

print(boreholes.header)

PointHeader instance containing 67 objects
          nr       x       y  surface      end               geometry
0   B31H0541  139585  456000    1.200   -9.900  POINT (139585 456000)
1   B31H0611  139600  455060    1.200  -23.000  POINT (139600 455060)
2   B31H0718  139950  455200    1.300 -271.200  POINT (139950 455200)
3   B31H0803  139675  455087    2.160   -4.840  POINT (139675 455087)
4   B31H0806  139684  455384    1.000  -49.500  POINT (139684 455384)
..       ...     ...     ...      ...      ...                    ...
62  B32C1889  140919  456000    1.728   -2.772  POINT (140919 456000)
63  B32C1890  140938  455878    3.948    0.248  POINT (140938 455878)
64  B32C1891  141025  455998    2.112   -1.888  POINT (141025 455998)
65  B32C1892  141091  455966    1.769   -2.731  POINT (141091 455966)
66  B32C1893  141266  455989    2.328   -2.172  POINT (141266 455989)

[67 rows x 6 columns]

We can see that the first number (at index 0) is “B31H0541”. Notice in the “surface” and “end” columns that the borehole is between +1.2 m NAP and -9.9 m NAP. We can select the first index from geotop_at_locations, select the depths between 1.2 and -9.9 m NAP and create a DataFrame to compare the data. In similar way, selection method for other use cases are available when working with GeoTop. See the API reference for available methods.

# Increase zmin and zmax by a little to make sure the complete depth of the borehole is covered
voxel_column = geotop_at_locations.sel(idx=0, z=slice(-10.5, 1.5))
voxel_column = voxel_column.to_dataframe().sort_index(ascending=False) # Make depth increase from top to bottom
voxel_column # Check out the result

	x	y	strat	lithok	kans_1	kans_2	kans_3	kans_4	kans_5	kans_6	kans_7	kans_8	kans_9	onz_lk	onz_ls	spatial_ref	idx
z
1.25	139550.0	456050.0	2010.0	2.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
0.75	139550.0	456050.0	2010.0	1.0	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
0.25	139550.0	456050.0	2010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	2.0	0	0
-0.25	139550.0	456050.0	2010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	11.0	0	0
-0.75	139550.0	456050.0	3030.0	5.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	17.0	0	0
-1.25	139550.0	456050.0	3030.0	5.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	7.0	0	0
-1.75	139550.0	456050.0	3030.0	5.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
-2.25	139550.0	456050.0	3030.0	5.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	4.0	0	0
-2.75	139550.0	456050.0	3030.0	5.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	13.0	0	0
-3.25	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	17.0	0	0
-3.75	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	13.0	0	0
-4.25	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	3.0	0	0
-4.75	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-5.25	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-5.75	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	3.0	0	0
-6.25	139550.0	456050.0	3100.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	16.0	0	0
-6.75	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	12.0	0	0
-7.25	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-7.75	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-8.25	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-8.75	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-9.25	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-9.75	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	14.0	77.0	9.0	0.0	0.0	36.0	0.0	0	0
-10.25	139550.0	456050.0	4010.0	6.0	0.0	0.0	0.0	0.0	27.0	68.0	5.0	0.0	0.0	39.0	0.0	0	0

Relabelling GeoTOP information

Notice that in the voxel_column result, the columns “strat” (i.e. stratigraphy) and “lithok” (i.e. lithology) contain numbers, instead of names which makes it difficult to understand what the unit and lithology at specific depths are. This makes a comparison with the borehole data difficult.

GeoST provides the StratGeotop class and Lithology class which make it easy to translate these numbers into more meaningful names. Let’s import them and check which stratigraphic unit belongs to the number 2010 in the “strat” column and what the number 6 in the “lithok” means:

from geost.bro import Lithology, StratGeotop

unit = StratGeotop.select_values(2010)
print(unit)

lith = Lithology.select_values(6)
print(lith)

[<HoloceneUnits.EC: 2010>]
[<Lithology.medium_sand: 6>]

As you can see, the printed result says 2010 in the “strat” column belongs to “HoloceneUnits.EC”. The first part relates to the main group the unit belongs to. StratGeotop contains four main groups of stratigraphic units as class attributes: “holocene”, “channel”, “older” and “antropogenic”. The second part is the abbreviated name of the unit which can be found in the stratigraphic nomenclature. Most units belong to the “Boven-Noordzee Groep” (prefix “NU”). For example, the unit “EC” from the selection result corresponds to the “Echteld formation” and can be found here.

The printed result of lith says 6 in the “lithok” column refers to “Lithology.medium_sand” which is medium coarse sand.

StratGeotop and Lithology can also easily be transformed into a dictionary so each value in the “strat” and “lithok” columns can be transformed into meaningful labels:

strat_dict = StratGeotop.to_dict(key="value")
lith_dict = Lithology.to_dict(key="value")

voxel_column["strat"] = voxel_column["strat"].replace(strat_dict)
voxel_column["lithok"] = voxel_column["lithok"].replace(lith_dict)
voxel_column # See the result

	x	y	strat	lithok	kans_1	kans_2	kans_3	kans_4	kans_5	kans_6	kans_7	kans_8	kans_9	onz_lk	onz_ls	spatial_ref	idx
z
1.25	139550.0	456050.0	EC	clay	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
0.75	139550.0	456050.0	EC	organic	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
0.25	139550.0	456050.0	EC	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	2.0	0	0
-0.25	139550.0	456050.0	EC	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	11.0	0	0
-0.75	139550.0	456050.0	BXWISIKO	fine_sand	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	17.0	0	0
-1.25	139550.0	456050.0	BXWISIKO	fine_sand	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	7.0	0	0
-1.75	139550.0	456050.0	BXWISIKO	fine_sand	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0.0	0	0
-2.25	139550.0	456050.0	BXWISIKO	fine_sand	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	4.0	0	0
-2.75	139550.0	456050.0	BXWISIKO	fine_sand	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	13.0	0	0
-3.25	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	17.0	0	0
-3.75	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	13.0	0	0
-4.25	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	3.0	0	0
-4.75	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-5.25	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-5.75	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	3.0	0	0
-6.25	139550.0	456050.0	BX	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	16.0	0	0
-6.75	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	12.0	0	0
-7.25	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-7.75	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-8.25	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-8.75	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-9.25	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	0.0	100.0	0.0	0.0	0.0	0.0	0.0	0	0
-9.75	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	14.0	77.0	9.0	0.0	0.0	36.0	0.0	0	0
-10.25	139550.0	456050.0	KRBXDE	medium_sand	0.0	0.0	0.0	0.0	27.0	68.0	5.0	0.0	0.0	39.0	0.0	0	0

Now, let’s select borehole “B31H0541” from the BoreholeCollection (this was the borehole at the first index in geotop_at_locations) and then compare lithology and stratigraphy information in the borehole data with the voxel column.

borehole = boreholes.get("B31H0541")
# Check out relevant information in the borehole data
borehole.data[['nr', 'surface', 'end', 'top', 'bottom', 'lith', 'strat_2003']]

	nr	surface	end	top	bottom	lith	strat_2003
0	B31H0541	1.2	-9.9	0.00	0.20	K	EC
1	B31H0541	1.2	-9.9	0.20	0.60	K	EC
2	B31H0541	1.2	-9.9	0.60	0.95	V	NI
3	B31H0541	1.2	-9.9	0.95	2.80	Z	EC
4	B31H0541	1.2	-9.9	2.80	4.20	Z	BXWI
5	B31H0541	1.2	-9.9	4.20	8.40	Z	BX
6	B31H0541	1.2	-9.9	8.40	9.60	Z	KR
7	B31H0541	1.2	-9.9	9.60	11.10	Z	KR

Combining with Collections

GeoST also makes it possible to add information from GeoTop directly to borehole data in a BoreholeCollection or CPT data in a CptCollection. This is shown in the examples section.