Estimation

Contents

Estimation#

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

pv.set_jupyter_backend('panel')


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

/tmp/ipykernel_3212/3172071959.py:6: PyVistaDeprecationWarning: `panel` backend is deprecated and is planned for future removal.
  pv.set_jupyter_backend('panel')

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

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

/opt/hostedtoolcache/Python/3.9.17/x64/lib/python3.9/site-packages/pyvista/core/grid.py:873: PyVistaDeprecationWarning: `UniformGrid` is deprecated. Use `ImageData` instead.
  warnings.warn(

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