steadfastasart/geoscience-skills/geological-modelling

gemgis              gempy / loopstructural       pyvista
[GIS Preprocessing] --> [Implicit Modelling]     --> [3D Visualization]
  |                      |                           |
  Shapefile parsing      Surface interpolation       Volume render
  Raster extraction      Fault modelling             Cross-sections
  Borehole to points     Unconformities              Mesh export
  CRS transforms         Scalar field solving        Interactive pick

import gemgis as gg
import geopandas as gpd
import rasterio
# Load geological map (shapefile)
contacts = gpd.read_file('geological_contacts.shp')
orientations = gpd.read_file('orientations.shp')
# Extract surface points from GIS contacts with DEM
with rasterio.open('dem.tif') as dem:
    # All inputs must share the DEM's projected CRS and metre units.
    contacts = contacts.to_crs(dem.crs)
    orientations = orientations.to_crs(dem.crs)
    surface_points = gg.vector.extract_xyz(contacts, dem=dem)
    orientation_pts = gg.vector.extract_xyz(orientations, dem=dem)
# Define model extent from data bounds
extent = gg.utils.set_extent(
    gdf=surface_points,
    z_min=-500, z_max=1000
)

# example: gempy-model
import numpy as np
import gempy as gp
from gempy.core.data.options import InterpolationOptionsType
# Synthetic contact: z = -227.5 - 0.125 * (x - 1000).
extent = np.array([1000., 1200., 2000., 2120., -300., -180.])
resolution = np.array([8, 6, 5])  # Number of CELLS, not corner points.
surface_table = gp.data.SurfacePointsTable.from_arrays(
    x=np.array([1020., 1180., 1020.]),
    y=np.array([2020., 2020., 2100.]),
    z=np.array([-230., -250., -230.]),
    names='Sandstone_base',
)
pole = np.array([0.125, 0., 1.])
pole /= np.linalg.norm(pole)
orientation_table = gp.data.OrientationsTable.from_arrays(
    x=np.array([1060.]), y=np.array([2060.]), z=np.array([-235.]),
    G_x=pole[[0]], G_y=pole[[1]], G_z=pole[[2]],
    names=['Sandstone_base'], name_id_map=surface_table.name_id_map,
)
frame = gp.data.StructuralFrame.from_data_tables(surface_table, orientation_table)
geo_model = gp.create_geomodel(
    project_name='SyntheticInclinedContact',
    extent=extent, resolution=resolution, structural_frame=frame,
    intpolation_options_tye=InterpolationOptionsType.DENSE_GRID,
)
gp.map_stack_to_surfaces(geo_model, mapping_object={'Stratigraphy': ['Sandstone_base']})
geo_model.interpolation_options.evaluation_options.mesh_extraction = False
sol = gp.compute_model(
    geo_model,
    engine_config=gp.data.GemPyEngineConfig(backend=gp.data.AvailableBackends.numpy),
)

from LoopStructural import GeologicalModel
import pandas as pd
# Prepare input DataFrames
data = pd.DataFrame({
    'X': x_coords, 'Y': y_coords, 'Z': z_coords,
    'feature_name': formation_names,
    'val': stratigraphic_values,     # Scalar values for interface position
    'gx': gradient_x, 'gy': gradient_y, 'gz': gradient_z  # Orientation
})
# Build model
model = GeologicalModel(
    origin=[extent[0], extent[2], extent[4]],
    maximum=[extent[1], extent[3], extent[5]]
)
model.data = data
# Add features
model.create_and_add_foliation('Stratigraphy', interpolatortype='FDI')
model.create_and_add_fault('MainFault', displacement=100)
# Update and access
model.update()
lithology = model.evaluate_model(model.regular_grid())

# example: gempy-export
import numpy as np
import pyvista as pv
regular = geo_model.grid.regular_grid
shape = np.asarray(regular.resolution, dtype=int)
bounds = np.asarray(regular.extent, dtype=float)
lith_block = np.asarray(sol.raw_arrays.lith_block)
if lith_block.size != np.prod(shape):
    raise ValueError('Export requires a complete dense-grid lithology block')
spacing = (bounds[1::2] - bounds[::2]) / shape
grid = pv.ImageData(
    dimensions=tuple(shape + 1), spacing=tuple(spacing), origin=tuple(bounds[::2]),
)
grid.cell_data['lithology'] = lith_block.reshape(tuple(shape), order='C').ravel(order='F')
grid.save('geological_model.vti')

import pyvista as pv
plotter = pv.Plotter()
sliced = grid.slice(normal='y', origin=grid.center)
plotter.add_mesh(sliced, scalars='lithology', cmap='tab10', categories=True)
plotter.add_mesh(grid.outline(), color='black')
plotter.show_axes()
plotter.show()

- [ ] Gather input data: geological map, DEM, borehole logs, structural measurements
- [ ] Load shapefiles and rasters with geopandas and rasterio
- [ ] Extract surface contact points with elevation using gemgis
- [ ] Extract orientation data (dip, azimuth) with gemgis
- [ ] Define model extent and resolution (cover data + buffer)
- [ ] Create GemPy GeoModel with extent and resolution
- [ ] Add surface points and orientations
- [ ] Define stratigraphic pile and structural relationships (ERODE, ONLAP)
- [ ] Add fault surfaces if present
- [ ] Configure interpolation options and compute model
- [ ] Validate with 2D cross-sections through known data points
- [ ] Visualize 3D result with pyvista
- [ ] Export to VTK or numpy for downstream use

- [ ] Load borehole data (collar, survey, lithology intervals)
- [ ] Convert lithology picks to surface contact points at formation boundaries
- [ ] Estimate orientations from multi-well dip calculation or assign regional dip
- [ ] Supply sufficient contact geometry and at least one orientation per structural group
- [ ] Validate: check model honours borehole intersections
- [ ] Iterate: add more data or adjust orientations to fix artifacts

- [ ] Load geological map polygons and structural measurements from shapefiles
- [ ] Reproject to common CRS with geopandas
- [ ] Extract formation boundary polylines from polygon contacts
- [ ] Sample points along polylines with gemgis
- [ ] Drape points onto DEM to get 3D coordinates
- [ ] Build GemPy model from extracted points and orientations
- [ ] Compare model surface traces with original geological map