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:
    surface_points = gg.vector.extract_xyz(contacts, dem=dem)

# Extract orientations with elevation
orientation_pts = gg.vector.extract_xyz(orientations, dem=dem)

# Extract borehole data
boreholes = gpd.read_file('boreholes.shp')
bh_points = gg.vector.extract_xyz_from_cross_sections(boreholes)

# Define model extent from data bounds
extent = gg.utils.set_extent(
    gdf=surface_points,
    z_min=-500, z_max=1000
)

import gempy as gp

# Create model
geo_model = gp.create_geomodel(
    project_name='RegionalModel',
    extent=extent,
    resolution=[50, 50, 30]
)

# Add data from GemGIS output
gp.add_surface_points(
    geo_model,
    x=surface_points['X'], y=surface_points['Y'],
    z=surface_points['Z'], surface=surface_points['formation']
)
gp.add_orientations(
    geo_model,
    x=orientation_pts['X'], y=orientation_pts['Y'],
    z=orientation_pts['Z'],
    dip=orientation_pts['dip'],
    azimuth=orientation_pts['azimuth'],
    surface=orientation_pts['formation']
)

# Define stratigraphic relationships
gp.map_stack_to_surfaces(geo_model, mapping={
    'Cover': ['Alluvium'],
    'Sedimentary': ['Sandstone', 'Shale', 'Limestone'],
    'Basement': ['Granite']
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.ERODE

# Add faults
gp.map_stack_to_surfaces(geo_model, mapping={
    'Fault_Series': ['MainFault'],
    'Cover': ['Alluvium'],
    'Sedimentary': ['Sandstone', 'Shale', 'Limestone'],
    'Basement': ['Granite']
})
geo_model.structural_frame.structural_groups[0].structural_relation = \
    gp.data.StackRelationType.FAULT

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

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())

import pyvista as pv

# GemPy model to PyVista
lith_block = sol.raw_arrays.lith_block
grid_3d = lith_block.reshape(geo_model.grid.regular_grid.resolution)

grid = pv.ImageData(dimensions=geo_model.grid.regular_grid.resolution)
grid.point_data['lithology'] = lith_block.flatten(order='F')

# 3D visualization with cross-section
plotter = pv.Plotter()
plotter.add_volume(grid, scalars='lithology', cmap='tab10', opacity='sigmoid')
sliced = grid.slice(normal='y', origin=grid.center)
plotter.add_mesh(sliced, scalars='lithology', cmap='tab10')
plotter.show()

# Export to VTK for external tools
grid.save('geological_model.vtk')

- [ ] 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
- [ ] Set interpolator 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
- [ ] Build model with GemPy (minimum 2 points + 1 orientation per surface)
- [ ] 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