# framework
using GeoStats
# IO modules
using GeoIO
# viz modules
using PairPlots
import CairoMakie as Mke
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:
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:
= "https://zenodo.org/record/7051975/files/drillholes.csv?download=1"
url
= download(url, tempname()*".csv")
csv
= GeoIO.load(csv, coords = ("X", "Y", "Z"))
dtable
|> Select("Cu ppm") |> viewer dtable
┌ 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
|> describe dtable
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.
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:
- Preliminary analysis and processing
- Definition of interpolation domain
- Multivariate geostatistical interpolation
- 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:
= Select("HOLEID")
selectholeid
= Select("Cu ppm" => "Cu",
selectgrades "Au ppm" => "Au",
"Ag ppm" => "Ag",
"S ppm" => "S") →
Functional(x -> 1e-4*x*u"percent") # 1 ppm = 1e-4 percent
= selectholeid ⊔ selectgrades dclean
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)
= dclean(dtable) dtable
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
:
|> Select("Cu", "Au", "Ag", "S") |> values |> pairplot dtable
┌ 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
= boundingbox(dtable.geometry)
bbox
# size of blocks in meters
= (25.0u"m", 25.0u"m", 12.5u"m")
bsize
# define Cartesian grid
= CartesianGrid(extrema(bbox)..., bsize) grid
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)
current_figure() Mke.
┌ 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")
= shadow.(dtable.geometry)
points
= convexhull(points) chull
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")
current_figure() Mke.
┌ 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 Hexahedron
s for which the projected centroid is inside the convexhull
:
= findall(h -> shadow(centroid(h)) ∈ chull, grid)
active
= view(grid, active) blocks
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 Hexahedron
s 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
= @chain dtable begin
ztable @groupby(:HOLEID)
@transform(:z = zcoord(:geometry), :geometry = shadow(:geometry))
@combine(:z = first(:z), :geometry = first(:geometry))
end
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:
= unique(shadow.(centroid.(blocks)))
centroids
= ztable |> Select("z") |> Interpolate(centroids, IDW()) ztable
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) |
⋮ | ⋮ |
|> viewer ztable
┌ 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
Finally, we can filter the blocks for which the z
coordinate is below the terrain:
p(h) = shadow(centroid(h))
z(h) = zcoord(centroid(h))
= Dict(ztable.geometry .=> ztable.z)
zdict
= findall(h -> z(h) < zdict[p(h)], blocks)
active
= view(blocks, active) blocks
110854 view(::CartesianGrid, [76, 77, 78, 79, ..., 337107, 337109, 337193, 342719])
├─ 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: 924.995 m, y: -682.9 m, z: 337.58 m), ..., (x: 924.995 m, y: -657.9 m, z: 350.08 m))
├─ Hexahedron((x: 874.995 m, y: -757.9 m, z: 350.08 m), ..., (x: 874.995 m, y: -732.9 m, z: 362.58 m))
├─ Hexahedron((x: 924.995 m, y: -757.9 m, z: 350.08 m), ..., (x: 924.995 m, y: -732.9 m, z: 362.58 m))
├─ Hexahedron((x: 899.995 m, y: -732.9 m, z: 350.08 m), ..., (x: 899.995 m, y: -707.9 m, z: 362.58 m))
└─ Hexahedron((x: 924.995 m, y: -757.9 m, z: 362.58 m), ..., (x: 924.995 m, y: -732.9 m, z: 375.08 m))
viz(blocks)
┌ 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
The filtered blocks constitute our domain of interpolation.
12.3.3 Interpolation of grades
We saw that the distribution of chemical elements in the drill hole samples is very skewed. This is always the case in the mining industry. Another issue is that metal and mineral grades are examples of compositional data (Aitchison 1982). The values in these variables are constrained to live in the interval \([0,1]\) and to sum up to 100% if all chemical elements are considered.
12.3.3.1 Preprocessing
In order to remove compositional data constraints, we will perform the centered log-ratio transform (CLR
) from the CoDa.jl module:
= dtable |> Select("Cu", "Au", "Ag", "S")
grades
|> CLR() |> values |> pairplot grades
┌ 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
After the transform, the variables are free to vary in the unbounded interval \([-\infty,\infty]\). The theory behind this transform is beyond the scope of this book. Nevertheless, it is a simple mathematical expression in terms of logarithms of ratios (e.g., Cu/S).
Next, we attempt to transform the multivariate distribution to a multivariate standard normal using the ProjectionPursuit
transform:
|> CLR() |> ProjectionPursuit() |> values |> pairplot grades
┌ 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
The ProjectionPursuit
is an advanced statistical transform that removes non-linear associations between variables using an iterative procedure (Friedman 1987). The result is a set of independent variables that can be interpolated separately.
In order to “undo” these transforms after the interpolation, we create a revertible pipeline:
= CLR() → ProjectionPursuit()
preproc
= apply(preproc, grades)
samples, cache
samples
PP1 | PP2 | PP3 | PP4 | geometry |
---|---|---|---|---|
Continuous | Continuous | Continuous | Continuous | Point |
[NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | 🖈 Cartesian{NoDatum} |
-0.630435 | 1.53296 | -1.26015 | 1.29223 | (x: 559.725 m, y: -513.31 m, z: 252.82 m) |
-1.49701 | -0.289225 | -1.91429 | 1.17947 | (x: 558.955 m, y: -515.23 m, z: 246.87 m) |
-0.299542 | 0.91393 | 0.268712 | 1.28714 | (x: 557.225 m, y: -519.61 m, z: 233.37 m) |
-1.39207 | 1.53604 | 0.113834 | -0.45409 | (x: 555.375 m, y: -524.3 m, z: 218.92 m) |
-0.545962 | 1.24846 | -0.476919 | -0.423197 | (x: 553.825 m, y: -528.21 m, z: 206.94 m) |
-0.905358 | 2.00377 | -0.0110115 | 0.115843 | (x: 552.075 m, y: -532.56 m, z: 193.91 m) |
-0.294819 | 2.00301 | -1.87776 | 1.52664 | (x: 550.305 m, y: -537.08 m, z: 180.8 m) |
-1.91203 | 0.373393 | -0.318449 | 0.688997 | (x: 547.495 m, y: -544.31 m, z: 160.08 m) |
-0.0516281 | -1.31772 | -1.10327 | 1.36823 | (x: 544.335 m, y: -552.43 m, z: 137.08 m) |
-1.37098 | -0.547381 | -0.193644 | -0.157526 | (x: 541.695 m, y: -559.19 m, z: 117.96 m) |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
12.3.3.2 Geospatial correlation
Let’s fit a theoretical variogram for all four (independent) variables up to a given maximum lag:
= 300.0u"m"
maxlag
= setdiff(names(samples), ["geometry"])
ns
= [EmpiricalVariogram(samples, n, maxlag = maxlag) for n in ns]
gs
= [GeoStatsFunctions.fit(Variogram, g, h -> 1 / h^2) for g in gs] γs
4-element Vector{Variogram}:
MaternVariogram(sill: 0.879918, nugget: 0.440821, order: 1.0, range: 177.724 m, distance: Euclidean)
GaussianVariogram(sill: 0.940239, nugget: 0.644895, range: 97.522 m, distance: Euclidean)
MaternVariogram(sill: 0.900413, nugget: 0.547532, order: 1.0, range: 161.489 m, distance: Euclidean)
MaternVariogram(sill: 0.889081, nugget: 0.487782, order: 1.0, range: 241.576 m, distance: Euclidean)
We performed the fit
of the variogram model using the weighting function h -> exp(-h/100)
that penalizes the lag distance h
with an exponential model. The constant 100
was chosen based on visual inspection of the EmpiricalVariogram
estimates below.
function gammaplot(n, g, γ)
varioplot(g, axis = (; title = n))
varioplot!(γ, maxlag = maxlag, color = "teal")
current_figure()
Mke.end
gammaplot(ns[1], gs[1], γs[1])
┌ 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
gammaplot(ns[2], gs[2], γs[2])
┌ 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
gammaplot(ns[3], gs[3], γs[3])
┌ 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
gammaplot(ns[4], gs[4], γs[4])
┌ 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
Assuming that the variogram models are adequate, we can proceed to interpolation.
12.3.3.3 Geostatistical interpolation
Given the domain of interpolation, the samples and the variogram models, we can perform interpolation with InterpolateNeighbors
:
= [n => Kriging(γ) for (n, γ) in zip(ns, γs)]
models
= samples |> InterpolateNeighbors(blocks, models...) interp
PP1 | PP2 | PP3 | PP4 | geometry |
---|---|---|---|---|
Continuous | Continuous | Continuous | Continuous | Hexahedron |
[NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | 🖈 Cartesian{NoDatum} |
-0.283672 | -0.283745 | 0.263632 | 0.436704 | 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)) |
0.111823 | -0.197538 | 0.675881 | 0.76356 | 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)) |
0.111809 | -0.197538 | 0.675895 | 0.763337 | 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)) |
0.111799 | -0.197538 | 0.675903 | 0.763182 | 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)) |
0.0576 | -0.128827 | 0.145369 | -0.203214 | 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)) |
-0.344709 | -0.039885 | 0.416538 | 0.43668 | Hexahedron((x: 674.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 674.995 m, y: -832.9 m, z: -387.42 m)) |
-0.283676 | -0.283745 | 0.263636 | 0.436985 | Hexahedron((x: 699.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 699.995 m, y: -832.9 m, z: -387.42 m)) |
-0.283662 | -0.283745 | 0.263677 | 0.43681 | Hexahedron((x: 724.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 724.995 m, y: -832.9 m, z: -387.42 m)) |
0.111792 | -0.197538 | 0.675909 | 0.763355 | Hexahedron((x: 749.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 749.995 m, y: -832.9 m, z: -387.42 m)) |
0.111775 | -0.197538 | 0.675926 | 0.763087 | Hexahedron((x: 774.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 774.995 m, y: -832.9 m, z: -387.42 m)) |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
Let’s confirm that the interpolated values follow the same standard normal distribution:
|> Sample(10000) |> values |> pairplot interp
┌ 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
12.3.3.4 Postprocessing
In order to get the interpolated values in the original compositional space, we need to revert
the preprocessing pipeline:
= revert(preproc, interp, cache) estim
Cu | Au | Ag | S | geometry |
---|---|---|---|---|
Continuous | Continuous | Continuous | Continuous | Hexahedron |
[NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | 🖈 Cartesian{NoDatum} |
0.742081 | 2.97343e-5 | 0.000197271 | 0.257692 | 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)) |
0.744244 | 2.11916e-5 | 0.000205128 | 0.255529 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 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)) |
0.754422 | 3.40118e-5 | 0.000199597 | 0.245344 | 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)) |
0.742265 | 3.03434e-5 | 0.000208429 | 0.257496 | Hexahedron((x: 674.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 674.995 m, y: -832.9 m, z: -387.42 m)) |
0.742068 | 2.97418e-5 | 0.000197247 | 0.257705 | Hexahedron((x: 699.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 699.995 m, y: -832.9 m, z: -387.42 m)) |
0.74207 | 2.97406e-5 | 0.000197249 | 0.257703 | Hexahedron((x: 724.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 724.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.255529 | Hexahedron((x: 749.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 749.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11914e-5 | 0.000205128 | 0.255529 | Hexahedron((x: 774.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 774.995 m, y: -832.9 m, z: -387.42 m)) |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
|> Select("Cu") |> viewer estim
┌ 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
12.3.4 Model of recovery
We introduce a simplistic model of metallurgical recovery using the grade of copper estimated at the mining blocks. We assume that the logistic function represents an ideal behavior for the recovery as the grade of copper increases:
= mean(estim.Cu) - 0.1
μ = std(estim.Cu)
σ
f(Cu) = 1 / (1 + exp(-(Cu - μ) / σ))
f (generic function with 1 method)
= estim |> Map("Cu" => f => "f") estim
Cu | Au | Ag | S | f | geometry |
---|---|---|---|---|---|
Continuous | Continuous | Continuous | Continuous | Continuous | Hexahedron |
[NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | 🖈 Cartesian{NoDatum} |
0.742081 | 2.97343e-5 | 0.000197271 | 0.257692 | 0.97787 | 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)) |
0.744244 | 2.11916e-5 | 0.000205128 | 0.255529 | 0.979592 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 0.979591 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 0.979591 | 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)) |
0.754422 | 3.40118e-5 | 0.000199597 | 0.245344 | 0.98608 | 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)) |
0.742265 | 3.03434e-5 | 0.000208429 | 0.257496 | 0.978022 | Hexahedron((x: 674.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 674.995 m, y: -832.9 m, z: -387.42 m)) |
0.742068 | 2.97418e-5 | 0.000197247 | 0.257705 | 0.977859 | Hexahedron((x: 699.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 699.995 m, y: -832.9 m, z: -387.42 m)) |
0.74207 | 2.97406e-5 | 0.000197249 | 0.257703 | 0.977861 | Hexahedron((x: 724.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 724.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.255529 | 0.979592 | Hexahedron((x: 749.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 749.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11914e-5 | 0.000205128 | 0.255529 | 0.979592 | Hexahedron((x: 774.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 774.995 m, y: -832.9 m, z: -387.42 m)) |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
Please check the paper by Hoffimann et al. (2022b) for a more elaborate model of metallurgical recovery in the locked-cycle-test.
12.3.5 Economic assessment
Given the block model with the grade of copper and metallurgical recovery, we can proceed and apply the formula of economic value stated in our objectives:
= 1000u"kg"
ton
= 2.75 * ton / 1u"m^3"
ρ = 4000 / ton
P = 4 / ton
Cₘ = 10 / ton
Cₚ
= @transform(estim,
estim :value = volume(:geometry) * ρ * ((:Cu / 100) * :f * P - (Cₘ + Cₚ))
)
Cu | Au | Ag | S | f | value | geometry |
---|---|---|---|---|---|---|
Continuous | Continuous | Continuous | Continuous | Continuous | Continuous | Hexahedron |
[NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | [NoUnits] | 🖈 Cartesian{NoDatum} |
0.742081 | 2.97343e-5 | 0.000197271 | 0.257692 | 0.97787 | 3.22832e5 | 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)) |
0.744244 | 2.11916e-5 | 0.000205128 | 0.255529 | 0.979592 | 3.25751e5 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 0.979591 | 3.25751e5 | 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)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.25553 | 0.979591 | 3.2575e5 | 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)) |
0.754422 | 3.40118e-5 | 0.000199597 | 0.245344 | 0.98608 | 3.38525e5 | 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)) |
0.742265 | 3.03434e-5 | 0.000208429 | 0.257496 | 0.978022 | 3.23083e5 | Hexahedron((x: 674.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 674.995 m, y: -832.9 m, z: -387.42 m)) |
0.742068 | 2.97418e-5 | 0.000197247 | 0.257705 | 0.977859 | 3.22814e5 | Hexahedron((x: 699.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 699.995 m, y: -832.9 m, z: -387.42 m)) |
0.74207 | 2.97406e-5 | 0.000197249 | 0.257703 | 0.977861 | 3.22816e5 | Hexahedron((x: 724.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 724.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11915e-5 | 0.000205128 | 0.255529 | 0.979592 | 3.25751e5 | Hexahedron((x: 749.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 749.995 m, y: -832.9 m, z: -387.42 m)) |
0.744244 | 2.11914e-5 | 0.000205128 | 0.255529 | 0.979592 | 3.25751e5 | Hexahedron((x: 774.995 m, y: -857.9 m, z: -399.92 m), ..., (x: 774.995 m, y: -832.9 m, z: -387.42 m)) |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
We can then visualize all blocks with a positive economic value:
|> Filter(x -> x.value > 0) |> Select("value") |> viewer estim
┌ 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
Or any criterion of interest such as positive economic value and small fraction of contaminants:
|> Filter(x -> x.value > 0 && x.S < 0.25) |> Select("value") |> viewer estim
┌ 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
12.4 Summary
In this chapter, we illustrated an application of the framework in the mining industry. Among other things, we learned how to
- Perform simple economic assessment based on grades and metallurgical recoveries estimated at mining blocks using simple interpolation of transformed variables from drill hole samples.
- Use the tools covered in previous chapters to localize regions of interest in the mineral deposit.
Although the mathematical model presented here is simple, it is what most mining companies do. There is opportunity to improve these types of estimates with more sophisticated geospatial data science pipelines.