1. Introduction to GemPy

GemPy is an open-source Python library for 3D structural geological modeling. It uses an implicit modeling approach based on universal cokriging to construct 3D scalar fields, enabling fast creation of complex models including folds, faults, and unconformities. It also natively supports probabilistic modeling and uncertainty analysis.

2. Minimum Required Datasets

To build a geological model in GemPy, two core input data types are required. These can be extracted from borehole logs, geological cross-sections, and outcrop measurements.

2.1 Surface Points

Surface points are 3D coordinates marking geological interfaces (e.g., layer boundaries, fault planes, unconformities). In GemPy, each surface represents the bottom of the unit it is assigned to.

2.2 Orientation Measurements

Orientation data describes the spatial dip of geological interfaces, constraining the geometry of interpolation. At least one orientation per structural group is required.

3. Extracting Data from Borehole Logs

Borehole logs are the most common source of data for geological modeling. Below is the complete workflow to convert borehole data into GemPy inputs:

Step 1: Collect Raw Borehole Materials

The user should provide the following:

1. Borehole location table: Borehole ID, collar coordinates (X, Y), collar elevation (Z), total depth, inclination and azimuth (for deviated holes)
2. Stratigraphic interval table: Bottom depth of each layer, bottom elevation, formation code, formation name, lithology description
3. Structural orientation data: Dip direction and dip angle measured within or near the borehole
4. Study area extent: X/Y/Z minimum and maximum values defining the modeling domain
5. Stratigraphic sequence: Depositional order, presence of unconformities or intrusions, fault information

Step 2: Compute Surface Point Coordinates

For each layer bottom in each borehole, calculate the 3D coordinates:

# Vertical borehole example: collar (X0, Y0, Z0), layer bottom depth D
X = X0
Y = Y0
Z = Z0 - D   # if Z is negative downward

# Deviated borehole conversion using zenith angle theta and azimuth alpha:
# X = X0 + D * sin(theta) * sin(alpha)
# Y = Y0 + D * sin(theta) * cos(alpha)
# Z = Z0 - D * cos(theta)

Step 3: Prepare CSV Files

Save extracted surface points as surface_points.csv:

X,Y,Z,formation
0,0,800,rock3
0,0,450,rock2
0,0,250,rock1
500,0,800,rock3
...

Save orientation data as orientations.csv:

X,Y,Z,azimuth,dip,polarity,formation
100,500,650,90,315,1,rock2
100,500,450,90,315,1,rock1
...

4. Basic GemPy Modeling Workflow (Demo)

4.1 Install GemPy

# Core library
pip install gempy

# Visualization module
pip install gempy-viewer

# 3D visualization dependency (recommended)
pip install pyvista

# For GPU acceleration (PyTorch backend)
pip install torch

4.2 Import Data from CSV and Create Model

import gempy as gp
import gempy_viewer as gpv

# Data paths
data_path = "https://raw.githubusercontent.com/cgre-aachen/gempy_data/master/"
path_to_data = data_path + "/data/input_data/video_tutorials_v3/"

# Create geological model
geo_model = gp.create_geomodel(
    project_name='tutorial_model',
    extent=[0, 2500, 0, 1000, 0, 1000],
    resolution=[100, 40, 40],
    importer_helper=gp.data.ImporterHelper(
        path_to_orientations=path_to_data + "tutorial_model_orientations.csv",
        path_to_surface_points=path_to_data + "tutorial_model_surface_points.csv"
    )
)

4.3 Define Structural Sequence

# Map formations to structural groups (youngest to oldest)
gp.map_stack_to_surfaces(
    gempy_model=geo_model,
    mapping_object={
        "Strat_Series2": ("rock3"),
        "Strat_Series1": ("rock2", "rock1")
    }
)

4.4 Compute the Model

import os
os.environ["DEFAULT_BACKEND"] = "PYTORCH"  # or "NUMPY"

gp.compute_model(geo_model)

4.5 Visualize Results

# 2D cross-section
gpv.plot_2d(geo_model, cell_number=20, show_data=True)

# 3D interactive visualization
gpv.plot_3d(geo_model, show_lith=True, show_boundaries=True)

4.6 Create from Scratch (Without CSV)

# Create empty model
geo_model = gp.create_geomodel(
    project_name='Model1',
    extent=[0, 780, -200, 200, -582, 0],
    resolution=(50, 50, 50),
    structural_frame=gp.data.StructuralFrame.initialize_default_structure()
)

# Add surface points manually
gp.add_surface_points(
    geo_model=geo_model,
    x=[225, 460, 617],
    y=[0, 0, 0],
    z=[-95, -100, -10],
    elements_names=['Limestone', 'Limestone', 'Limestone']
)

# Add orientation manually
gp.add_orientations(
    geo_model=geo_model,
    x=[350], y=[0], z=[-120],
    elements_names=['Limestone'],
    pole_vector=[[0, 0, 1]]  # Horizontal bedding, normal vector up
)

# Compute and visualize
gp.compute_model(geo_model)
gpv.plot_2d(geo_model)

5. Complete Data Requirements Checklist

To successfully build a 3D geological model with GemPy, please ensure the following materials are provided:

6. Official Documentation & Help Links