Model data#

Geological model data typically consists of structured, multidimensional data representing two types of models: voxel models and layer models. These models are often stored in NetCDF format, which can be handled efficiently using the popular Python library Xarray.

Reading model data#

GeoST provides several readers for general and specific model data, such as GeoTOP and REGIS II. These return Xarray objects. Modeldata can be read from local NetCDF files or from OPeNDAP services. The table below lists the different readers for the different types of models.

Model	Read function	Returns	Description
Generic modeldata	`read_model_netcdf` `read_model_from_opendap`	`xarray.Dataset` or `xarray.DataArray`	Generic reader for any kind of modeldata available in NetCDF format or from an OPeNDAP service
GeoTOP (BRO)	`read_geotop_netcdf` `read_geotop_from_opendap`	`xarray.Dataset` or `xarray.DataArray`	Reader for GeoTOP from a local NetCDF or the Dinodata OPeNDAP service
REGIS II (BRO)	`read_regis_netcdf` `read_regis_from_opendap`	`xarray.Dataset` or `xarray.DataArray`	Reader for REGIS II from a local NetCDF or the Dinodata OPeNDAP service

The .gst accessor for model data#

GeoST extends Xarray with a .gst accessor that provides GeoST-specific functionality for Xarray’s core Dataset and DataArray objects. The accessor becomes available when GeoST is imported as follows:

import geost

Once imported, the .gst accessor can be used on Dataset and DataArray objects. In general, voxel models and layer models can be distinguished based on their dimensions:

voxelmodel: "x", "y", and "z" dimensions
layermodel: "x", "y", and "layer" dimensions, with "top" and "bottom" information for each "layer"

For the accessor to work, all three dimensions required for either a voxel model or a layer model must be present. Note that the distinction between a voxel model and a layer model is determined by the third dimension: "z" or "layer". This dimension determines whether the data is handled as a voxel model or a layer model.

The dimension names listed above are not required to match exactly, similar to the positional columns of survey data. For example, latitude and longitude coordinates are also valid "x" and "y" dimensions. The table below shows the possible names for each dimension that are recognized by the accessor.

Dimension	Model type	Names
x	voxel/layer	“x”, “xco”, “xcoord”, “longitude”, “lon”, “easting”
y	voxel/layer	“y”, “yco”, “ycoord”, “latitude”, “lat”, “northing”
z	voxel	“z”, “depth”, “elevation”
z	layer	“layer”, “unit”, “horizon”, “stratunit”
top	layer	“top”, “tv_top_nap”, “top_diepte”, “top_depth”, “upperboundary”
bottom	layer	“bottom”, “tv_bottom_nap”, “basis_diepte”, “bottom_depth”, “lowerboundary”

Note

The third dimension in both a voxel model and a layer model is referred to as the z-dimension. In a voxel model, the z-dimension generally represents depth or elevation with respect to NAP. In a layer model, it generally represents a unit (e.g. a stratigraphic unit), while the associated depths are stored in "top" and "bottom" variables.

Note that GeoTOP is an xarray.Dataset instance. In the coordinates section, we can see the names of the dimensions: “x”, “y”, and “z”. We can verify that these dimensions are indeed recognized by GeoST:

print(f"x-dimension: {geotop.gst.x_dim}")
print(f"y-dimension: {geotop.gst.y_dim}")
print(f"z-dimension: {geotop.gst.z_dim}")

x-dimension: x
y-dimension: y
z-dimension: z

Attempting to call the accessor when not all required dimensions can be found will raise an error. An error will also be raised when the z-dimension of a Dataset or DataArray is ambiguous, meaning that it could represent either a voxel model or a layer model: for example, the Dataset contains both a “z” and “layer” dimension.

import numpy as np
import xarray as xr

# Create an example Dataset with "z" and "layer" dimensions.
invalid_model = xr.Dataset(
    data_vars={"data": (("y", "x", "z", "layer"), np.zeros((3, 3, 2, 1)))},
    coords={"x": [0, 1, 2], "y": [2, 1, 0], "z": [0, 1], "layer": [0]},
)
try:
    invalid_model.gst
except geost.exceptions.InvalidModelError as e:
    print(e)

Ambiguous vertical dimension: voxel=z, layer=layer

Working with model data#

We will continue this example with the previously loaded GeoTOP data for the Utrecht Science Park. As said, all GeoST functionality is available through the .gst accessor. For example, we can check the horizontal and vertical bounds of the GeoTOP data:

print(geotop.gst.bounds())
print(geotop.gst.vertical_bounds())

(139500.0, 454700.0, 141500.0, 456300.0)
(-50.0, 10.5)

This prints the bounding box (xmin, ymin, xmax, ymax) and the depth extent (zmin, zmax) of the data. Because we only extend Xarray objects, familiar selection methods .sel and .isel, to select specific coordinates or indices, can be used directly.

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)

<xarray.Dataset> Size: 202kB
Dimensions:      (y: 16, x: 2, z: 121)
Coordinates:
  * y            (y) float64 128B 4.548e+05 4.548e+05 ... 4.562e+05 4.562e+05
  * x            (x) float64 16B 1.4e+05 1.404e+05
  * z            (z) float32 484B -49.75 -49.25 -48.75 ... 9.25 9.75 10.25
    spatial_ref  int64 8B 0
Data variables: (12/13)
    strat        (y, x, z) float32 15kB ...
    lithok       (y, x, z) float32 15kB ...
    kans_1       (y, x, z) float32 15kB ...
    kans_2       (y, x, z) float32 15kB ...
    kans_3       (y, x, z) float32 15kB ...
    kans_4       (y, x, z) float32 15kB ...
    ...           ...
    kans_6       (y, x, z) float32 15kB ...
    kans_7       (y, x, z) float32 15kB ...
    kans_8       (y, x, z) float32 15kB ...
    kans_9       (y, x, z) float32 15kB ...
    onz_lk       (y, x, z) float32 15kB ...
    onz_ls       (y, x, z) float32 15kB ...
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

Using GeoST, we can select within a bounding box very easily.

xmin, ymin, xmax, ymax = 140_000, 455_500, 141_000, 456_000
selection = geotop.gst.select_within_bbox(xmin, ymin, xmax, ymax)
print(selection.gst.bounds())

(140000.0, 455500.0, 141000.0, 456000.0)

If your coordinates are in a different CRS, it also works by specifying the CRS of your bounding box coordinates.

from geost.utils import projections

# Transform coordinates from EPSG:28992 (Amersfoort / RD New) to EPSG:32631 (WGS 84 / UTM zone 31N)
transformer = projections.horizontal_reference_transformer(28992, 32631)
utm_xmin, utm_ymin = transformer.transform(xmin, ymin)
utm_xmax, utm_ymax = transformer.transform(xmax, ymax)

selection = geotop.gst.select_within_bbox(
    utm_xmin, utm_ymin, utm_xmax, utm_ymax, crs=32631
)
print(selection.gst.bounds())

(140000.0, 455400.0, 141100.0, 456100.0)

Overview of attributes and methods#

Many attributes and methods are implemented in the .gst accessor for analyses and spatial selections besides the ones demonstrated above. A short summary of the available attributes and methods is presented here. A complete overview can be found in the ModelDataArray and ModelDataset sections of the API reference.

All methods listed below work with both voxel model and layer model Dataset/DataArray objects.

Attributes

x: x-coordinates of the model
y: y-coordinates of the model
z: z-coordinates of the model
crs: Coordinate reference system (CRS) of the model
resolution: Resolution of the model: (x,y,z)-resolution for a voxelmodel and (x,y)-resolution for a layermodel
bounds: xmin, ymin, xmax, ymax bounding box of the model
vertical_bounds: Depth min, depth max bounds of the model
surface_level: 2D grid of the surface level of the voxel- or layermodel model

Analysis/selection

get_thickness: Calculate the thickness where a model matches a condition
most_common: Determine the most common value (mode) in a model
slice_depth_interval: Slice a specific depth interval from a model

Spatial methods

select_within_bbox: Select data within an xmin, ymin, xmax, ymax bounding box from the model
mask_geometries: Mask specific cells that overlap with point/line/polygon geometries in the model
select_points: Select specific locations from the model