3D viewing and VTK export#

GeoST makes use of PyVista to allow the user to export data to Visualization Toolkit (VTK) format as well as interactively view data in 3D. In this example we create a simple 3D plot of some boreholes and the GeoTOP model in the area of the Utrecht Science Park (USP).

Let’s import GeoST and PyVista and load the data we are going to work with:

import pyvista as pv

import geost

# Load the boreholes from the Utrecht Science Park (USP) area
boreholes_usp = geost.data.boreholes_usp()

# Load the GeoTOP model from the OpenDAP server based on the bounding box of the boreholes
# Load = True because the area is small, so we prefer loading over streaming data
geotop_usp = geost.read_geotop_from_opendap(
    bbox=boreholes_usp.header.total_bounds, load=True
)

# Load REGIS from the OpenDAP server based on the bounding box of the boreholes
regis_usp = geost.read_regis_from_opendap(
    bbox=boreholes_usp.header.total_bounds, load=True
)

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

print("Boreholes:\n")
boreholes_usp

Boreholes:

Collection
  header (rows, columns) : (67, 6)
  data (rows, columns)   : (1398, 32)
crs: RD
vertical datum: NAP height

We have a Collection containing 67 boreholes, an xarray.Dataset of the GeoTOP voxel model and an xarray.Dataset representing the REGIS layer model. GeoST offers a few options to turn these into PyVista objects:

Object	Method	PyVista object	Description
Collection	to_pyvista_cylinders	PyVista.MultiBlock	Fancy cylinder geometries of boreholes
Collection	to_pyvista_grid	PyVista.UnstructuredGrid	Simple geometries, useful for large datasets
xarray `Dataset` or `DataArray` (Voxel model)	gst.to_pyvista_grid(structured=True)	PyVista.ImageData	Structured grid, useful for large datasets
xarray `Dataset` or `DataArray` (Voxel model)	gst.to_pyvista_grid(structured=False)	PyVista.UnstructuredGrid	Individual rendering of each voxel
xarray `Dataset` or `DataArray` (Layer model)	gst.to_pyvista_grid	PyVista.UnstructuredGrid	Individual rendering of each voxel

We will apply to_pyvista_grid for our boreholes to get a PyVista.UnstructuredGrid and the accessor method to_pyvista_grid with structured=True (default) for the GeoTOP model to get a PyVista.ImageData object. The loaded REGIS layer model will be converted to a PyVista.UnstructuredGrid.

Note

If your xarray Dataset or DataArray represents a layer model, the ‘structured’ argument will be ignored as a layer model can only be converted to an unstructured grid

# The boreholes are converted to an unstructured grid of the lithology data. The borehole
# geometries are given a radius of 10 meters.
boreholes_usp_pyvista = boreholes_usp.to_pyvista_grid("lith", radius=10)

# The GeoTOP model is converted to an structured grid of the stratigraphy data
geotop_usp_pyvista = geotop_usp.gst.to_pyvista_grid("strat")

# The REGIS model is converted to an unstructured grid. We first slice it between -80
# and -120 m NAP for better visual separation with GeoTOP later on.
regis_usp = regis_usp.gst.slice_depth_interval(upper=-80, lower=-120)
regis_usp_pyvista = regis_usp.gst.to_pyvista_grid(["kh", "kv"])

# Let's see what the objects look like
print("Boreholes PyVista grid:\n", boreholes_usp_pyvista)
print("\n")
print("GeoTOP PyVista grid:\n", geotop_usp_pyvista)
print("\n")
print("REGIS PyVista grid:\n", regis_usp_pyvista)

Boreholes PyVista grid:
 UnstructuredGrid (0x7a2a58e6cee0)
  N Cells:    1398
  N Points:   11184
  X Bounds:   1.395e+05, 1.413e+05
  Y Bounds:   4.548e+05, 4.560e+05
  Z Bounds:   -2.712e+02, 8.037e+00
  N Arrays:   1


GeoTOP PyVista grid:
 ImageData (0x7a2a58e6f6a0)
  N Cells:      73242
  N Points:     83524
  X Bounds:     1.396e+05, 1.414e+05
  Y Bounds:     4.548e+05, 4.560e+05
  Z Bounds:     -4.975e+01, 1.068e+02
  Dimensions:   19, 14, 314
  Spacing:      1.000e+02, 1.000e+02, 5.000e-01
  N Arrays:     1


REGIS PyVista grid:
 UnstructuredGrid (0x7a2a58e42ce0)
  N Cells:    491
  N Points:   1604
  X Bounds:   1.395e+05, 1.413e+05
  Y Bounds:   4.547e+05, 4.560e+05
  Z Bounds:   -1.200e+02, -8.000e+01
  N Arrays:   2

You can use the save method of these PyVista objects to export them as VTK files and load them into your viewer of choice.

Note

Tip: if an UnstructuredGrid is large (several millions of points and faces) run the clean method with a tolerance of e.g. 1e-5, which merges points that are close to each other. This decreases file size and RAM usage during rendering. Also save the model as vtu instead of vtk file, which is a compressed vtk-format to further save disk-space.

We typically export the PyVista objects to VTK and use ParaView to view data in 3D.

# Example of saving the PyVista objects to files
# boreholes_usp_pyvista.save("boreholes.vtu")
# geotop_usp_pyvista.save("geotop.vtk")
# regis_usp_pyvista.save("regis.vtu")

Paraview visualization example:

GeoTOP and IC model