Kriging

Kriging#

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

pv.set_jupyter_backend('html')
geo.Project().set_crs(CRS("EPSG:20350"))

Read Surface and Drillholes files#

dh = geo.read_file("../data/domained_drillholes.geo")

bm_panel = geo.read_file("../data/block_model_panel.geo")
bm_smu = geo.read_file("../data/block_model_smu.geo")

domain_solid = geo.datasets.load("rocklea_dome/domain_mesh.dxf")
domain_solid.name = "channel_hg"

Estimation#

domain_solid.contains(dh)

dh.set_region("channel_hg", "in_channel_hg & Fe_pct.notna()")

dh.regions()

['in_channel_hg', 'channel_hg']

Read Covariance Model from json file

with open("../data/cov_model.json") as json_data:
    cov_dict = json.load(json_data)

cov = cov_dict["nugget"] * geo.Nugget()
for i in range(len(cov_dict["structures"])):
    cov += cov_dict["structures"][i]["contribution"] * getattr(
        geo, 
        cov_dict["structures"][i]["structure_type"]
    )(
        metric=geo.Metric(
            scales=list(cov_dict["structures"][0]["anisotropy"]["ranges"].values()),
            angles=list(cov_dict["structures"][0]["anisotropy"]["rotation"].values()))
    )

For this estimation a MinMaxPointsNeighborhood is used. It allows to specify an ellipsoid of research and also a required minimum and allowed maximum of neighbors in the ellipsoid. If the minimum is not reached, the block is not estimated and if the maximum is exceeded, only the N first will be used in the computation.

neigh = geo.MinMaxPointsNeighborhood(
    dimension=3,
    metric=geo.Metric(
        angles=[0, 0, 0],
        scales=[3000, 3000, 20]
    ),
    min_n=1,
    max_n=20,
)

estimator = geo.OrdinaryKriging(cov, neigh)

estimator.solve(
    obj=dh,
    obj_region="channel_hg",
    obj_attribute="Fe_pct",
    support=bm_panel,
    support_region="channel_hg",
    support_attribute="Fe_kriged",
    kriging_slope=False,
    kriging_efficiency=False,
    discr=np.array([2, 2, 1])

)

estimator.solve(
    obj=dh,
    obj_region="channel_hg",
    obj_attribute="Fe_pct",
    support=bm_smu,
    support_region="channel_hg",
    support_attribute="Fe_kriged",
    kriging_slope=False,
    kriging_efficiency=False,
    discr=np.array([2, 2, 1])

)

bm_panel.to_file("../data/block_model_panel")
bm_smu.to_file("../data/block_model_smu")

The results can be visualised in 2D on a map. As the 3D data will be displayed on a 2D map, an aggregation operation is needed.

For the drillholes we will display the maximum Iron percentage per drillholes if there is at least one composite per hole in the “channel_hg” region. Following the same idea, for the BlockModel we will display the mean Iron percentage per column of block (being the blocks which share the same X and Y coordinates but have a different elevation).

dh.plot_2d(property="Fe_pct", agg_method="max", interactive_map=True, region='channel_hg', region_agg='any')
Make this Notebook Trusted to load map: File -> Trust Notebook

And in the same manner we will display each column which have at least one block in the “channel_hg” region. region_agg

bm_smu.plot_2d(property="Fe_kriged", agg_method="mean", interactive_map=True, region="channel_hg", region_agg="any")
Make this Notebook Trusted to load map: File -> Trust Notebook
dh_pv = dh.to_pyvista('Fe_pct', region='channel_hg')
bm_pv = bm_smu.to_pyvista('Fe_kriged', region='channel_hg')

p = pv.Plotter()

p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(bm_pv,show_edges=True, opacity=0.7)
p.remove_scalar_bar(title='Fe_pct')
p.set_scale(zscale=20)
p.show()
Contents