Geological Modelling#

import geolime as geo
from pyproj import CRS
import numpy as np
import pyvista as pv

pv.set_jupyter_backend('panel')


geo.Project().set_crs(CRS("EPSG:20350"))
/tmp/ipykernel_2888/3544329244.py:6: PyVistaDeprecationWarning: `panel` backend is deprecated and is planned for future removal.
  pv.set_jupyter_backend('panel')
# Function to graphically project data on X or Y plane
def project_section(
    dh: geo.Drillholes,
    coord: geo.Coord,
    var: str,
    region: str = None,
    flatten = False
):
    Z = "Z_M"
    yaxis = True

    if flatten:
        Z = "FROM"
        yaxis = "reversed"
    
    f = geo.scatter(geo_object=dh, property_x=f"{coord}_M", property_y=Z, region=region, yaxis_autorange=yaxis)
    
    # color options
    f["data"][0]["marker"]["color"] = dh.data(var, region)
    f["data"][0]["marker"]["colorscale"] = "rainbow"
    f["data"][0]["marker"]["line"]["width"] = 0
    
    return f
# reload drillhole data
dh = geo.read_file("../data/dh_hyper.geo")

Create Channel Geometry#

# center of the ellipsoid defining the zone of interest
cx = 547895.7
cy = 7475196.5
cz = np.max(dh["Z_B"])

# ellispoid range geometry
rx = 1000
ry = 2400
rz = 75

# compute the distance of each drillhole point wrt to the enveloppe of the ellipsoid
dh.set_property_expr(
    name="ellipsoidal_distance",
    expr=f"-1. * ((X_M - {cx}) * (X_M - {cx}) / ({rx} * {rx}) + (Y_M - {cy}) * (Y_M - {cy}) / ({ry} * {ry}) + (Z_M - {cz}) * (Z_M - {cz}) / ({rz} * {rz}))"
)

# define Zone Of Interest (OreZone) as anything inside ellipsoid enveloppe
dh.set_property_expr(
    name="OreZone",
    expr="(ellipsoidal_distance > -1) * (-ellipsoidal_distance)"
)

# map it to QC areal extent of enveloppe
dh.plot_2d(property="ellipsoidal_distance", agg_method="mean", interactive_map=True)
Make this Notebook Trusted to load map: File -> Trust Notebook
project_section(dh=dh, coord="X", var="OreZone", flatten=True)
project_section(dh=dh, coord="X", var="Fe_pct", flatten=True)
project_section(dh=dh, coord="Y", var="OreZone", flatten=True)

Create Domain based on distance to channel#

# set domain codes as follows
# 1: outside OreZone
# 2: topo surface = drillholes collars
# 3: OreZone
dh.set_property_expr(name="domain_code", expr="1")
dh.set_property_expr(name="domain_code", expr="3", region="OreZone > 0")
dh.set_property_expr(name="domain_code", expr="2", region="FROM < 0.01")
project_section(dh=dh, coord="X", var="domain_code")
project_section(dh=dh, coord="Y", var="domain_code")
# translate domain code as string for human-readability
dh.set_property_expr(name="domain", expr="'OreZone'")
dh.set_property_expr(name="domain", expr="'Lower'", region="domain_code == 1")
dh.set_property_expr(name="domain", expr="'Upper'", region="domain_code == 2")

Surfaces Creation#

Creation of Gridded Surface for interpolation#

# define cell size for the surface
size_x = 50
size_y = 50

# create regular 2D Voxel covering the drillholes areal extent (bounding box)
bbox = dh.bounds()
vx = geo.Voxel(
    name="Grid2D",
    shape=np.round((bbox[1] - bbox[0]) / np.array([size_x, size_y, 1]))[:2].astype(int) + np.array([4, 4]),
    origin=bbox[0] - np.array([50, 50, 0]),
    axis=[
        [size_x, 0, 0],
        [0, size_y, 0]
    ]
)

gs = geo.GriddedSurface(name="geol_surfaces", xyz=vx.coords())

print(f"Gridded Surface Shape: {gs.shape}")
Gridded Surface Shape: [ 64 140]

Creation of lithological sequence#

lithos = ["Upper", "OreZone", "Lower"]

Estimate Surface#

estimator = geo.ClosestPointValue(
    method=geo.SurfaceEstimationMethod.THICKNESS,
    interpolator=geo.DrillholesInterpolationMethod.RBF
)
estimator.build_constraints(data=dh, domain_col="domain", domain_sequence=lithos)
estimator.fill_constraints()
estimator.surface_creation(grid_surf=gs, interpolator=geo.ExternalSurfaceInterpolationMethod.RBF)

Visual check of created surfaces#

# plot to QC
gs.plot_2d(property="hard_thick_OreZone", agg_method="max", interactive_map=True, interactive_max_sample=9000)
Make this Notebook Trusted to load map: File -> Trust Notebook
gs_pv_low = gs.to_pyvista('hard_top_Lower')
gs_pv_up = gs.to_pyvista('hard_top_Upper')
gs_pv_ore = gs.to_pyvista('hard_top_OreZone')
dh_pv = dh.to_pyvista('domain_code')
p = pv.Plotter()

p.add_mesh(gs_pv_low, color='r')
p.add_mesh(gs_pv_up, color='g', style='wireframe')
p.add_mesh(gs_pv_ore, color='b')

p.add_mesh(dh_pv.tube(radius=20))

p.set_scale(zscale=20)

p.show()
# export Top/Base OreZone surfaces as CSV and Drillholes
gs.to_csv("../data/hard_top_Lower.csv", ["X", "Y", "hard_top_Lower"], index=False)
gs.to_csv("../data/hard_top_OreZone.csv", ["X", "Y", "hard_top_OreZone"], index=False)
gs.to_csv("../data/hard_top_Upper.csv", ["X", "Y", "hard_top_Upper"], index=False)

gs.to_file("../data/surfaces")
dh.to_file("../data/domained_drillholes")