3.1: Simple example of kriging in gempy

In this notebook it will be shown how to create a kriged or simulated field in a simple geological model in gempy. We start by creating a simple model with three horizontally layered units, as shown in the gempy examples.

Importing GemPy

import gempy as gp
import gempy_viewer as gpv

# Importing auxiliary libraries
import numpy as np
import matplotlib.pyplot as plt
import os

# new for this
from gempy_plugins.kriging import kriging

np.random.seed(5555)

Creating the model by importing the input data and displaying it:

data_path = os.path.abspath('../../')

geo_data: gp.data.GeoModel = gp.create_geomodel(
    project_name='kriging',
    extent=[0, 1000, 0, 50, 0, 1000],
    resolution=[50, 10, 50],
    refinement=4,
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=data_path + "/data/input_data/jan_models/model1_orientations.csv",
        path_to_surface_points=data_path + "/data/input_data/jan_models/model1_surface_points.csv",
    )
)

Setting and ordering the units and series:

gp.map_stack_to_surfaces(
    gempy_model=geo_data,
    mapping_object={
        "Strat_Series": ('rock2', 'rock1'),
        "Basement_Series": ('basement')
    }
)

Could not find element 'basement' in any group.

Structural Groups:	StructuralGroup: Name: Strat_Series Structural Relation: StackRelationType.ERODE Elements: StructuralElement: Name: rock2 StructuralElement: Name: rock1
Fault Relations:	Strat_Seri... Strat_Series
	True False

Calculating the model:

no mesh computed as basically 2D model

sol = gp.compute_model(geo_data)

Setting Backend To: AvailableBackends.numpy
/home/leguark/gempy/gempy/core/data/geo_model.py:164: UserWarning: You are using refinement and passing a regular grid. The resolution of the regular grid will be overwritten
  warnings.warn(

So here is the very simple, basically 2D model that we created:

gpv.plot_2d(geo_data, cell_number=0, show_data=False)

<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f28bcdeadd0>

1) Creating domain

Let us assume we have a couple of measurements in a domain of interest within our model. In our case the unit of interest is the central rock layer (rock1). In the kriging module we can define the domain by handing over a number of surfaces by id - in this case the id of rock1 is 2. In addition we define four input data points in cond_data, each defined by x,y,z coordinate and a measurement value.

conditioning data (data measured at locations)

cond_data = np.array([[100, .5, 500, 2], [900, .5, 500, 1],
                      [500, .5, 550, 1], [300, .5, 400, 5]])

creating a domain object from the gempy solution, a defined domain conditioning data

domain = kriging.Domain(
    model_solutions=sol,
    transform=geo_data.transform,
    domain=[2],
    data=cond_data
)

2) Creating a variogram model

variogram_model = kriging.VariogramModel(
    theoretical_model='exponential',
    range_=200,
    sill=np.var(cond_data[:, 3])
)

variogram_model.plot(type_='both', show_parameters=True)
plt.show()

3) Kriging interpolation

In the following we define an object called kriging_model and set all input parameters. Finally we generate the kriged field.

kriging_solution = kriging.create_kriged_field(domain, variogram_model)

The result of our calculation is saved in the following dataframe, containing an estimated value and the kriging variance for each point in the grid:

kriging_solution.results_df.head()

	X	Y	Z	estimated value	estimation variance
0	31.2564	1.5689	406.2564	2.296721	1.975385
1	31.2564	1.5689	468.7564	2.109911	1.509387
2	31.2564	1.5689	531.2564	2.004282	1.518538
3	31.2564	1.5689	593.7564	2.008109	2.018219
4	31.2564	4.6939	406.2564	2.296786	1.976096

It is also possible to plot the results in cross section similar to the way gempy models are plotted.

from gempy_viewer.modules.plot_2d.visualization_2d import Plot2D

plot_2d: Plot2D = gpv.plot_2d(
    model=geo_data,
    cell_number=0,
    show_data=False,
    show=False,
    kwargs_lithology={ 'alpha': 0.5 }
)
kriging.plot_kriging_results(
    geo_data=geo_data,
    kriging_solution=kriging_solution,
    plot_2d=plot_2d,
    title='Kriging interpolation: Estimated values',
    result_column=['estimated value']
)

plot_2d_both = gpv.plot_2d(
    model=geo_data,
    cell_number=[0, 0],
    show_data=False,
    show=False,
    kwargs_lithology={ 'alpha': 0.5 }
)

kriging.plot_kriging_results(
    geo_data=geo_data,
    kriging_solution=kriging_solution,
    plot_2d=plot_2d_both,
    title='Kriging interpolation: Estimated values',
    result_column=['estimated value', 'estimation variance']
)

4) Simulated field

Based on the same objects (domain and varigoram model) also a simulated field (stationary Gaussian Field) can be generated. A Sequential Gaussian Simulation approach is applied in this module:

solution_sim = kriging.create_gaussian_field(domain, variogram_model)

solution_sim.results_df.head()

	X	Y	Z	estimated value	estimation variance
256	31.2564	1.5689	406.2564	3.687710	0.228314
652	31.2564	1.5689	468.7564	3.203486	0.153882
460	31.2564	1.5689	531.2564	1.588148	0.281716
108	31.2564	1.5689	593.7564	3.825871	0.746155
745	31.2564	4.6939	406.2564	3.744558	0.041987

solution_sim.results_df['estimated value']

256    3.687710
652    3.203486
460    1.588148
108    3.825871
745    3.744558
         ...
237    3.287792
25     3.305687
946    2.686856
62     2.840951
243    3.332732
Name: estimated value, Length: 1024, dtype: float64

plot_2d: Plot2D = gpv.plot_2d(
    model=geo_data,
    cell_number=0,
    show_data=False,
    show=False,
    kwargs_lithology={ 'alpha': 0.5 }
)
kriging.plot_kriging_results(
    geo_data=geo_data,
    kriging_solution=solution_sim,
    plot_2d=plot_2d,
    title='Kriging interpolation: Estimated values',
    result_column=['estimated value']
)

plot_2d_both = gpv.plot_2d(
    model=geo_data,
    cell_number=[0, 0],
    show_data=False,
    show=False,
    kwargs_lithology={ 'alpha': 0.5 }
)

kriging.plot_kriging_results(
    geo_data=geo_data,
    kriging_solution=solution_sim,
    plot_2d=plot_2d_both,
    title='Kriging interpolation: Estimated values',
    result_column=['estimated value', 'estimation variance']
)
# sphinx_gallery_thumbnail_number = 3