12-mining.qmd

---
engine: julia
---

# Mineral deposits

```{julia}
#| echo: false
#| output: false
import Pkg
Pkg.activate(".")
```

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:**

```{julia}
# framework
using GeoStats

# IO modules
using GeoIO

# viz modules
using PairPlots
import CairoMakie as Mke
```

```{julia}
#| echo: false
#| output: false
Mke.activate!(type = "png")
```

::: {.callout-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.

:::

## Data

The [GeoMet](https://zenodo.org/record/7051975) dataset [@Hoffimann2022_1]
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:

```{julia}
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
```

```{julia}
dtable |> describe
```

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.

::: {.callout-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](https://github.com/JuliaEarth/DrillHoles.jl).

:::

## 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$

## 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

### 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:

```{julia}
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
```

```{julia}
dtable = dclean(dtable)
```

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

```{julia}
dtable |> Select("Cu", "Au", "Ag", "S") |> values |> pairplot
```

We can observe that the distribution is very skewed.

### 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:

```{julia}
# 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)
```

```{julia}
viz(dtable.geometry, color = "black")
viz!(grid, alpha = 0.2)
Mke.current_figure()
```

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

```{julia}
shadow(point) = point |> Shadow("xy")

points = shadow.(dtable.geometry)

chull = convexhull(points)
```

```{julia}
viz(chull)
viz!(points, color = "black")
Mke.current_figure()
```

We can filter the grid to retain `Hexahedron`s for which the projected centroid
is inside the `convexhull`:

```{julia}
active = findall(h -> shadow(centroid(h)) ∈ chull, grid)

blocks = view(grid, active) 
```

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:

```{julia}
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
```

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

```{julia}
centroids = unique(shadow.(centroid.(blocks)))

ztable = ztable |> Select("z") |> Interpolate(centroids, model=IDW())
```

```{julia}
ztable |> viewer
```

Finally, we can filter the blocks for which the `z` coordinate is below the terrain:

```{julia}
p(h) = shadow(centroid(h))
z(h) = zcoord(centroid(h))

zdict = Dict(ztable.geometry .=> ztable.z)

active = findall(h -> z(h) < zdict[p(h)], blocks)

blocks = view(blocks, active)
```

```{julia}
viz(blocks)
```

The filtered blocks constitute our domain of interpolation.

### 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** [@Aitchison1982]. 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.

#### Preprocessing

In order to remove compositional data constraints, we will perform the centered log-ratio
transform (`CLR`) from the [CoDa.jl](https://github.com/JuliaEarth/CoDa.jl) module:

```{julia}
grades = dtable |> Select("Cu", "Au", "Ag", "S")

grades |> CLR() |> values |> pairplot
```

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:

```{julia}
grades |> CLR() |> ProjectionPursuit() |> values |> pairplot
```

The `ProjectionPursuit` is an advanced statistical transform that removes non-linear associations
between variables using an iterative procedure [@Friedman1987]. 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:

```{julia}
preproc = CLR() → ProjectionPursuit()

samples, cache = apply(preproc, grades)

samples
```

#### Geospatial correlation

Let's fit a theoretical variogram for all four (independent) variables
up to a given maximum lag:

```{julia}
maxlag = 300.0u"m"

vs = setdiff(names(samples), ["geometry"])

gs = [EmpiricalVariogram(samples, v, maxlag = maxlag) for v in vs]

γs = [GeoStatsFunctions.fit(Variogram, g, h -> 1 / h^2) for g in gs]
```

::: {.callout-note}

We performed the `fit` of the variogram model using the weighting function
`h -> 1 / h^2` that penalizes the lag distance `h` with the inverse of its
square.

:::

```{julia}
function gammaplot(n, g, γ)
  fig = Mke.Figure()
  Mke.Axis(fig[1,1], title = n)
  funplot!(fig, g, maxlag = maxlag)
  funplot!(fig, γ, maxlag = maxlag)
end

gammaplot(vs[1], gs[1], γs[1])
```

```{julia}
gammaplot(vs[2], gs[2], γs[2])
```

```{julia}
gammaplot(vs[3], gs[3], γs[3])
```

```{julia}
gammaplot(vs[4], gs[4], γs[4])
```

Assuming that the variogram models are adequate, we can proceed to interpolation.

#### Geostatistical interpolation

Given the domain of interpolation, the samples and the variogram models, we
can perform interpolation with `InterpolateNeighbors`:

```{julia}
interps = map(vs, γs) do v, γ
  samples |> Select(v) |> InterpolateNeighbors(blocks, model=Kriging(γ))
end

interp = reduce(hcat, interps)
```

Let's confirm that the interpolated values follow the same standard normal distribution:

```{julia}
interp |> Sample(10000) |> values |> pairplot
```

#### Postprocessing

In order to get the interpolated values in the original compositional space, we need to
`revert` the preprocessing pipeline:

```{julia}
estim = revert(preproc, interp, cache)
```

```{julia}
estim |> Select("Cu") |> viewer
```

### 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:

```{julia}
μ = mean(estim.Cu) - 0.1
σ = std(estim.Cu)

f(Cu) = 1 / (1 + exp(-(Cu - μ) / σ))
```

```{julia}
estim = estim |> Map("Cu" => f => "f")
```

Please check the paper by @Hoffimann2022_2 for a more elaborate model of metallurgical
recovery in the locked-cycle-test.

### 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](#objectives):

```{julia}
ton = 1000u"kg"

ρ = 2.75 * ton / 1u"m^3"
P = 4000 / ton
Cₘ = 4 / ton
Cₚ = 10 / ton

estim = @transform(estim,
  :value = volume(:geometry) * ρ * ((:Cu / 100) * :f * P - (Cₘ + Cₚ))
)
```

We can then visualize all blocks with a positive economic value:

```{julia}
estim |> Filter(x -> x.value > 0) |> Select("value") |> viewer
```

Or any criterion of interest such as positive economic value and
small fraction of contaminants:

```{julia}
estim |> Filter(x -> x.value > 0 && x.S < 0.25) |> Select("value") |> viewer
```

## 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.