12  Mineral deposits

In the mining industry, resource estimation consists of interpolating measurements of metal and mineral grades from drill hole samples to 3D grids known as “block models”. Due to highly skewed distributions, several pre-processing steps need to be performed before the actual interpolation. In this chapter, we will cover simple steps for resource estimation and economic assessment of a real mineral deposit.

TOOLS COVERED: @groupby, @transform, @combine, CLR, ProjectionPursuit, EmpiricalVariogram, Kriging, Interpolate, InterpolateNeighbors, Shadow, Map, Filter, boundingbox, convexhull, viewer

MODULES:

# framework
using GeoStats

# IO modules
using GeoIO

# viz modules
using PairPlots
import CairoMakie as Mke
Note

Although we use CairoMakie.jl in this book, many of the 3D visualizations in this chapter demand a more performant Makie.jl backend. Consider using GLMakie.jl if you plan to reproduce the code locally.

12.1 Data

The GeoMet dataset (Hoffimann et al. 2022a) consists of three geospatial tables stored as CSV files. In this chapter, we will only use the drillholes.csv table.

Drill hole samples are always available in mining projects. They contain chemical information for each rock sample (a cylinder) along the drill hole trajectories. In this case, the data has been processed, and only the Cartesian “X”, “Y”, “Z” coordinates of the centroids of the cylinders were stored:

url = "https://zenodo.org/record/7051975/files/drillholes.csv?download=1"

csv = download(url, tempname()*".csv")

dtable = GeoIO.load(csv, coords = ("X", "Y", "Z"))

dtable |> Select("Cu ppm") |> viewer
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227
dtable |> describe
Table with 6 columns and 19 rows:
      variable  mean      minimum  median    maximum   nmissing
    ┌──────────────────────────────────────────────────────────
 1  │ HOLEID    nothing   1        nothing   119       0
 2  │ Ag ppm    1.86199   0.01     1.14      15.38     0
 3  │ Al ppm    52374.8   2400.0   56100.0   103300.0  0
 4  │ Au ppm    0.470965  0.0      0.21      5.89      0
 5  │ C ppm     1117.7    100.0    700.0     20700.0   0
 6  │ Ca ppm    9900.6    100.0    6400.0    73700.0   0
 7  │ Cl ppm    2303.75   5.48     1756.23   14273.8   0
 8  │ Cu ppm    7540.4    0.0      5000.0    60200.0   0
 9  │ F ppm     2528.58   20.89    1610.47   19786.3   0
 10 │ Fe ppm    245452.0  9900.0   246000.0  500000.0  0
 11 │ K ppm     13610.9   100.0    12450.0   45500.0   0
 12 │ Mg ppm    17673.0   200.0    14000.0   88800.0   0
 13 │ Mn ppm    4282.96   88.95    3685.35   24812.8   0
 14 │ Na ppm    4404.8    100.0    1400.0    58900.0   0
 15 │ P ppm     720.896   38.08    630.705   5576.14   0
 16 │ Pb ppm    9.3938    0.8      6.005     247.83    0
 17 │ S ppm     2712.95   100.0    1700.0    21500.0   0
 18 │ Th ppm    8.37926   0.37     4.08      94.81     0
 19 │ U ppm     19.0894   0.21     10.825    521.49    0

There are 18 chemical elements in the table, all measured in parts per million (ppm). The table also stores an integer identifier for each hole trajectory in the “HOLEID” column. There are 119 such trajectories as shown in the “maximum” column of the describe output.

Note

In most mining projects, the drill hole samples are available as “SURVEY”, “COLLAR” and “INTERVAL” tables, which can be desurveyed and composited with DrillHoles.jl.

12.2 Objectives

Our main objective is to estimate the economic value associated with each mining block in a 3D block model, i.e. a CartesianGrid with Hexahedron geometries (the blocks). This economic value in U$ dollars is estimated in terms of various other geospatial variables:

\[ Value = \underbrace{V \times \rho \times Cu \times f \times P}_{\text{revenue}} - \underbrace{V \times \rho \times (C_m + C_p)}_{\text{cost}} \]

where

  • \(V\) is the volume of the block in \(m^3\)
  • \(\rho\) is the rock density in \(ton/m^3\)
  • \(Cu\) is the grade of copper in \([0,1]\)
  • \(f\) is the recovery of copper in \([0,1]\)
  • \(P\) is the selling price in \(U\$/ton\)
  • \(C_m\) is the mining cost in \(U\$/ton\)
  • \(C_p\) is the plant cost in \(U\$/ton\)

Secondary objectives include the localization (through 3D visualization) of blocks with high economic value, high grades of Au and Ag, and low grade of S.

For simplicity, we assume the following constants:

  • \(\rho = 2.75\ ton / m^3\)
  • \(P = 4000\ U\$ / ton\)
  • \(C_m = 4\ U\$ / ton\)
  • \(C_p = 10\ U\$ / ton\)

12.3 Methodology

In order to estimate the economic value of each mining block, we need to interpolate the grade of Cu. Because we also want to localize the blocks with high grades of Au and Ag, and low grade of S, we will perform multivariate geostatistical interpolation of Cu, Au, Ag and S.

The proposed methodology has the following steps:

  1. Preliminary analysis and processing
  2. Definition of interpolation domain
  3. Multivariate geostatistical interpolation
  4. Economic assessment and visualizations

12.3.1 Preliminary analysis

We recommend to start any application discarding all information that is not relevant for the stated objectives. In this case, the geotable contains measurements of various chemical elements that are not used in the economic assessment. We define a cleaning pipeline that selects and renames columns of interest, and adds units to the measurements:

selectholeid = Select("HOLEID")

selectgrades = Select("Cu ppm" => "Cu",
                      "Au ppm" => "Au",
                      "Ag ppm" => "Ag",
                      "S ppm"  => "S") →
               Functional(x -> 1e-4*x*u"percent") # 1 ppm = 1e-4 percent

dclean = selectholeid ⊔ selectgrades
ParallelTableTransform
├─ Select(selector: [:HOLEID], newnames: nothing)
└─ SequentialTransform
   ├─ Select(selector: [Symbol("Cu ppm"), Symbol("Au ppm"), Symbol("Ag ppm"), Symbol("S ppm")], newnames: [:Cu, :Au, :Ag, :S])
   └─ Functional(selector: all, fun: #11)
dtable = dclean(dtable)
2000×6 GeoTable over 2000 PointSet
HOLEID Cu Au Ag S geometry
Categorical Continuous Continuous Continuous Continuous Point
[NoUnits] [%] [%] [%] [%] 🖈 Cartesian{NoDatum}
1 0.72 % 0.000104 % 6.4e-5 % 0.61 % (x: 559.725 m, y: -513.31 m, z: 252.82 m)
1 0.1 % 6.0e-6 % 2.0e-5 % 0.06 % (x: 558.955 m, y: -515.23 m, z: 246.87 m)
1 2.18 % 0.000202 % 0.000344 % 1.17 % (x: 557.225 m, y: -519.61 m, z: 233.37 m)
1 2.29 % 0.000343 % 0.000488 % 1.1 % (x: 555.375 m, y: -524.3 m, z: 218.92 m)
1 0.59 % 5.0e-5 % 0.000125 % 0.22 % (x: 553.825 m, y: -528.21 m, z: 206.94 m)
1 0.56 % 9.0e-5 % 9.0e-5 % 0.28 % (x: 552.075 m, y: -532.56 m, z: 193.91 m)
1 2.0 % 0.000356 % 9.3e-5 % 1.96 % (x: 550.305 m, y: -537.08 m, z: 180.8 m)
1 0.28 % 2.1e-5 % 7.6e-5 % 0.14 % (x: 547.495 m, y: -544.31 m, z: 160.08 m)
1 0.61 % 2.0e-5 % 0.000112 % 0.28 % (x: 544.335 m, y: -552.43 m, z: 137.08 m)
1 0.25 % 1.6e-5 % 7.7e-5 % 0.09 % (x: 541.695 m, y: -559.19 m, z: 117.96 m)
⋮ ⋮ ⋮ ⋮ ⋮ ⋮

In order to better understand the multivariate distribution of chemical elements, we visualize the values of the drill hole samples with the pairplot:

dtable |> Select("Cu", "Au", "Ag", "S") |> values |> pairplot
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

We can observe that the distribution is very skewed.

12.3.2 Domain of interpolation

Before we can interpolate these variables, we need to define our domain of interpolation. In this application, we will define a 3D CartesianGrid in terms of the drill hole trajectories alone. Some of the Hexahedron geometries will be disabled whenever they are outside the convexhull of the points.

First, let’s create our full CartesianGrid using the boundingbox of the trajectories:

# compute bounding box
bbox = boundingbox(dtable.geometry)

# size of blocks in meters
bsize = (25.0u"m", 25.0u"m", 12.5u"m")

# define Cartesian grid
grid = CartesianGrid(extrema(bbox)..., bsize)
85×66×66 CartesianGrid
├─ minimum: Point(x: -1150.004825000069 m, y: -882.9003199972212 m, z: -399.92 m)
├─ maximum: Point(x: 974.9951749999309 m, y: 767.0996800027788 m, z: 425.08 m)
└─ spacing: (25.0 m, 25.0 m, 12.5 m)
viz(dtable.geometry, color = "black")
viz!(grid, alpha = 0.2)
Mke.current_figure()
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

Second, let’s compute the convexhull of the Shadow of all points on the xy plane:

shadow(point) = point |> Shadow("xy")

points = shadow.(dtable.geometry)

chull = convexhull(points)
PolyArea
  outer
  └─ Ring((x: -1150.0 m, y: 641.23 m), ..., (x: -1120.56 m, y: 701.22 m))
viz(chull)
viz!(points, color = "black")
Mke.current_figure()
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227

We can filter the grid to retain Hexahedrons for which the projected centroid is inside the convexhull:

active = findall(h -> shadow(centroid(h)) ∈ chull, grid)

blocks = view(grid, active) 
136422 view(::CartesianGrid, [76, 77, 78, 79, ..., 370102, 370103, 370104, 370186])
├─ Hexahedron((x: 724.995 m, y: -882.9 m, z: -399.92 m), ..., (x: 724.995 m, y: -857.9 m, z: -387.42 m))
├─ Hexahedron((x: 749.995 m, y: -882.9 m, z: -399.92 m), ..., (x: 749.995 m, y: -857.9 m, z: -387.42 m))
├─ Hexahedron((x: 774.995 m, y: -882.9 m, z: -399.92 m), ..., (x: 774.995 m, y: -857.9 m, z: -387.42 m))
├─ Hexahedron((x: 799.995 m, y: -882.9 m, z: -399.92 m), ..., (x: 799.995 m, y: -857.9 m, z: -387.42 m))
├─ Hexahedron((x: 824.995 m, y: -882.9 m, z: -399.92 m), ..., (x: 824.995 m, y: -857.9 m, z: -387.42 m))
⋮
├─ Hexahedron((x: -900.005 m, y: 717.1 m, z: 412.58 m), ..., (x: -900.005 m, y: 742.1 m, z: 425.08 m))
├─ Hexahedron((x: -875.005 m, y: 717.1 m, z: 412.58 m), ..., (x: -875.005 m, y: 742.1 m, z: 425.08 m))
├─ Hexahedron((x: -850.005 m, y: 717.1 m, z: 412.58 m), ..., (x: -850.005 m, y: 742.1 m, z: 425.08 m))
├─ Hexahedron((x: -825.005 m, y: 717.1 m, z: 412.58 m), ..., (x: -825.005 m, y: 742.1 m, z: 425.08 m))
└─ Hexahedron((x: -900.005 m, y: 742.1 m, z: 412.58 m), ..., (x: -900.005 m, y: 767.1 m, z: 425.08 m))

We would also like to filter Hexahedrons that are above the terrain. Let’s create a simple terrain elevation model by interpolating the vertical z coordinate of the first point of each trajectory:

zcoord(point) = coords(point).z

ztable = @chain dtable begin
  @groupby(:HOLEID)
  @transform(:z = zcoord(:geometry), :geometry = shadow(:geometry))
  @combine(:z = first(:z), :geometry = first(:geometry))
end
119×3 GeoTable over 119 PointSet
HOLEID z geometry
Categorical Continuous Point
[NoUnits] [m] 🖈 Cartesian{NoDatum}
1 252.82 m (x: 559.725 m, y: -513.31 m)
2 280.94 m (x: 745.345 m, y: -280.6 m)
3 281.09 m (x: -304.025 m, y: 318.13 m)
4 256.16 m (x: -416.575 m, y: 570.35 m)
5 324.65 m (x: 198.815 m, y: -524.78 m)
6 317.1 m (x: 493.605 m, y: -530.83 m)
7 245.13 m (x: -367.965 m, y: 394.78 m)
8 250.67 m (x: 17.6752 m, y: 33.4297 m)
9 254.49 m (x: -514.305 m, y: 594.68 m)
10 294.44 m (x: -674.415 m, y: 596.84 m)
⋮ ⋮ ⋮

We perform the interpolation of the z coordinate on the projected centroids of the blocks:

centroids = unique(shadow.(centroid.(blocks)))

ztable = ztable |> Select("z") |> Interpolate(centroids, IDW())
2067×2 GeoTable over 2067 PointSet
z geometry
Continuous Point
[m] 🖈 Cartesian{NoDatum}
315.273 m (x: 737.495 m, y: -870.4 m)
318.783 m (x: 762.495 m, y: -870.4 m)
321.816 m (x: 787.495 m, y: -870.4 m)
324.655 m (x: 812.495 m, y: -870.4 m)
327.816 m (x: 837.495 m, y: -870.4 m)
309.922 m (x: 687.495 m, y: -845.4 m)
314.0 m (x: 712.495 m, y: -845.4 m)
318.755 m (x: 737.495 m, y: -845.4 m)
323.911 m (x: 762.495 m, y: -845.4 m)
327.361 m (x: 787.495 m, y: -845.4 m)
⋮ ⋮ 
ztable |> viewer
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/We6MY/src/scenes.jl:227