Juyter Logo


Bane Sullivan
%matplotlib inline
from pyvista import set_plot_theme

#Kriging with GSTools

This example utilizes data available from the FORGE geothermal reserach site <https://utahforge.com>'s 2019 Geothermal Design Challenge <https://utahforge.com/studentcomp/>. In this example, the data is archived in the Open Mining Format v1 (OMF) <https://github.com/gmggroup/omf>_ specification and the omfvista <https://opengeovis.github.io/omfvista/>_ software is leverage to load those data into a PyVista MultiBlock data structure.

The goal of this workflow is to create a 3D temperature model by kriging the Observed Temperature data (sparse observational data). The open-source, Python software GSTools <https://geostat-framework.github.io/>_ is used to perform variogram analysis and kriging of the temperature data onto a PyVista mesh to create the 3D model.

# sphinx_gallery_thumbnail_number = 4
import pyvista as pv
from pyvista import examples
import numpy as np
import omfvista
import PVGeo
import matplotlib.pyplot as plt
from gstools import Exponential, krige, vario_estimate_unstructured
from gstools.covmodel.plot import plot_variogram

#Load the Data

For this project, we have two data archives in the Open Mining Format (OMF) <https://github.com/gmggroup/omf>_ specification and we will use one of PyVista's companion projects, omfvista <https://opengeovis.github.io/omfvista/>_ to load those data archives into PyVista a MultiBlock dataset.

url = "https://dl.dropbox.com/s/3cuxvurj8zubchb/FORGE.omf?dl=0"
path, _ = examples.downloads._retrieve_file(url, "FORGE.omf")

project = omfvista.load_project(path)

#Initial Inspection

Now we can go ahead and create an integrated visualization of all of the data available to us.

p = pv.Plotter(window_size=np.array([1024, 768]) * 2)
    project["Site Boundary"], color="yellow", render_lines_as_tubes=True, line_width=10
p.add_mesh(project["Terrain"], texture="geo_aer", opacity=0.7, lighting=False)
p.add_mesh(project["Opal Mound Fault"], color="brown", opacity=0.7)
p.add_mesh(project["Negro Mag Fault"], color="lightblue", opacity=0.7)
    project["Observed Temperature"],
    clim=[10, 270],

p.camera_position = [
    (315661.9406719345, 4234675.528454831, 15167.291249498076),
    (337498.00521202036, 4260818.504034578, -1261.5688408692681),
    (0.2708862567924439, 0.3397398234107863, 0.9006650255615491),


Now we will use an external library, gstools, to krige the temperature probe data into a 3D temperature model of the subsurface.

#Choosing a Model Space

We start out by creating the 3D model space as a PyVista UniformGrid that encompasses the project region. The following were chosen manually:

# Create the kriging grid
grid = pv.UniformGrid()
# Bottom south-west corner
grid.origin = (329700, 4252600, -2700)
# Cell sizes
grid.spacing = (250, 250, 50)
# Number of cells in eaxh direction
grid.dimensions = (60, 75, 100)

Visually inspect the kriging grid in relation to data

p = pv.Plotter(window_size=np.array([1024, 768]) * 2)
p.add_mesh(grid, opacity=0.5, color=True)
p.add_mesh(project["Terrain"], texture="geo_aer", opacity=0.75)
    project["Observed Temperature"],
    stitle="Temperature (C)",
p.camera_position = [
    (303509.4197523619, 4279629.689766085, 8053.049483835099),
    (336316.405356571, 4261479.748583805, -1756.358124546427),
    (0.22299463811939535, -0.11978828465250713, 0.9674317331109259),

#Variogram Analysis

Next, we need to compute a variogram for the temperature probe data and fit that variogram to an exponential model

bins = np.linspace(0, 12300, 1000)
bin_center, gamma = vario_estimate_unstructured(
    project["Observed Temperature"].points.T,
    project["Observed Temperature"]["temperature (C)"],

fit_model = Exponential(dim=3)
fit_model.fit_variogram(bin_center, gamma, nugget=False)

plt.figure(figsize=(10, 5))
plt.plot(bin_center, gamma)
plot_variogram(fit_model, x_max=bins[-1], ax=plt.gca())
plt.xlabel("Lag Distance")
/opt/homebrew/Caskroom/miniforge/base/lib/python3.9/site-packages/gstools/covmodel/plot.py:112: UserWarning: Matplotlib is currently using module://matplotlib_inline.backend_inline, which is a non-GUI backend, so cannot show the figure.
<Figure size 720x360 with 1 Axes>

#Perfroming the Kriging

Then we pass the fitted exponential model when instantiating the kriging operator from GSTools.

# Create the kriging model
krig = krige.Ordinary(
    project["Observed Temperature"].points.T,
    project["Observed Temperature"]["temperature (C)"],

After instantiating the kriging operator, we can have it operate on the nodes of our 3D grid that we created earlier and collect the results back onto the grid.

    name="temperature (C)",
(array([116.47318011, 117.61463724, 118.68608621, ..., 38.14076571, 39.54938107, 40.95250078]), array([8097.97043138, 7671.53933859, 7233.85868221, ..., 8755.23073853, 9108.5391645 , 9445.15771335]))

And now the grid model has the temperature scalar field and kriging variance as data arrays.

project["Kriged Temperature Model"] = grid


We can use a filter from PVGeo to extract the region of the temperature model that is beneath the topography surface as the kriging did not account for the surface boundary.

# Instantiate the algorithm
extractor = PVGeo.grids.ExtractTopography(tolerance=10, remove=True)
# Apply the algorithm to the PyVista grid using the topography surface
subsurface = extractor.apply(grid, project["Terrain"])

#Resulting Visualization

# Create the main figure
def add_contents(p):
    """A helper to add data to scenes."""
        project["Site Boundary"].tube(50),
    p.add_mesh(project["Terrain"], texture="geo_aer", opacity=0.7, lighting=False)

    p.add_mesh(project["Opal Mound Fault"], color="brown", opacity=0.75)
    p.add_mesh(project["Negro Mag Fault"], color="lightblue", opacity=0.75)

        subsurface.ctp().contour([175, 225]),
        name="the model",
        scalars="temperature (C)",
        clim=[10, 270],
        stitle="Temperature (C)",
        project["Observed Temperature"],
        clim=[10, 270],
        stitle="Temperature (C)",


p = pv.Plotter()


p.camera_position = [
    (319663.46494887985, 4294870.4704494225, -8787.973684799075),
    (336926.4650625, 4261892.29012, 103.0),
    (-0.09983848586767283, 0.20995262898057362, 0.9726007250273854),

/opt/homebrew/Caskroom/miniforge/base/lib/python3.9/site-packages/pyvista/plotting/plotting.py:2052: PyvistaDeprecationWarning: 
"stitle" is a depreciated keyword and will be removed in a future

Use ``scalar_bar_args`` instead.  For example:

scalar_bar_args={'title': 'Scalar Bar Title'}

  warnings.warn(USE_SCALAR_BAR_ARGS, PyvistaDeprecationWarning)
ValueError                                Traceback (most recent call last)
Input In [12], in <cell line: 37>()
     32     return
     35 p = pv.Plotter()
---> 37 add_contents(p)
     39 p.camera_position = [
     40     (319663.46494887985, 4294870.4704494225, -8787.973684799075),
     41     (336926.4650625, 4261892.29012, 103.0),
     42     (-0.09983848586767283, 0.20995262898057362, 0.9726007250273854),
     43 ]
     45 p.show()

Input In [12], in add_contents(p)
     10 p.add_mesh(project["Opal Mound Fault"], color="brown", opacity=0.75)
     11 p.add_mesh(project["Negro Mag Fault"], color="lightblue", opacity=0.75)
---> 13 p.add_mesh(
     14     subsurface.ctp().contour([175, 225]),
     15     name="the model",
     16     scalars="temperature (C)",
     17     cmap="coolwarm",
     18     clim=[10, 270],
     19     opacity=0.9,
     20     stitle="Temperature (C)",
     21 )
     22 p.add_mesh(
     23     project["Observed Temperature"],
     24     cmap="coolwarm",
     28     stitle="Temperature (C)",
     29 )
     31 p.enable_depth_peeling(20)

File /opt/homebrew/Caskroom/miniforge/base/lib/python3.9/site-packages/pyvista/plotting/plotting.py:2183, in BasePlotter.add_mesh(self, mesh, color, style, scalars, clim, show_edges, edge_color, point_size, line_width, opacity, flip_scalars, lighting, n_colors, interpolate_before_map, cmap, label, reset_camera, scalar_bar_args, show_scalar_bar, multi_colors, name, texture, render_points_as_spheres, render_lines_as_tubes, smooth_shading, split_sharp_edges, ambient, diffuse, specular, specular_power, nan_color, nan_opacity, culling, rgb, categories, silhouette, use_transparency, below_color, above_color, annotations, pickable, preference, log_scale, pbr, metallic, roughness, render, component, **kwargs)
   2178     mesh, scalars = prepare_smooth_shading(
   2179         mesh, scalars, texture, split_sharp_edges, feature_angle, preference
   2180     )
   2182 if mesh.n_points < 1:
-> 2183     raise ValueError('Empty meshes cannot be plotted. Input mesh has zero points.')
   2185 # set main values
   2186 self.mesh = mesh

ValueError: Empty meshes cannot be plotted. Input mesh has zero points.