gempy

GemPy - 3D Geological Modelling

Quick Reference

import gempy as gp

# Create model
geo_model = gp.create_geomodel(
    project_name='Model',
    extent=[0, 1000, 0, 1000, 0, 500],  # [xmin, xmax, ymin, ymax, zmin, zmax]
    resolution=[50, 50, 25]
)

# Add surface points and orientations
gp.add_surface_points(geo_model, x=[100, 500, 900], y=[500, 500, 500],
                      z=[400, 350, 400], surface='TopFormation')
gp.add_orientations(geo_model, x=[500], y=[500], z=[375],
                    pole_vector=[0, 0, 1], surface='TopFormation')

# Compute and visualize
gp.set_interpolator(geo_model)
sol = gp.compute_model(geo_model)
gp.plot_2d(geo_model, cell_number=[25], direction='y')

Key Classes

Class	Purpose
GeoModel	Main container - holds all model data
StructuralFrame	Manages geological relationships
Grid	Computation mesh for interpolation
Solutions	Model results (lithology, scalar fields)

Essential Operations

Define Geological Relationships

gp.map_stack_to_surfaces(geo_model, mapping={
    'Strata1': ['TopFormation', 'BaseFormation'],
    'Basement': ['Basement']
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.ERODE

Add Faults

gp.add_surface_points(geo_model, x=[300, 300, 300], y=[200, 500, 800],
                      z=[100, 250, 400], surface='Fault1')
gp.add_orientations(geo_model, x=[300], y=[500], z=[250],
                    pole_vector=[1, 0, 0.5], surface='Fault1')
gp.map_stack_to_surfaces(geo_model, mapping={
    'Fault_Series': ['Fault1'],
    'Strata1': ['Layer1', 'Layer2'],
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.FAULT

Load Data from Files

import pandas as pd

points_df = pd.read_csv('surface_points.csv')  # Columns: X, Y, Z, surface
gp.add_surface_points(geo_model, x=points_df['X'], y=points_df['Y'],
                      z=points_df['Z'], surface=points_df['surface'])

ori_df = pd.read_csv('orientations.csv')  # Columns: X, Y, Z, dip, azimuth, surface
gp.add_orientations(geo_model, x=ori_df['X'], y=ori_df['Y'], z=ori_df['Z'],
                    dip=ori_df['dip'], azimuth=ori_df['azimuth'],
                    surface=ori_df['surface'])

Visualization

gp.plot_2d(geo_model, cell_number=[25], direction='y', show_data=True)
gp.plot_3d(geo_model, show_data=True, show_surfaces=True)  # Requires PyVista

Access Results

sol = gp.compute_model(geo_model)
lithology = sol.raw_arrays.lith_block
lith_3d = lithology.reshape(geo_model.grid.regular_grid.resolution)

Structural Relation Types

Type	Description
ERODE	Younger surface erodes older (unconformity)
ONLAP	Younger onlaps onto older
FAULT	Surface is a fault plane
INTRUSION	Intrusive body

When to Use vs Alternatives

Scenario	Recommendation
3D implicit geological modelling from surface data	GemPy - purpose-built, Python-native
Complex fold modelling with structural frames	LoopStructural - better fold support
Commercial-grade subsurface modelling	SKUA-GOCAD - industry standard, proprietary
Quick 2D cross-sections only	GemPy works but may be overkill; consider manual interpolation
Need gravity/magnetics forward model from geology	GemPy - has built-in potential field support

Choose GemPy when: You need implicit 3D geological modelling from surface contacts and orientations with fault/unconformity relationships, especially when integrated with Python workflows. It has a gentler learning curve than LoopStructural for standard cases.

Avoid GemPy when: Your geology is dominated by complex folding (use LoopStructural), or you need production-grade reservoir modelling (use commercial tools).

Common Workflows

Build 3D geological model from surface data

• Load or define surface contact points (min 2 per surface) and orientations (min 1 per surface)
• Create GeoModel with appropriate extent and resolution
• Add surface points with gp.add_surface_points()
• Add orientations with gp.add_orientations()
• Map geological stack with gp.map_stack_to_surfaces()
• Set structural relations (ERODE, ONLAP, FAULT)
• Set interpolator with gp.set_interpolator()
• Compute model with gp.compute_model()
• Validate with 2D cross-sections using gp.plot_2d()
• Visualize 3D result with gp.plot_3d() or export to VTK

Common Issues

Issue	Solution
Model not computing	Check min 2 points + 1 orientation per surface
Artifacts at edges	Extend model extent beyond data
Wrong fault offset	Check pole_vector direction
Memory errors	Reduce grid resolution

References

• Structural Relations - Detailed guide on ERODE, ONLAP, FAULT
• Data Requirements - Input data formats and constraints
• Troubleshooting - Common problems and solutions

Scripts

• scripts/validate_input.py - Validate input data files
• scripts/export_model.py - Export model to VTK, numpy, or CSV