Localised Uniform Conditioning

Localised Uniform Conditioning#

Library Import#

import numpy as np
import geolime as geo
import pyvista as pv
from pyproj import CRS
import matplotlib.pyplot as plt

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

Read Surface and Drillholes files#

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

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

gs = geo.read_file("../data/surfaces.geo")
nu, nv = np.array(gs.shape) // 2
size_x, size_y = np.array(gs.axis()) * 2
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 Panel 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_panel = geo.BlockModel(name="Block3D", xyz=vx.coords())
bm_panel.element_count()

181560

Flag panel cells position compared to surfaces#

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

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

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

geo.compute_block_proportion_inside_solid(
    grid_object=bm_panel, 
    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_panel.set_property_expr(name="domain_code", expr="1")
bm_panel.set_property_expr(name="domain_code", expr="3", region="(Z > Z_top_Lower) & (Z < Z_top_OreZone)")
bm_panel.set_property_expr(name="domain_code", expr="2", region="(Z >= Z_top_OreZone)")

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

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

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

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

# bm_panel.delete_cells('OVB')

Create SMU Estimation Grid#

smu_size = [size_x / 2, size_y / 2, 1]

smu_size

[50.0, 50.0, 1]

smu_coords = np.vstack(
    (
        np.hstack((bm_panel.coords() + [-25, -25, 0], bm_panel[['domain_code', 'Ore']])),
        np.hstack((bm_panel.coords() + [-25, 25, 0], bm_panel[['domain_code', 'Ore']])),
        np.hstack((bm_panel.coords() + [25, 25, 0], bm_panel[['domain_code', 'Ore']])),
        np.hstack((bm_panel.coords() + [25, -25, 0], bm_panel[['domain_code', 'Ore']])),
    )
)

bm_smu = geo.BlockModel('', smu_coords[:,:3].astype(float))
bm_smu.set_property('domain_code', smu_coords[:,3])
bm_smu.set_region('Ore', smu_coords[:,4])

Estimation#

domain_solid.contains(dh)

dh.set_region_condition("Ore", "in_OreZone and Fe_pct.notna()")

Change Of Support#

anam = geo.Anamorphosis(
    data=dh.data("Fe_pct", "Fe_pct.notna()"), 
    coefficient_number=50
)


x = np.linspace(-5, 5, 100)
fig, ax = plt.subplots()
ax.plot(x, anam.theoretical_transform(x), color ='b')
ax.set_title("Hermite expansion evaluation")
plt.show()
../_images/a7fee2a7efdc67783e92b4ebf767cc5ab394dfdb2f9b0c84e520130a2d2a6c91.png

Covariance & Neighborhood settings#

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

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

Localised Uniform Conditioning#

uc_estimator = geo.LocalisedUniformConditioning(
    covariance_model=covariance, 
    neighborhood_model=neighborhood,
)

ore = uc_estimator.solve(
    obj=dh, 
    obj_region='Ore',
    obj_attribute='Fe_pct',
    anamorphosis=anam, 
    support=bm_smu, 
    support_region='Ore', 
    panel_support=bm_panel, 
    panel_support_region='Ore', 
    cutoff_grades=np.arange(0., 60., 5, dtype=np.float64), 
    support_attribute='Fe_LUC', 
    discr=np.array([2, 2, 1],dtype=np.int64),
    panel_discr=np.array([2, 2, 1],dtype=np.int64)
)

Results Validation#

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

bm_smu_pv = bm_smu.to_pyvista(properties="Fe_LUC")

p = pv.Plotter()
p.add_mesh(bm_smu_pv)
p.set_scale(zscale=20)
p.show()

QA-QC#

dh.keep_composites("Ore")

geo.swath_plot(
    [
        {"obj":bm_smu, "attribute":"Fe_LUC", "region":None, "swath_interval":50, "axis":"Y", "color": "blue"},
        {"obj":dh, "attribute":"Fe_pct", "region":None, "swath_interval":50, "axis":"Y_M", "color": "red"}
    ]
)
Contents