Estimation

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

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

Read Surface and Drillholes files

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

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

gs = geo.read_file("../data/surfaces.geo")
nu, nv = gs.shape
size_x, size_y = gs.axis()
size_z = 1
gs.properties()

['X',
 'Y',
 'Z',
 'U',
 'V',
 'hard_top_Upper',
 'hard_thick_Upper',
 'hard_top_OreZone',
 'hard_thick_OreZone',
 'hard_top_Lower']

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

(478.6, 389.2)

Create Estimation Grid

vx = geo.Voxel(
    name="Grid3D",
    shape=[nu - 2, nv - 2, (z_max - z_min) // size_z],
    origin=[
        gs.origin[0].round() + size_x,
        gs.origin[1].round() + size_y,
        z_min.round()
    ],
    axis=[
        [size_x, 0, 0],
        [0, size_y, 0],
        [0, 0, size_z],
    ]
)

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

761484

Flag grid cells position compared to surfaces

geo.compute_surface_elevation_intersection(
    obj=bm,
    grid_object=gs,
    z_prop="hard_top_OreZone",
    new_property_name="Z_top_OreZone"
)

geo.compute_surface_elevation_intersection(
    obj=bm,
    grid_object=gs,
    z_prop="hard_top_Lower",
    new_property_name="Z_top_Lower"
)

geo.compute_surface_elevation_intersection(
    obj=bm,
    grid_object=gs,
    z_prop="hard_top_Upper",
    new_property_name="Z_top_Upper"
)

geo.compute_block_proportion_inside_solid(
    grid_object=bm, 
    solid=domain_solid, 
    block_proportion_name='domain_prop', # name of the new property
    block_discretization=(3,3,3) # proportion will be calculated on 81ths
)

bm.set_property_expr(name="domain_code", expr="1")
bm.set_property_expr(name="domain_code", expr="3", region="(Z > Z_top_Lower) & (Z < Z_top_OreZone)")
bm.set_property_expr(name="domain_code", expr="2", region="(Z >= Z_top_OreZone)")

bm.set_region_condition(
    name="OreZone_Domain",
    condition="(Z > Z_top_Lower) & (Z < Z_top_OreZone)"
)

bm.set_region_condition(
    name="OVB",
    condition="(Z >= Z_top_OreZone)"
)

bm.set_region_condition(
    name="BDR",
    condition="(Z <= Z_top_Lower) and (domain_prop < 0.1)"
)

bm.set_region_condition(
    name="Ore",
    condition="domain_prop >= 0.1"
)

dh_pv = dh.to_pyvista('Fe_pct')
bm_pv = bm.to_pyvista('domain_code')
domain_solid_pv = domain_solid.to_pyvista()

p = pv.Plotter()

p.add_mesh(dh_pv.tube(radius=10))
p.add_mesh(bm_pv, opacity=0.3)
p.add_mesh(domain_solid_pv, style='wireframe')

p.set_scale(zscale=10)

p.remove_scalar_bar(title='domain_code')

p.show()

Estimation

domain_solid.contains(dh)

dh.set_region_condition("HighGrade", "in_OreZone and (Fe_pct == Fe_pct)")

cov = (
    16.2 * geo.Nugget() 
    + 53.2 * geo.Spherical(metric=geo.Metric(scales=[150., 60., 4.]))
    + 22.8 * geo.Spherical(metric=geo.Metric(scales=[210, 110, 14]))
)

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

estimator = geo.OrdinaryKriging(cov, neigh)

estimator.solve(
    obj=dh,
    obj_region="HighGrade",
    obj_attribute="Fe_pct",
    support=bm,
    support_region="Ore",
    support_attribute="Fe_kriged",
    kriging_slope=False,
    kriging_efficiency=False
)

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

dh.set_region_condition('valid_holes', "HighGrade")

dh.keep_composites('valid_holes')

bm.keep_only_cells('Ore')

dh_pv = dh.to_pyvista('Fe_pct')
bm_pv = bm.to_pyvista('Fe_kriged')

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