Model the Earth
with Python

# By default every series erodes the one below it.
# To make a series onlap instead (drape rather than cut):
geo_model.set_bottom_relation('Younger_Series', 'Onlap')

# Option A: conda (recommended)
conda create -n gempy_env python=3.10
conda activate gempy_env
pip install gempy

# Option B: pip (plain venv)
python -m venv .venv
source .venv/bin/activate   # Windows: .venv\Scripts\activate
pip install gempy

# Optional but useful for visualisation
pip install pyvista panel matplotlib

# Verify installation
import gempy as gp
print(gp.__version__)

# Expected: 2.x.x

import gempy as gp
import numpy as np

# ── Step 1: Create a GeoModel ──────────────────────────────
# The extent defines the bounding box of your model (x_min, x_max,
# y_min, y_max, z_min, z_max) in metres (or any consistent unit)
geo_model = gp.create_model('My_First_Model')

gp.init_data(
    geo_model,
    extent=[0, 2000,   # X: 2 km wide
            0, 2000,   # Y: 2 km deep
            0, 1500],  # Z: 1.5 km tall
    resolution=[50, 50, 50]  # 3-D grid cells for interpolation
)

# ── Step 2: Define surfaces (geological units) ─────────────
# These are the names of layers from TOP to BOTTOM
# (youngest to oldest, following superposition law)
gp.surfaces.add_surfaces(
    geo_model,
    ['Topsoil', 'Limestone', 'Shale', 'Sandstone']
)

# ── Step 3: Assign surfaces to series ──────────────────────
# A 'series' groups surfaces that were deposited together
# (same tectonic event / same time period)
gp.map_series_to_surfaces(
    geo_model,
    {
        'Sedimentary_Series': ('Topsoil', 'Limestone', 'Shale'),
        'Basement_Series':    ('Sandstone',)
    }
)

print(geo_model.surfaces)

# Surface points: WHERE does a geological surface pass through?
# Each row: X,  Y,    Z,    surface_name
# Z is elevation (positive = above a datum, usually sea level)

gp.add_surface_points(
    geo_model,
    X = [500, 1000, 1500,   # Limestone top: 3 points
         500, 1000, 1500,   # Shale top: 3 points
         500, 1000, 1500],  # Sandstone top: 3 points
    Y = [1000] * 9,          # All at Y=1000 (mid-model cross-section)
    Z = [1100, 1050, 1000,   # Limestone dips slightly east
         900,  870,  840,   # Shale
         600,  600,  600],  # Sandstone: flat
    surface = [
        'Limestone', 'Limestone', 'Limestone',
        'Shale',     'Shale',     'Shale',
        'Sandstone', 'Sandstone', 'Sandstone'
    ]
)

# Orientations: HOW is the surface tilted at a point?
# G_x, G_y, G_z = components of the surface NORMAL vector
# A flat horizontal surface has normal (0, 0, 1) — pointing straight up
# A surface dipping east has a normal tilted toward west

gp.add_orientations(
    geo_model,
    X       = [1000, 1000],
    Y       = [1000, 1000],
    Z       = [1050,  870],
    surface = ['Limestone', 'Shale'],
    pole_vector = [
        (0.087, 0, 0.996),  # Limestone: 5° east dip  (sin5°≈0.087, cos5°≈0.996)
        (0.087, 0, 0.996),  # Shale: same dip
    ]
)

# ── Tip: convert dip angle → pole vector ───────────────────
# import numpy as np
# dip_deg = 15   # degrees
# pole = (np.sin(np.radians(dip_deg)), 0, np.cos(np.radians(dip_deg)))

# ── Step 1: Set interpolation options ──────────────────────
gp.set_interpolator(
    geo_model,
    theano_optimizer='fast_run',  # 'fast_compile' for first run
    verbose=[]
)

# ── Step 2: Compute ────────────────────────────────────────
# This runs the kriging interpolation on the input data.
# Returns a Solution object containing scalar fields & surfaces
sol = gp.compute_model(geo_model)

print(sol)
# Solution: lith_block  — integer label per voxel (which layer is it?)
#           scalar_field_matrix — continuous scalar field values
#           vertices/edges — the extracted surface meshes

# ── Access the lithology block ──────────────────────────────
# lith_block[i] = layer index at voxel i
lith_block = sol.lith_block
print(lith_block.shape)  # (50*50*50,) → flatten of 50³ voxels

# Reshape to 3-D grid
lith_3d = lith_block.reshape(geo_model.grid.resolution)

# ── 2-D Cross-section (matplotlib) ─────────────────────────
gp.plot_2d(
    geo_model,
    direction='y',      # cross-section perpendicular to Y axis
    section_names=['y'],
    show_data=True,     # show surface points & orientations
    show_lith=True      # colour-fill by lithology
)

# ── 3-D interactive render (PyVista) ───────────────────────
p = gp.plot_3d(
    geo_model,
    show_surfaces=True,
    show_data=True,
    show_topography=False,
    image=False           # set True for a static PNG
)
# In Jupyter: this opens an interactive PyVista window

# ── Export to VTK for external viewers ─────────────────────
gp.plot_3d(geo_model, image=True)   # renders PNG inline in notebook

# Export surfaces as mesh files
for surf in geo_model.solutions.meshes:
    surf.mesh.save(f"output_{surf.surface}.vtk")

# ── 1. Add a fault surface ─────────────────────────────────
gp.surfaces.add_surfaces(geo_model, ['MainFault'])

# ── 2. Create a Fault series ───────────────────────────────
# Fault series MUST come BEFORE the geological series they cut
gp.map_series_to_surfaces(
    geo_model,
    {
        'Fault_Series':       ('MainFault',),   # <-- first
        'Sedimentary_Series': ('Limestone', 'Shale'),
        'Basement_Series':    ('Sandstone',)
    },
    remove_unused_series=True
)

# ── 3. Mark the series as a fault ─────────────────────────
geo_model.set_is_fault(['Fault_Series'])

# ── 4. Add fault observations ──────────────────────────────
gp.add_surface_points(
    geo_model,
    X=[1000, 1000], Y=[500, 1500], Z=[1200, 800],
    surface=['MainFault', 'MainFault']
)
gp.add_orientations(
    geo_model,
    X=[1000], Y=[1000], Z=[1000],
    surface=['MainFault'],
    pole_vector=[(0.7, 0, 0.7)]   # 45° dip
)

sol = gp.compute_model(geo_model)

# ── Option A: synthetic random topography ─────────────────
geo_model.set_topography(source='random', fd=1.9, d_z=200)
# fd  = fractal dimension (roughness): 1.5–2.5
# d_z = elevation range in model units

# ── Option B: from a GeoTIFF DEM ──────────────────────────
geo_model.set_topography(
    source='gdal',
    filepath='path/to/dem.tif'
)

# ── Option C: from a numpy array ──────────────────────────
# topo_array shape: (res_x, res_y)  values = elevation (Z)
topo_array = np.load('topo.npy')
geo_model.set_topography(source='numpy', array=topo_array)

# After setting topography, re-compute
sol = gp.compute_model(geo_model)

# Visualise with topo overlay
gp.plot_3d(geo_model, show_topography=True)

import gempy as gp
import numpy as np

# ── 1. CREATE ──────────────────────────────────────────────
geo_model = gp.create_model('ModelName')
gp.init_data(geo_model,
    extent=[x_min, x_max, y_min, y_max, z_min, z_max],
    resolution=[nx, ny, nz])

# ── 2. DEFINE STRUCTURE ────────────────────────────────────
gp.surfaces.add_surfaces(geo_model, ['S1', 'S2', 'S3'])  # youngest→oldest
gp.map_series_to_surfaces(geo_model, {'Series1': ('S1', 'S2'), 'Base': ('S3',)})

# ── 3. ADD DATA ────────────────────────────────────────────
gp.add_surface_points(geo_model, X=[...], Y=[...], Z=[...], surface=[...])
gp.add_orientations(geo_model, X=[...], Y=[...], Z=[...],
                   surface=[...], pole_vector=[...])

# ── 4. (OPTIONAL) FAULTS ──────────────────────────────────
geo_model.set_is_fault(['FaultSeries'])

# ── 5. (OPTIONAL) TOPOGRAPHY ──────────────────────────────
geo_model.set_topography(source='random')

# ── 6. COMPUTE ─────────────────────────────────────────────
gp.set_interpolator(geo_model)
sol = gp.compute_model(geo_model)

# ── 7. VISUALISE ───────────────────────────────────────────
gp.plot_2d(geo_model, direction='y', show_lith=True, show_data=True)
gp.plot_3d(geo_model, show_surfaces=True, show_topography=True)

# ── 8. INSPECT ─────────────────────────────────────────────
sol.lith_block           # voxel lithology IDs
sol.scalar_field_matrix  # raw scalar field
geo_model.solutions.meshes  # PyVista mesh objects

import gempy as gp

geo_model = gp.create_model('Beginner_FlatLayers')
gp.init_data(geo_model,
    extent=[0,1000,0,1000,0,800],
    resolution=[30,30,30])

gp.surfaces.add_surfaces(geo_model, ['Sand', 'Clay'])
gp.map_series_to_surfaces(geo_model, {'Strat': ('Sand', 'Clay')})

gp.add_surface_points(geo_model,
    X=[100,900], Y=[500,500], Z=[400,400],
    surface=['Sand','Sand'])
gp.add_orientations(geo_model,
    X=[500], Y=[500], Z=[400],
    surface=['Sand'], pole_vector=[(0,0,1)])  # flat → normal points straight up

gp.set_interpolator(geo_model)
sol = gp.compute_model(geo_model)
gp.plot_2d(geo_model, show_lith=True)

import numpy as np

dip_deg = 12
pole = (np.sin(np.radians(dip_deg)), 0, np.cos(np.radians(dip_deg)))

gp.add_orientations(geo_model,
    X=[500], Y=[500], Z=[450],
    surface=['Sand'], pole_vector=[pole])
# Re-run compute_model + plot_2d to see the tilt

geo_model = gp.create_model('Intermediate_Unconformity')
gp.init_data(geo_model, extent=[0,2000,0,2000,0,1000], resolution=[50,50,50])

gp.surfaces.add_surfaces(geo_model,
    ['Recent_Sand', 'Old_Limestone', 'Old_Shale', 'Basement'])

# Two series = two independent deformation events separated by an erosion surface
gp.map_series_to_surfaces(geo_model, {
    'Younger_Series': ('Recent_Sand',),
    'Older_Series':   ('Old_Limestone', 'Old_Shale'),
    'Basement_Series': ('Basement',)
})

# ... add surface_points + orientations for each surface ...
# Tip: give 'Older_Series' a steeper dip than 'Younger_Series'
# to visually show the angular unconformity

gp.surfaces.add_surfaces(geo_model, ['NormalFault'])
gp.map_series_to_surfaces(geo_model, {
    'Fault_Series': ('NormalFault',),
    'Strat_Series': ('Limestone', 'Shale', 'Sandstone')
}, remove_unused_series=True)
geo_model.set_is_fault(['Fault_Series'])
# ... add fault surface_points + orientations, then compute

geo_model.set_topography(source='gdal', filepath='region_dem.tif')
sol = gp.compute_model(geo_model)

# Get the geological map (which unit outcrops at the surface)
geo_map = sol.geological_map
gp.plot_2d(geo_model, show_topography=True, show_lith=True)
gp.plot_3d(geo_model, show_topography=True, show_surfaces=True)

gp.surfaces.add_surfaces(geo_model, ['Fault_A', 'Fault_B'])
gp.map_series_to_surfaces(geo_model, {
    'FaultA_Series': ('Fault_A',),   # older fault — cut first
    'FaultB_Series': ('Fault_B',),   # younger fault — cuts Fault_A too
    'Strat_Series':  ('Limestone', 'Shale')
})
geo_model.set_is_fault(['FaultA_Series', 'FaultB_Series'])

# Define that Fault_B offsets Fault_A (fault-fault relation matrix)
geo_model.set_fault_relation(
    np.array([[0,1,1],
              [0,0,1],
              [0,0,0]])
)

import pymc as pm
import gempy as gp

with pm.Model() as model:
    # Prior: true Z of a borehole pick is normally distributed
    # around the measured value (accounting for survey error)
    z_uncertain = pm.Normal('limestone_top_z', mu=1050, sigma=15)

    def run_gempy(z_val):
        geo_model.modify_surface_points(
            0, Z=z_val)  # update point 0's Z
        sol = gp.compute_model(geo_model)
        return sol.lith_block

    # Wrap as a PyTensor Op for gradient-free sampling (or use
    # gempy's built-in pymc backend in v3 for full HMC support)
    likelihood = pm.Potential(
        'observed_thickness',
        custom_log_likelihood(z_uncertain, observed_well_data)
    )

    trace = pm.sample(1000, tune=500, chains=4)

# Result: a distribution of plausible geological models,
# not just one "best guess" — quantifies subsurface uncertainty

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

# Map each lithology ID to a physical property
perm_lookup = {1: 50e-15, 2: 5e-15, 3: 200e-15}  # m²
permeability_field = np.vectorize(perm_lookup.get)(lith_3d)

# Export as VTK structured grid for a reservoir simulator
import pyvista as pv
grid = pv.UniformGrid()
grid.dimensions = np.array(permeability_field.shape) + 1
grid.cell_data['permeability'] = permeability_field.flatten(order='F')
grid.save('reservoir_perm.vtk')

# Now consumable by simulators like MRST, OPM, or custom solvers