Kriging

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

pv.set_jupyter_backend("html")

geo.Project().set_crs(CRS("EPSG:26909"))

dh_tom = geo.read_file("../data/dh_tom_comp_capped.geo")

tmw_sol = geo.solid_from_dxf("TMW", "../data/tom_west_blob_small.dxf")

Grid Creation

To create the interpolation grid called BlockModel in GeoLime, we first create a Voxel wich is the plain parallelepiped. We then convert it into a BlockModel allowing us to remove the block outside the domain.

The Voxel is created from the drillholes bounds, they can be accessed using the minimum and maximum of each coordinates, or also using the bounds() methods.

z_max = np.max(dh_tom["Z_B"])
z_min = np.min(dh_tom["Z_E"])

z_min = dh_tom.bounds()[0][2]
z_max = dh_tom.bounds()[1][2]

z_min, z_max

(883.8832451628475, 1545.2686973526997)

As the bounds are floating values, it might be best to round them. Using the numpy round method, we can even round to the lower decaded using (-1), lower hundred using (-2), etc ..

z_min.round(-1)

880.0

For safety we also shift a bit the origin of the Voxel to ensure total coverage of the data.

Here we create a Voxel with a cell size of 15x15x15.

vx = geo.Voxel(
    name="Grid3D",
    shape=[50, 100, (z_max - z_min) // 5],
    origin=[
        dh_tom.bounds()[0][0].round(-1) - 5,
        dh_tom.bounds()[0][1].round(-1) - 5,
        z_min.round(-1) - 5
    ],
    axis=[
        [15, 0, 0],
        [0, 15, 0],
        [0, 0, 5],
    ]
)

bm = geo.BlockModel(name="Block3D", xyz=vx.coords())
bm.element_count()

660000

geo.plot_2d(
    [
        {"georef_object": vx, "property": "Z", "agg_method":"first"},
        {"georef_object": dh_tom, "property": "Zn_pct", "agg_method":'sum'},
    ],
    width=400, 
    height=650,
    colorscale='viridis'
)

Unnecessary blocks removal

To flag each cell of the BlockModel in the solid we can use the method

tmw_sol.contains(bm)

which will create a region in_TMW in the BM. The cell will be inside the solid if the centroid is in the sold.

We can also use the function compute_block_proportion_inside_solid which will compute the proportion of each cell in the solid based on the given discretisation. A discretisation of (1, 1, 1) is equivalent to the method explained above but the results will be a new property with 0 and 1 instead of a new region.

geo.compute_block_proportion_inside_solid(
    grid_object=bm, 
    solid=tmw_sol, 
    block_proportion_name='domain_prop', # name of the new property
    block_discretization=(2, 2, 1) # proportion will be calculated on centroid only
)

set(bm['domain_prop'])

{0.0, 0.25, 0.5, 0.75, 1.0}

We can then transform the property into a region using condition such as equality or threshold.

bm.set_region(
    name="Ore",
    data="domain_prop > 0"
)

Cells in a specific region can then be kept or removed.

bm.keep_only_cells('Ore')

Plotting for visual QC

geo.plot_2d(
    [
        {"georef_object": bm, "property": "Z", "agg_method":"first"},
        {"georef_object": dh_tom, "property": "Z_M", "agg_method":'first'},
    ],
    width=400, 
    height=650,
    colorscale='viridis'
)

bm.set_property(name="cell_size_z", data=f"{bm.axis()[2]}")

geo.plot_2d(
    [
        {"georef_object": bm, "property": "cell_size_z", "agg_method":"sum"},
    ],
    width=400, 
    height=650,
    colorscale='viridis'
)

dh_pv = dh_tom.to_pyvista("Zn_pct")
bm_pv = bm.to_pyvista()
tmw_sol_pv = tmw_sol.to_pyvista()

p = pv.Plotter()

p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(bm_pv,show_edges=True, opacity=0.7)
p.add_mesh(tmw_sol_pv, style='wireframe')
p.show()

2025-12-18 16:40:15.953 (  22.339s) [    7F7229FF4B80]vtkXOpenGLRenderWindow.:1458  WARN| bad X server connection. DISPLAY=

Covariance and Neighborhood setup

Covariance and neighborhood are necessary to initialise an Kriging estimator, here we defined a covariance based on the variography analysis and a appropriate Neighborhood strategy.

cov = (
    0.13 * geo.Nugget() 
    + 0.51 * geo.Spherical(metric=geo.Metric(scales=[60., 85., 10.], angles=[145.0, 45.0, 20.0]))
    + 0.34 * geo.Spherical(metric=geo.Metric(scales=[120., 90., 30.], angles=[145.0, 45.0, 20.0]))
)

neigh = geo.MaxPerCategoryNeighborhood(
    dimension=3,
    metric=geo.Metric(
        angles=[145, 45, 10],
        scales=[60, 60, 15]
    ),
    min_n=2,
    max_n=12,
    max_per_category=4,
    category=dh_tom.holeids()
)

Estimator Setup

GeoLime estimators are first initalised with the neccessary parameters. For the Kriging, covariance and neighborhod are needed.

estimator = geo.OrdinaryKriging(covariance_model=cov, neighborhood_model=neigh)

Now let’s imagine we have a list of grades sharing the same covariance model, we can use a loop to perform multiple kriging at once.

You can complexify and make nested loops if you have multiple domains and also multiple kriging pass.

for attrs in ['Zn_pct', 'Pb_pct']:
    estimator.solve(
        obj=dh_tom,
        obj_region=f"{attrs}.notna()",
        obj_attribute=attrs,
        support=bm,
        support_region=None,
        support_attribute=f"{attrs}_krig",
        discr=(3, 3, 3)
    )

Kriging Visual QC

geo.swath_plot(
    [
        {"obj":bm, "attribute":"Zn_pct_krig", "region":"Zn_pct_krig.notna()", "swath_interval":30, "axis":"Y", "color": "blue"},
        {"obj":dh_tom, "attribute":"Zn_pct", "region":"Zn_pct.notna()", "swath_interval":30, "axis":"Y_M", "color": "red"}
    ]
)

Using a rough estimation of the density block model we can then compute the grade tonnage curves.

bm.set_property(name='density', data='3')

grades = np.linspace(start=0, stop=15, num=10)

geo.gtc(grid=bm, attribute='Zn_pct_krig', density='density', grades=grades)

Tonnage and mean grades values can be extracted using the following functions.

tonnage = geo.gtc_plotting.compute_selectivity_tonnage(grid=bm, attribute='Zn_pct_krig', density='density', grades=grades)

tonnage

array([63655875., 43709625., 32538375., 19824750.,  9237375.,  3931875.,
        1397250.,   685125.,   371250.,   148500.])

mean_grade = geo.gtc_plotting.compute_selectivity_mean_grade(grid=bm, attribute='Zn_pct_krig', density='density', grades=grades)

They can be merged together in a DataFrame for a nicer display.

df = pd.DataFrame({'grades':grades, 'tonnage':tonnage, 'mean grade':mean_grade})

df

	grades	tonnage	mean grade
0	0.000000	63655875.0	3.657375
1	1.666667	43709625.0	5.081157
2	3.333333	32538375.0	5.959645
3	5.000000	19824750.0	7.076705
4	6.666667	9237375.0	8.567434
5	8.333333	3931875.0	10.139321
6	10.000000	1397250.0	12.231742
7	11.666667	685125.0	13.768232
8	13.333333	371250.0	15.072404
9	15.000000	148500.0	16.752916

For reporting, decimal values can be rounded.

df['grades'] = df['grades'].round(2)
df['tonnage'] = df['tonnage'].round(-3).astype(int)

df

	grades	tonnage	mean grade
0	0.00	63656000	3.657375
1	1.67	43710000	5.081157
2	3.33	32538000	5.959645
3	5.00	19825000	7.076705
4	6.67	9237000	8.567434
5	8.33	3932000	10.139321
6	10.00	1397000	12.231742
7	11.67	685000	13.768232
8	13.33	371000	15.072404
9	15.00	148000	16.752916

IDW for sanity check

Most of the time we also perform IDW for sanity check. In the same manner as the Kriging we first initialise the IDW estimator with the corresponding ellipsoid ans power and we solve it.

estimator = geo.IDW(
    metric=geo.Metric(
        angles=[145.0, 45.0, 20.0],
        scales=[60, 60, 15]
    ),
    power=2,
    neighbors_number=10
)

estimator.solve(
    obj=dh_tom,
    obj_region="Zn_pct.notna()",
    obj_attribute="Zn_pct",
    support=bm,
    support_region=None,
    support_attribute="Zn_idw",
)

metric=geo.Metric(
        angles=[145.0, -45.0, 20.0],
        scales=[3 * 60, 2* 60, 3 * 15]
    )

dh_pv = dh_tom.to_pyvista("Zn_pct")
bm_pv = bm.to_pyvista("Zn_pct_krig", "Zn_pct_krig.notna()")
ellipsoid = metric.to_pyvista(center=dh_tom.bounds()[1])

p = pv.Plotter()
p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(bm_pv,show_edges=True, opacity=1)
p.add_mesh(ellipsoid, style='wireframe')
p.remove_scalar_bar("Zn_pct_krig")
p.show()

slices = bm_pv.slice_orthogonal(x=442050, y=7003791, z=1380)
p = pv.Plotter()
p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(slices,show_edges=True, opacity=1)
p.remove_scalar_bar("Zn_pct_krig")
p.show()

bm.to_file("../data/block_model.geo")

Reporting

Results can then store in a nested dictionnary for resource classification

tonnage_zn = geo.gtc_plotting.compute_selectivity_tonnage(
    grid=bm, 
    attribute='Zn_pct_krig', 
    density='density', 
    grades=grades
)

tonnage_pb = geo.gtc_plotting.compute_selectivity_tonnage(
    grid=bm, 
    attribute='Pb_pct_krig', 
    density='density', 
    grades=grades
)

resource_dict = {}
resource_dict['TMZ'] = {}
resource_dict['TMZ']['Zn_krig'] = {}
resource_dict['TMZ']['Pb_krig'] = {}
for i, elt in enumerate(grades):
    resource_dict['TMZ']['Zn_krig'][elt] = tonnage_zn[i]
    resource_dict['TMZ']['Pb_krig'][elt] = tonnage_pb[i]    

Dictionnary can be converted in DataFrame and exported to excel file to keep fusionated cells.

final_classif = pd.DataFrame.from_dict(
    {(i,j): resource_dict[i][j] 
    for i in resource_dict.keys() 
    for j in resource_dict[i].keys()},
).T

final_classif

		0.000000	1.666667	3.333333	5.000000	6.666667	8.333333	10.000000	11.666667	13.333333	15.000000
TMZ	Zn_krig	63655875.0	43709625.0	32538375.0	19824750.0	9237375.0	3931875.0	1397250.0	685125.0	371250.0	148500.0
	Pb_krig	63595125.0	18096750.0	9018000.0	4556250.0	2585250.0	1471500.0	759375.0	337500.0	175500.0	81000.0

final_classif.to_excel('../data/gtc_report.xlsx', engine='xlsxwriter')