12. Prediction, blocks and surfaces

A model’s output is not a map, it is a container full of columns. This chapter is the road from those columns to the objects a mining workflow actually exchanges: a refined block model, a grade shell or domain boundary as a triangulated surface, volumes and tonnages cut against geometry, and files other software opens.

The running example rebuilds chapter 8’s ore pod, refined.

import os
import tempfile

import geoml
import numpy as np

geoml.set_seed(1234)
rng = np.random.default_rng(1234)

xyz = rng.uniform(0, 160, size=[500, 3])
radius = np.linalg.norm(xyz - np.array([80.0, 80.0, 80.0]), axis=1)

point = geoml.data.PointData.from_array(xyz)
point.add_continuous_variable(
    "au", 4.0 * np.exp(-(radius / 28.0) ** 2) + 0.02)
point.get("au").set_cutoffs([1.0])

inducing = geoml.data.inducing.from_kmeans(point, 150, seed=0)

root = geoml.latent.BasicInput(
    inducing,
    transform=geoml.transform.Isotropic(22.0))

gp = geoml.latent.BasicGP(
    root,
    size=1,
    kernel=geoml.kernels.Gaussian())

warping = geoml.warping.ChainedWarping(
    geoml.warping.Softplus(1),
    geoml.warping.ZScore(1))

model = geoml.models.VGPNetwork(
    point, "au",
    geoml.likelihood.Gaussian(warping),
    gp,
    options=geoml.models.GPOptions(verbose=False))

model.train_full(max_iter=150)

blocks = geoml.data.BlockSet3D(
    start=[10, 10, 5],
    n=[8, 8, 16],
    step=[20.0, 20.0, 10.0],
    discretization=(2, 2, 2),
    max_levels=2)

blocks = geoml.models.refine(model, blocks, n_sim=20)
print(blocks.n_data, "blocks after refinement")

12.1 From a field to a surface

Any contoured value is a surface: a grade shell at a cut-off, a domain boundary at a categorical’s log-odds zero (chapter 6), a topography at an elevation. On a block model the extraction is careful in a way worth knowing about. A coarse block beside finer neighbours would tear the surface along their shared face, so the mesh handed to the contouring is locally cut down to the finest size, in the export only and never in the model. The result is smoother and closer to the true level set than a uniformly fine model several times the size.

shell = blocks.get_contour("au", 1.0)

print(type(shell).__name__,
      "closed:", shell.closed,
      "volume:", round(float(shell.volume), 0))

What comes back is a real mesh object with invariants rather than a soup of triangles. Surface3D never closes, Solid3D always closes and knows its volume, and a mesh that satisfies neither is refused rather than passed along quietly. Solids support union, intersection and difference, with a robust fallback where exact geometry fails, which warns and states its resolution when it is used. simplify(max_error) decimates against a geometric budget it actually verifies, and from_dxf and export_dxf exchange geometry with the rest of the mining world.

12.2 Geometry cutting data

The opposite direction, geometry deciding about locations, is one call on any container: assign_from_surface for above or below a sheet, and assign_from_solid for inside or outside a body, each writing a metadata column. Blocks add fraction=, the share of each block inside, measured on the sub-blocks the model already discretizes into.

blocks.assign_from_solid(shell, "shell", fraction="share")

share = blocks.get_metadata("share")
inside = float((blocks.block_volume * share).sum())

print("volume inside the shell:", round(inside, 0),
      "of", round(float(shell.volume), 0))

The two volumes agree to a couple of percent, which is the resolution of the test rather than a disagreement about the shape: the block-side number asks eight sub-blocks per block whether they are inside, while the mesh knows its own volume exactly. Refine further and the two converge. The second number is the one to quote when the geometry is the reference, and the point of computing both is that they check each other. The same calls put a topography or a lease boundary into the workflow: assign, then hand the mask to refine(..., where=...) so that ground nobody asked about is never predicted at all.

12.3 Leaving the package

Three doors, by destination:

frame = blocks.as_data_frame()
print(list(frame.columns)[:6])

mesh = blocks.as_pyvista()
print(mesh.n_cells, "hexahedra for a 3D viewer")

out = os.path.join(tempfile.mkdtemp(), "shell.dxf")
shell.export_dxf(out)
print("dxf written:", os.path.getsize(out) > 0)

as_data_frame gives one row per block with a size per row, which reaches spreadsheets and general software. as_pyvista reaches 3D viewers and VTK-based pipelines, with every column travelling as cell data and its path carried alongside. DXF reaches CAD and the mine-planning packages, geometry only. And to_zarr (chapter 10) remains the lossless hand-off between geoML scripts: the others are exports, this one is the container itself.

In the code. BlockSet3D.get_contour and the mesh hierarchy in data/meshes.py, with topo.clip_meshes for batch topography cuts; assignments in data/containers.py and data/blocks.py; the pure mesh mathematics in math/geometry.py. Design record: docs/variable-block-models.md.

12. Prediction, blocks and surfaces

12.1 From a field to a surface

12.2 Geometry cutting data

12.3 Leaving the package

Further reading