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")
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[3], line 1
----> 1 dh_tom = geo.read_file("../data/dh_tom_comp_capped.geo")
File ~/work/fireweed-report/fireweed-report/.venv/lib/python3.12/site-packages/geolime/objects/object_io.py:239, in read_file(filepath, attributes)
217 def read_file(
218 filepath: str,
219 attributes: list[str] | None = None,
220 ):
221 """
222 Read the object from its folder path. If it does not end with `.geo`,
223 the suffix will automatically be appended.
(...) 237
238 """
--> 239 manifest = ObjectIO.read_manifest(filepath)
240 object_type = get(manifest, ["meta", "object"])
242 module = sys.modules[__package__]
File ~/work/fireweed-report/fireweed-report/.venv/lib/python3.12/site-packages/geolime/objects/object_io.py:77, in ObjectIO.read_manifest(filepath)
75 manifest_path = os.path.join(filepath, "MANIFEST.toml")
76 if not os.path.exists(manifest_path):
---> 77 raise ValueError("Invalid GeoLime object path")
79 with open(manifest_path, "rb") as f:
80 return tomllib.load(f)
ValueError: Invalid GeoLime object path
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
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)
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()
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'])
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",
"show_deviated_traces": False,
},
],
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", color="black", line_width=2)
p.show()
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]))
)
Estimator Setup#
GeoLime estimators are first initalised with the neccessary parameters. For the Kriging, covariance and neighborhod are needed.
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"]:
region = f"{attrs}.notna()"
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.data("HOLEID", region),
)
estimator = geo.OrdinaryKriging(covariance_model=cov, neighborhood_model=neigh)
estimator.solve(
obj=dh_tom,
obj_region=region,
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
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
For reporting, decimal values can be rounded.
df['grades'] = df['grades'].round(2)
df['tonnage'] = df['tonnage'].round(-3).astype(int)
df
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 a neighborhood (search ellipsoid and number of neighbors) and a power, and we solve it.
estimator = geo.IDW(
neighborhood_model=geo.MinMaxPointsNeighborhood(
dimension=3,
metric=geo.Metric(
angles=[145.0, 45.0, 20.0],
scales=[60, 60, 15]
),
min_n=1,
max_n=10,
),
power=2,
)
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
final_classif.to_excel('../data/gtc_report.xlsx', engine='xlsxwriter')