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)