About me

Dr Kane Maxwell
Director & Principal Consultant
Matrix Geoscience

PhD, Geology & Machine Learning
MAusIMM (CP Geology)


https://matrixgeoscience.com/

kane.maxwell@matrixgeoscience.com

https://au.linkedin.com/in/maxwellk

  • Geologist and Data Scientist with more than 20 years of experience in the minerals and coal mining sector.

  • Director and Principal Consultant at Matrix Geoscience, applying advanced analytics and machine learning to geological and resource estimation problems.

  • Specialist in resource estimation, geostatistics, spatial modelling and machine learning.

Spatial data

What changes when observations have locations?

What is spatial data?

Spatial data are observations connected to a location or geometry.

In R, this may look like:

  • points: samples, drill holes, monitoring sites
  • lines: roads, rivers, faults, transects
  • polygons: boundaries, domains, catchments
  • rasters: grids, images, prediction surfaces

For this lecture, we mostly focus on spatial point/sample data.

A map is not just a picture of a table. It can reveal relationships that the table alone hides.

Why does location matter?

The same set of values can tell a very different story once we look at where they occur.

The question changes from:

what values do we have?

To:

how are those values arranged in space?

That is the starting point for spatial EDA.

Same distribution, different spatial structure: location changes what the data can tell us.

Where spatial data are used

Spatial methods are useful whenever observations are sparse, locations carry information, and decisions are made between samples.

Environmental and earth systems

  • groundwater contamination: where should monitoring occur?
  • soil chemistry and agriculture: where are nutrients deficient?
  • rainfall and climate surfaces: what happens between stations?
  • ecology: where are species likely to occur?
  • geology and mineral exploration: where should we drill next to improve resource confidence?

Human and engineered systems

  • disease surveillance: where are clusters emerging?
  • housing and urban data: how do neighbourhood effects matter?
  • infrastructure risk: where is exposure highest?
  • transport and access: which areas are poorly connected?

The shared structure is the same: location helps explain variation, and we often need to predict or decide between observations.

Spatial data in R

How do we store, map and analyse spatial objects?

From a table to a spatial object

A point dataset often starts as a normal table:

hole Easting Northing Ash
DH01 630120 7565320 10.4
DH02 630540 7565480 11.2

After conversion to an sf object, coordinates become a geometry column with a coordinate reference system.

Simple feature collection with 2 features
geometry type: POINT
CRS: AGD84 / AMG Zone 55

hole   Ash     geometry
DH01   10.4    POINT (630120 7565320)
DH02   11.2    POINT (630540 7565480)
R code 
holes_sf <- sf::st_as_sf(
  holes,
  coords = c("Easting", "Northing"),
  crs = 20355
)

The geometry column lets R understand location, not just store coordinates as ordinary numbers.

A practical spatial R toolkit

Common packages

  • sf: vector data, points, lines, polygons
  • ggplot2: mapping with geom_sf()
  • terra: raster and gridded data
  • gstat: variograms, interpolation, kriging

Useful resources

  • package vignettes and reference documentation
  • course examples and reproducible notebooks
  • RSpatial.org: practical spatial R tutorials and reference material
    https://rspatial.org/

The details change by problem, but the core workflow is usually similar.

A minimal spatial workflow in R

read data
  ↓
convert to spatial object
  ↓
map and explore
  ↓
model spatial structure
  ↓
predict / validate / decide

This is the workflow the lecture keeps returning to.

Compact code example 
holes <- readr::read_csv("data/holes.csv")

holes_sf <- sf::st_as_sf(
  holes,
  coords = c("Easting", "Northing"),
  crs = 20355
)

ggplot2::ggplot(holes_sf) +
  ggplot2::geom_sf(ggplot2::aes(colour = Ash))

vgm_emp <- gstat::variogram(Ash ~ 1, holes_sf)
vgm_fit <- gstat::fit.variogram(
  vgm_emp,
  gstat::vgm("Sph")
)

pred <- gstat::krige(
  Ash ~ 1,
  holes_sf,
  newdata = grid_sf,
  model = vgm_fit
)

Code is not the focus today, but the spatial workflow is reproducible and programmable.

Conventional EDA and spatial trend

Start with distributions, relationships, sampling pattern and trend.

Coordinates are not enough

A dataset can have coordinates without behaving like the kind of spatial process we want to model.

Coordinates attached

Households in a suburb:

  • address
  • number of people
  • income
  • dwelling type

Nearby houses are not necessarily more similar in household size than far houses.

Spatially structured variable

Coal quality, soil chemistry, grade, rainfall:

  • location matters
  • neighbourhood matters
  • continuity may exist
  • process context matters

Nearby samples may be more similar than distant samples.

For this lecture, spatial data means more than coordinates: it means there is evidence of spatial dependence worth modelling.

The working assumption

For the methods in this lecture, the central assumption is:

near samples tend to be more similar than far samples.

If this is supported

  • spatial EDA is meaningful
  • variograms are useful
  • interpolation may be defensible
  • spatial CV becomes important

If this is not supported

  • coordinates may not help
  • spatial prediction may be misleading
  • a different modelling strategy may be needed

Explore before modelling spatial dependence

Before modelling spatial dependence, ask:

Conventional EDA

  • distributions and skew
  • missingness and quality
  • unusual values
  • known groups or domains
  • relationships between variables

Spatial EDA

  • clustering and gaps
  • spatial trend
  • local anomalies
  • boundaries and changing processes
  • candidate groups/domains

The goal is to identify large-scale structure, sampling artefacts and possible groups/domains before choosing a spatial model.

Conventional EDA is still the starting point

What it does well

Conventional EDA helps us understand:

  • distributions and skew
  • unusual values
  • relationships between variables
  • missingness and data quality
  • differences between known groups/domains

The limitation is not EDA itself. The limitation is using EDA without asking where samples are, and what they are near.

A histogram or boxplot can flag a value as unusual; it cannot tell you whether it is isolated, supported by nearby samples, or part of a separate group or process.

EDA can show spatial trend

Trend is part of EDA. If the mean changes through space, we need to recognise that broad structure before modelling local spatial dependence.

A spatial trend can be two-dimensional

A simple trend surface is not the final spatial model.

It is a diagnostic question:

  • is there broad spatial drift?
  • do residuals still have structure?
  • are apparent anomalies partly trend effects?

This is why EDA comes before spatial modelling.

Spatial relationships can depend on direction

So far, we have mostly talked about separation as distance. But two observations the same distance apart may not have the same relationship.

  • Isotropic: dependence is mainly a function of distance
  • Anisotropic: dependence depends on distance and direction

Euclidean distance is one definition of separation, not a universal law. Direction, barriers, networks and time can all change what “near” means.

The variogram

How do we test whether near observations tend to be more similar than far observations?

Why construct a variogram?

At the start we made the key spatial assumption:

near samples tend to be more similar than far samples.

The variogram is one way to check and quantify that idea.

It asks:

  • does dissimilarity increase with distance?
  • how quickly does similarity decay?
  • over what distance does one sample help predict another?

The variogram cloud: every pair tells us something

Variogram clouds are useful for seeing scatter, short-distance discontinuities, possible outliers, boundaries or domain mixing before binning hides individual pairs.

Cloud to empirical variogram

Binning turns a noisy cloud of pairwise comparisons into an average dissimilarity trend by lag distance.

Why fit a variogram model?

The empirical variogram is useful, but it is still noisy and only exists at chosen lag distances.

A fitted model gives us a continuous, interpretable summary of spatial dependence over distance.

The fitted model is used to:

  • smooth the empirical relationship
  • estimate nugget, sill and range
  • provide a valid distance-to-dependence function
  • support model-based interpolation and uncertainty calculations

Anatomy of the fitted variogram

Semivariance: the calculation behind the plot

For pairs of samples separated by approximately the same lag distance \(h\):

\[ \gamma(h) = \frac{1}{2N(h)} \sum_{(i,j) \in h} \left[ Z(s_i) - Z(s_j) \right]^2 \]

  • \(Z(s_i)\): value at location \(i\)
  • \(Z(s_j)\): value at location \(j\)
  • \(h\): separation, or lag distance
  • \(N(h)\): number of pairs in that lag class
  • squared difference: pair dissimilarity
  • one half: geostatistical convention

Semivariance is half the average squared difference between observations separated by approximately the same distance.

What the terms mean for decisions

Nugget

Short-range variability, measurement/support effects, or structure below sampling scale.

Range

Approximate distance over which one sample informs another.

Sill

The variance level reached once samples are effectively no longer spatially correlated.

The range connects directly to sampling redundancy: two very close observations may add less independent information than observations separated by a meaningful fraction of the range.

Stationarity: the quiet assumption

Reasonable comparison

Pairs come from broadly the same process or population.

The mean and variance are reasonably stable across the area being modelled.

Potential problem

Pairs span strong trends, distinct environments, changing populations or boundaries between processes.

The variogram may still compute, but it may not answer the right question.

A variogram can be calculated even when the assumptions are poor. The question is whether the spatial relationship it summarises is meaningful.

Spatial EDA

What changes when exploration explicitly uses location?

Spatial EDA asks different questions

Spatial EDA is not only about finding extreme values. It asks how values relate to location, spatial regions and nearby samples.

Global spatial pattern

Is there a broad trend, gradient or clustering pattern across the whole study area?

Examples: value maps, smoothers, spatial clustering diagnostics.

Spatial domains

Do mapped regions, zones or environments behave differently?

Examples: compare values within spatially defined areas rather than pooling everything.

Local neighbourhoods

Is a sample surprising relative to nearby samples?

Example: nearest-neighbour residuals.

None of these automatically means bad data. They identify questions to investigate.

Global spatial EDA example: broad spatial pattern

Global spatial EDA asks whether the whole map shows a broad trend, gradient, cluster or gap pattern.

Spatial domain EDA example: compare mapped regions

Domain EDA becomes spatial when the groups are mapped regions, environments or zones. The same value may be ordinary in one domain and unusual in another.

Local spatial EDA example: nearest-neighbour residuals

A nearest-neighbour residual is one simple local spatial EDA method. It asks whether a value is unusual relative to the closest available spatial support, not just relative to the full dataset.

\[ r_i^{NN}=Z(s_i)-Z(s_{NN(i)}) \]

Tiny example

  • sample value: 15.8
  • nearest-neighbour value: 10.4
  • residual: 15.8 - 10.4 = +5.4

A large positive or negative value is locally surprising.

R code 
d <- sf::st_distance(samples)
diag(d) <- Inf

nn_id <- apply(d, 1, which.min)

samples$nn_residual <-
  samples$value - samples$value[nn_id]

A nearest-neighbour residual is a screen for investigation, not a deletion rule.

Spatial outlier screen: nearest-neighbour residual

A spatial outlier is not just an extreme value. It is a value that disagrees with its local support.

Spatial EDA also examines sampling geometry

Spatial EDA is not only about measured values.

It also asks how the samples are arranged:

Sampling density

  • where are samples clustered?
  • where are there gaps?
  • is coverage uniform?

Sampling direction

  • are samples aligned along traverses?
  • are some directions better informed?
  • are gaps process-related, logistical or accidental?

Sampling pattern affects what the data can tell us before we even look at the measured values.

What does the sampling pattern tell us?

A cluster may give detailed local information while adding little coverage elsewhere; a gap may dominate uncertainty.

Spatial prediction

Use spatial structure to estimate values between observations.

Why spatial prediction?

Spatial prediction algorithms try to infill between observations using the assumption:

near samples are more informative than far samples.

Common choices

  • IDW: distance-weighted average
  • Kriging: distance plus variogram model
  • Trend + residual spatial model: covariates first, residual dependence second

The method changes, but the spatial construct is the same: neighbourhood matters.

IDW: the simplest spatial interpolator

\[ \hat Z(s_0)=\frac{\sum_i w_i Z(s_i)}{\sum_i w_i},\qquad w_i=\frac{1}{d_i^p} \]

  • \(s_0\): unsampled prediction location
  • \(Z(s_i)\): observed value at nearby sample \(i\)
  • \(d_i\): distance from sample \(i\) to \(s_0\)
  • \(p\): power controlling distance decay
R idea 
weights <- 1 / distance^power
prediction <- sum(weights * value) / sum(weights)

IDW is easy to understand, but it does not estimate spatial covariance or model-based uncertainty.

IDW by hand: closer samples count more

The closer a sample is to the prediction point, the more influence it has. Increasing the power makes this local influence stronger.

Exact interpolation as a visual check

Why exact IDW can help EDA

  • exact at sample locations
  • highlights local influence
  • reveals bullseyes
  • exposes candidate domains

Why it is not the final answer

  • bullseyes are not automatically errors
  • no hard-boundary awareness
  • no model-based uncertainty

What does kriging add?

IDW

  • weights depend on distance
  • simple and exact
  • useful visual diagnostic
  • no variogram model
  • no kriging standard error

Kriging

  • weights depend on spatial covariance
  • uses a variogram model
  • can smooth noisy local effects
  • provides model-based prediction variance
  • still conditional on assumptions

Kriging adds a spatial covariance model, but it does not add contextual knowledge by itself.

Kriging is conditional on its assumptions

What kriging uses

  • sample values
  • sample locations
  • variogram model
  • local neighbourhood

What it does not know by itself

  • hard process boundaries
  • inappropriate domains or groups
  • invalid target areas
  • model misspecification

A smooth surface can hide a bad process assumption.

Prediction is not the same as certainty

Kriging standard error

Model-based uncertainty from:

  • data spacing
  • variogram
  • local support
  • edges and gaps

Contextual uncertainty

Interpretive uncertainty from:

  • hard boundaries
  • groups/domains
  • unmeasured drivers
  • target existence

Low kriging error does not guarantee low decision risk.

Spatial cross-validation

Validation should match the spatial prediction task.

Spatial cross-validation: random split is not the same task

Random 70/30 split

  • 70% train, 30% test
  • train and test samples are mixed through space
  • usually tests prediction among nearby samples
  • can be optimistic or pessimistic for the real task

Spatial cross-validation

  • test samples are withheld by area/block
  • the model predicts into less-supported space
  • closer to step-out or new-area prediction
  • matches the task more directly when location matters

A random split and a spatial split are not just different sampling choices. They test different prediction tasks.

What the splits look like in space

If the operational task is step-out prediction, the test set should look like step-out prediction.

Predictive performance depends on the validation design

Same synthetic spatial field, same IDW interpolation algorithm:

  • Random 70/30 split: repeated random train/test partitions
  • Spatial CV: leave-one-spatial-block-out folds
  • metric shown: MAE

The exact numbers are not the point. The validation design defines the task.

A random 70/30 score can be optimistic or pessimistic. The deeper problem is that it may answer the wrong question.

Three validation questions

The split should match the decision you care about.

1. Interpolation

Can I predict between nearby samples?

Random CV may be acceptable if this is truly the task.

2. Step-out

Can I predict into nearby but less-drilled space?

Use spatial blocks or leave-area-out tests.

3. Extrapolation

Can I predict into a new area or new domain?

This is the hardest problem and may not be supported by the data.

Spatial residual diagnostics

After fitting a model, is there spatial structure left behind?

Why inspect residuals?

After fitting a model, ask whether what remains is just random noise:

\[ Y_i=f(X_i)+\varepsilon_i \]

What the model explains

  • trend
  • covariates
  • groups
  • interactions

What we still need to check

  • residual maps
  • residual trend
  • residual variogram
  • domain or boundary effects

Spatially structured residuals do not automatically mean the model is useless. They mean the remaining errors may not be independent noise, and there may be spatial structure still left to explain.

Fit the trend first

A trend model can explain broad structure, but the residuals may still have spatial dependence.

Map the residuals

Residual maps ask whether model errors are spatially random or spatially organised.

Residual variogram

A residual variogram asks whether remaining errors are still spatially structured.

If residuals are spatially structured, what next?

A useful way to think about the workflow is:

\[ \text{observed pattern} = \text{explained trend} + \text{remaining spatial structure} \]

First improve the explainable part

  • check the trend model
  • add relevant covariates
  • reconsider groups or spatial domains
  • check whether boundaries matter

Then inspect what remains

  • map the residuals
  • check the residual variogram
  • ask whether nearby errors are still similar
  • model residual spatial structure if useful

The key question is simple: after the model, is the leftover pattern just noise, or does location still help explain it?

From lecture to worksheet

Use the workflow to make a spatial decision.

Case study: Mavis Downs

For the tutorial, we switch from generic examples to a synthetic mineral-exploration dataset based on the Mavis Downs deposit in Queensland.

  • existing drill-hole observations across a study area
  • uneven sampling coverage
  • mapped geological structure
  • 12 candidate drill sites

The detailed geological context is introduced in the tutorial, when it becomes decision-relevant.

Tutorial task

You will use the Mavis Downs case study to choose six new drill-hole locations.

Your decision must use:

  • sampling geometry
  • spatial continuity
  • variogram range
  • anomalies and domains
  • prediction surfaces
  • model uncertainty
  • geological boundaries
  • decision value

The highest-uncertainty location is not always the best drill target.

Take-home workflow

SPATIAL DATA → EDA & TREND → VARIOGRAM → SPATIAL EDA → PREDICT → VALIDATE & DIAGNOSE → DECIDE

  1. What is the spatial object, coordinate system and decision context?
  2. What do conventional EDA, trends, domains, clustering and gaps show?
  3. Does the variogram support the assumption that nearby samples are more similar?
  4. Do global, domain or local spatial EDA checks reveal anomalies or boundaries?
  5. Which prediction method is appropriate, and what uncertainty does it show?
  6. Does validation match the spatial prediction task, and do residuals still have structure?
  7. What action is justified by the data, uncertainty and context?