gempy

name	gempy
description	3D structural geological modeling using implicit methods. Create geological models with faults, folds, and unconformities from surface points and orientations. Use when Claude needs to: (1) Build 3D geological models from surface contacts and orientations, (2) Model faults, unconformities, or intrusions, (3) Compute and visualize subsurface geology, (4) Export models to VTK or numpy arrays, (5) Generate gravity forward models, (6) Create cross-sections or 3D visualizations.
version	1.0.0
author	Geoscience Skills
license	MIT
tags	["Geological Modelling","3D","Implicit Surfaces","Faults","Structural Geology"]
dependencies	["gempy>=2.3.0","numpy","pyvista"]
complements	["gemgis","loopstructural","pyvista"]
workflow_role	modelling

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