1.1 -Basics of geological modeling with GemPy

import numpy as np

Importing Necessary Libraries

import gempy as gp
import gempy_viewer as gpv

Importing and Defining Input Data

gempy.core.data.GeoModel GemPy uses Python objects to store the data that builds the geological model. The main data classes include:

You can also create data from raw CSV files (comma-separated values). This could be useful if you are exporting model data from a different program or creating it in a spreadsheet software like Microsoft Excel or LibreOffice Calc.

In this tutorial, we’ll use CSV files to generate input data. You can find these example files in the gempy data repository on GitHub. The data consists of x, y, and z positional values for all surface points and orientation measurements. Additional data includes poles, azimuth and polarity (or the gradient components). Surface points are assigned a formation, which can be a lithological unit (like “Sandstone”) or a structural feature (like “Main Fault”).

It’s important to note that, in GemPy, interface position points mark the bottom of a layer. If you need points to represent the top of a formation (for example, when modeling an intrusion), you can define an inverted orientation measurement.

While generating data from CSV files, we also need to define the model’s real extent in x, y, and z. This extent defines the area used for interpolation and many of the plotting functions. We also set a resolution to establish a regular grid right away. This resolution will dictate the number of voxels used during modeling. We’re using a medium resolution of 50x50x50 here, which results in 125,000 voxels. The model extent should enclose all relevant data in a representative space. As our model voxels are prisms rather than cubes, the resolution can differ from the extent. However, it is recommended to avoid going beyond 100 cells in each direction (1,000,000 voxels) to prevent excessive computational costs.

data_path = 'https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/'

geo_model: gp.data.GeoModel = gp.create_geomodel(
    project_name='Tutorial_ch1_1_Basics',
    extent=[0, 2000, 0, 2000, 0, 750],
    resolution=[20, 20, 20],  # * Here we define the resolution of the voxels
    refinement=4,  # * Here we define the number of octree levels. If octree levels are defined, the resolution is ignored.
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=data_path + "/data/input_data/getting_started/simple_fault_model_orientations.csv",
        path_to_surface_points=data_path + "/data/input_data/getting_started/simple_fault_model_points.csv",
        hash_surface_points="4cdd54cd510cf345a583610585f2206a2936a05faaae05595b61febfc0191563",
        hash_orientations="7ba1de060fc8df668d411d0207a326bc94a6cdca9f5fe2ed511fd4db6b3f3526"
    )
)
Surface points hash:  4cdd54cd510cf345a583610585f2206a2936a05faaae05595b61febfc0191563
Orientations hash:  7ba1de060fc8df668d411d0207a326bc94a6cdca9f5fe2ed511fd4db6b3f3526

New in GemPy 3!

GemPy 3 has introduced the ImporterHelper class to streamline importing data from various sources. This class simplifies the process of passing multiple arguments needed for importing data and will likely see further extensions in the future. Currently, one of its uses is to handle pooch arguments for downloading data from the internet.

# The input data can be reviewed using the properties `surface_points` and `orientations`. However, note that at this point,
# the sequence of formations and their assignment to series are still arbitrary. We will rectify this in the subsequent steps.
XYZidnugget
700.001000.00300.00869021620.00
600.001000.00200.00869021620.00
500.001000.00100.00869021620.00
1000.001000.00600.00869021620.00
1100.001000.00700.00869021620.00
1000.0050.00350.001751782670.00
1000.00150.00333.331751782670.00
1000.00300.00333.331751782670.00
1000.00500.00366.671751782670.00
1000.001000.00433.331751782670.00
...............
1100.001700.00473.334106404740.00
1100.001950.00490.004106404740.00
0.001000.00433.334106404740.00
300.001000.00400.004106404740.00
600.001000.00366.674106404740.00
1300.001000.00566.674106404740.00
1600.001000.00550.004106404740.00
1900.001000.00566.674106404740.00
1700.00500.00533.334106404740.00
1700.001500.00516.674106404740.00


XYZG_xG_yG_zidnugget
500.001000.00300.00-0.95-0.000.32869021620.01
400.001000.00420.000.320.000.952090083980.01
1000.001000.00300.000.320.000.953256836050.01


Declaring the Sequential Order of Geological Formations

In our model, we want the geological units to appear in the correct chronological order. Such order could be determined by a sequence of stratigraphic deposition, unconformities due to erosion, or other lithological genesis events like igneous intrusions. A similar age-related order is declared for faults in our model. In GemPy, we use the function gempy.map_stack_to_surfaces to assign formations or faults to different sequential series by declaring them in a Python dictionary.

The correct ordering of series is crucial for model construction! It’s possible to assign several surfaces to one series. The order of units within a series only affects the color code, so we recommend maintaining consistency. The order can be defined by simply changing the order of the lists within gempy.core.data.StructuralFrame.structural_groups and gempy.core.data.StructuralGroups.elements attributes.

Faults are treated as independent groups and must be younger than the groups they affect. The relative order between different faults defines their tectonic relationship (the first entry is the youngest).

For a model with simple sequential stratigraphy, all layer formations can be assigned to one series without an issue. All unit boundaries and their order would then be determined by interface points. However, to model more complex lithostratigraphical relations and interactions, separate series definition becomes important. For example, modeling an unconformity or an intrusion that disrupts older stratigraphy would require declaring a “newer” series.

By default, we create a simple sequence inferred from the data:

Structural Groups: StructuralGroup:
Name:default_formation
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Main_Fault

StructuralElement:
Name:Sandstone_1

StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Shale

StructuralElement:
Name:Siltstone
Fault Relations:
default_fo...
default_formation
True
False


Our example model comprises four main layers (plus an underlying basement that is automatically generated by GemPy) and one main normal fault displacing those layers. Assuming a simple stratigraphy where each younger unit was deposited onto the underlying older one, we can assign these layer formations to one series called “Strat_Series”. For the fault, we declare a respective “Fault_Series” as the first key entry in the mapping dictionary. We could give any other names to these series, the formations however have to be referred to as named in the input data.

gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object=  # TODO: This mapping I do not like it too much. We should be able to do it passing the data objects directly
    {
        "Fault_Series": 'Main_Fault',
        "Strat_Series": ('Sandstone_2', 'Siltstone', 'Shale', 'Sandstone_1')
    }
)

geo_model.structural_frame  # Display the resulting structural frame
Structural Groups: StructuralGroup:
Name:Fault_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Main_Fault

StructuralGroup:
Name:Strat_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Siltstone

StructuralElement:
Name:Shale

StructuralElement:
Name:Sandstone_1
Fault Relations:
Fault_Seri...Strat_Seri...
Fault_Series
Strat_Series
True
False


gp.set_is_fault(
    frame=geo_model.structural_frame,
    fault_groups=['Fault_Series']
)
Structural Groups: StructuralGroup:
Name:Fault_Series
Structural Relation:StackRelationType.FAULT
Elements:
StructuralElement:
Name:Main_Fault

StructuralGroup:
Name:Strat_Series
Structural Relation:StackRelationType.ERODE
Elements:
StructuralElement:
Name:Sandstone_2

StructuralElement:
Name:Siltstone

StructuralElement:
Name:Shale

StructuralElement:
Name:Sandstone_1
Fault Relations:
Fault_Seri...Strat_Seri...
Fault_Series
Strat_Series
True
False


Now, all surfaces have been assigned to a series and are displayed in the correct order (from young to old).

Returning Information from Our Input Data

Our model input data, named “geo_model”, contains all essential information for constructing our model. You can access different types of information by accessing the attributes. For instance, you can retrieve the coordinates of our modeling grid as follows:

<gempy.core.data.grid.Grid object at 0x7f131d34ff70>

Visualizing input data

We can also visualize our input data. This might for example be useful to check if all points and measurements are defined the way we want them to. Using the function gempy_viewer.plot2d, we attain a 2D projection of our data points onto a plane of chosen direction (we can choose this attribute to be either \(x\), \(y\) or \(z\)).

plot = gpv.plot_2d(geo_model, show_lith=False, show_boundaries=False)
Cell Number: mid Direction: y

Using gempy_viewer.plot_3d, # we can also visualize this data in 3D. Note that direct 3D visualization in GemPy requires the Visualization Toolkit (VTK) to be installed.

gpv.plot_3d(geo_model, image=False, plotter_type='basic')
ch1 1 basics
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f131d39d720>

Model Generation

Once we’ve correctly defined all our primary information in our gempy.core.data.GeoModel object (referred to as geo_model in these tutorials), we can proceed to the next step: preparing the input data for interpolation.

New in GemPy 3! Numpy and TensorFlow backend

Unlike previous versions, GemPy 3 doesn’t rely on theano or asera. Instead, it utilizes numpy or tensorflow. Consequently, we no longer need to recompile all theano functions (TensorFlow uses eager execution; we found no notable speed difference after profiling the XLA compiler).

The parameters used for the interpolation are stored in gempy.core.data.GeoModel.interpolation_options. These parameters have sensible default values that you can modify if necessary. However, we advise caution when changing these parameters unless you fully understand their implications.

Display the current interpolation options

geo_model.interpolation_options
InterpolationOptions
kernel_options{'range': 5, 'c_o': 10, 'uni_degree': 1, 'i_res': 4, 'gi_res': 2, 'number_dimensions': 3, 'kernel_function': , derivative_div_r=, second_derivative=, consume_sq_distance=False)>, 'compute_condition_number': False, 'kernel_solver': }
number_octree_levels4
current_octree_level0
compute_scalar_gradientFalse
mesh_extractionTrue
mesh_extraction_masking_optionsMeshExtractionMaskingOptions.INTERSECT
mesh_extraction_fancyTrue
debugTrue
debug_water_tightFalse
sigmoid_slope50000
_number_octree_levels_surface4


With all our prerequisites in place, we can now compute our complete geological model using gempy.compute_model(). This function returns a gempy.core.data.Solutions object.

The following sections illustrate these different model solutions and how to utilize them.

Compute the geological model and get the solutions

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(
Solutions: 4 Octree Levels, 5 DualContouringMeshes


Solutions are also stored within the gempy.core.data.GeoModel object, for future reference.

geo_model.solutions
Solutions: 4 Octree Levels, 5 DualContouringMeshes


Direct model visualization in GemPy

Model solutions can be easily visualized in 2D sections in GemPy directly. Let’s take a look at our lithology block:

gpv.plot_2d(geo_model, show_data=True, cell_number="mid", direction='y')
Cell Number: mid Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f12b2c80400>

With cell_number=mid, we have chosen a section going through the middle of our block. We have moved in direction='y', the plot thus depicts a plane parallel to the \(x\)- and \(y\)-axes. Setting show_data=True, we could plot original data together with the results. Changing the values for cell_number and direction, we can move through our 3D block model and explore it by looking at different 2D planes.

We can do the same with the underlying scalar-field solution:

gpv.plot_2d(
    model=geo_model,
    series_n=0,  # This will plot the scalar field used for the fault
    show_data=False,
    show_scalar=True,
    show_lith=False
)
Cell Number: mid Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f1326e820b0>
gpv.plot_2d(
    model=geo_model,
    series_n=1,  # This will plot the scalar field used for the stratigraphy
    show_data=False,
    show_scalar=True,
    show_lith=False
)
Cell Number: mid Direction: y
<gempy_viewer.modules.plot_2d.visualization_2d.Plot2D object at 0x7f131ad60970>

Dual Contouring and vtk visualization

In addition to 2D sections we can extract surfaces to visualize in 3D renderers. Surfaces can be visualized as 3D triangle complexes in VTK (see function plot_surfaces_3D below). To create these triangles, we need to extract respective vertices and simplices from the potential fields of lithologies and faults. This process is automatized in GemPy using dual contouring in the gempy_engine.

New in GemPy 3! Dual Contouring

GemPy 3 uses dual contouring to extract surfaces from the scalar fields. The method is completely coded in gempy_engine what also enables further improvements in the midterm. This method is more efficient to use together with octrees and suited better the new capabilities of gempy3.

gpv.plot_3d(geo_model, show_data=False, image=False, plotter_type='basic')
ch1 1 basics
<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f1326e820b0>

Adding topography

In gempy we can add more grid types for different purposes. We will explore this concept in more detail in the next tutorials (1.3a: Grids.). For now, we will just add a topography grid to our model. This grid allows us to intersect the surfaces as well as compute a high resolution geological mal.

gp.set_topography_from_random(
    grid=geo_model.grid,
    fractal_dimension=1.2,
    d_z=np.array([350, 750]),
    topography_resolution=np.array([50, 50]),
)

gp.compute_model(geo_model)
gpv.plot_2d(geo_model, show_topography=True)

gpv.plot_3d(
    model=geo_model,
    plotter_type='basic',
    show_topography=True,
    show_surfaces=True,
    show_lith=True,
    image=False
)
ch1 1 basicsCell Number: mid Direction: y
Active grids: ['topography']
Setting Backend To: AvailableBackends.numpy

<gempy_viewer.modules.plot_3d.vista.GemPyToVista object at 0x7f12d5852950>

Compute at a given location

This is done by modifying the grid to a custom grid and recomputing.

x_i = np.array([[1000, 1000, 1000]])
lith_values_at_coords: np.ndarray = gp.compute_model_at(
    gempy_model=geo_model,
    at=x_i
)
lith_values_at_coords
Active grids: ['custom' 'topography']
Setting Backend To: AvailableBackends.numpy

array([2.])

Therefore if we just want the value at x_i:

geo_model.solutions.raw_arrays.custom
array([2.])

Work in progress

GemPy3 model serialization is currently being redisigned. Therefore, at the current version, there is not a build in method to save the model. However, since now the data model should be completely robust, you should be able to save the gempy.core.data.GeoModel and all its attributes using the standard python library [pickle](https://docs.python.org/3/library/pickle.html)

sphinx_gallery_thumbnail_number = -2

Total running time of the script: (0 minutes 6.253 seconds)

Gallery generated by Sphinx-Gallery