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