Gempy is an open-source library for modeling geology written in Python. The library is capable of creating complex 3D geological models including structures, fault networks, and unconformities and it can be coupled with uncertainty analysis.

We have created a tutorial of geological modeling based on geological contacts and surface orientations. The tutorial was developed in a Gempy container that runs under Docker in Windows 10; the tutorial covers the software installation and the geological modeling scripting.

Tutorial

Part 1: Tutorial to install Docker in Windows 10

Part2: Tutorial for the installation of Gempy in Windows 10 with Docker

Part3: Basic tutorial of geological modeling with Python and Gempy

Input data

You can download the input data from this link.

Scripting

This is the complete Python code for the tutorial:

#import required packages
import gempy as gp
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import os

WARNING (theano.tensor.blas): Using NumPy C-API based implementation for BLAS functions.

#create gempy model
geo_model = gp.create_model('Model1')

/opt/conda/lib/python3.7/site-packages/theano/gpuarray/dnn.py:184: UserWarning: Your cuDNN version is more recent than Theano. If you encounter problems, try updating Theano or downgrading cuDNN to a version >= v5 and <= v7.
  warnings.warn("Your cuDNN version is more recent than "

#create grid
gp.init_data(geo_model, [0, 2000., 0, 2000., 0, 750.],[50, 50, 50],
            path_o = 'simple_fault_model_orientations.csv',
            path_i = 'simple_fault_model_points.csv',
            default_values = True)

Active grids: ['regular']





Model1  2020-10-02 18:58

#get unsorted surfaces
geo_model.surfaces

	X	Y	Z	smooth	surface
0	1000	50	450.000000	0.000002	Shale
1	1000	150	433.333333	0.000002	Shale
2	1000	300	433.333333	0.000002	Shale
3	1000	500	466.666667	0.000002	Shale
4	1000	1000	533.333333	0.000002	Shale

#get information about orientations
gp.get_data(geo_model, 'orientations').head()

	X	Y	Z	G_x	G_y	G_z	smooth	surface
0	1000	1000	300	0.316229	1.000019e-12	0.948683	0.01	Shale
1	400	1000	420	0.316229	1.000019e-12	0.948683	0.01	Sandstone_2
2	500	1000	300	-0.948683	9.998257e-13	0.316229	0.01	Main_Fault

#define stack from layers
gp.map_stack_to_surfaces(geo_model,
                        {"Fault_Series":"Main_Fault",
                         "Strat_Series":('Sandstone_2','Siltstone','Shale','Sandstone_1','basement')},
                        remove_unused_series=True)

	Fault_Series	Strat_Series
Fault_Series	False	True
Strat_Series	False	False

#review fault information
geo_model.faults

	order_series	BottomRelation	isActive	isFault	isFinite
Fault_Series	1	Fault	True	True	False
Strat_Series	2	Erosion	True	False	False

#get grid geometry as array
geo_model.grid

Grid Object. Values: 
array([[  20. ,   20. ,    7.5],
       [  20. ,   20. ,   22.5],
       [  20. ,   20. ,   37.5],
       ...,
       [1980. , 1980. ,  712.5],
       [1980. , 1980. ,  727.5],
       [1980. , 1980. ,  742.5]])

#plot surface points and orientations
plot = gp.plot_2d(geo_model, show_lith=False, show_boundaries=False)
plt.show()

#set interpolator
gp.set_interpolator(geo_model, compile_theano=True, theano_optimizer='fast_compile',)

Setting kriging parameters to their default values.
Compiling theano function...
Level of Optimization:  fast_compile
Device:  cpu
Precision:  float64
Number of faults:  1
Compilation Done!
Kriging values: 
                   values
range            2926.17
$C_o$             203869
drift equations   [3, 3]





<gempy.core.interpolator.InterpolatorModel at 0x7f61c11db110>

#get variables of interpolation
gp.get_data(geo_model, 'kriging')

	values
range	2926.17
$C_o$	203869
drift equations	[3, 3]

#perform interpolation
sol = gp.compute_model(geo_model)

#show interpolated geological model 
gp.plot_2d(geo_model, show_data=True,figsize=(12,6))
plt.tight_layout()
plt.show()

Blog

Basic tutorial of geological modeling with Python and Gempy