Data Break 8: Kriging the Meuse River

BIOS 737 Spring 2004

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.1/27
 Meuse River: Reminder

• library(gstat)
• Data included in gstat library.
• data(meuse) makes the data available.
• Soil samples from flood plain of Meuse River, near
village of Stein (Belgium).
• Geostatistical data: locations s and attributes Z(s).

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.2/27
 Kriging steps

• Empirical semivariogram.
• Fit theoretical semivariogram to empirical
semivariogram values.
• Krige with best fitting semivariogram values.

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.3/27
 What do we have? (Zinc)
zinc concentrations (ppm)

333000

332000
100
200
400
y

800
1600

331000

330000

178500 179000 179500 180000 180500 181000 181500

x

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.4/27
 What do we want? (Predictions here).

Prediction Grid

333000
332000
y

331000
330000

178500 179000 179500 180000 180500 181000 181500

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.5/27
 Step 1: Data Gaussian?

• Recall from Data Break 6, log(zinc) closer to

Gaussian than zinc.
• We will krige the log(zinc) values.
• (Note: biased prediction of exp(log(zinc)) = zinc,
but we ignore this in the example).

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.6/27
 Step 2: Empirical Semivariogram(s)

logzinc.var <- variogram(log(zinc)˜1,

loc=˜x+y,
data=meuse)
plot(logzinc.var$dist,
logzinc.var$gamma,
xlab="distance",
ylab="semivariance")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.7/27
 Resulting plot

0.6

0.4
semivariance

0.2

0.0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.8/27
 Some other cool plots

• The variogram cloud consists of the plot of all

contrasts Z(si ) − Z(sj ) plotted versus ||si − sj ||

plot(variogram(log(zinc)˜1,
loc=˜x+y,
data=meuse, cloud=TRUE))
• Also some need interactive commands here....

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.9/27
 Variogram cloud

3
semivariance

0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.10/27
 Step 2: Fit Theoretical
Semivariograms

• fit.variogram function in gstat

• Not entirely clear what optimization rule, based on
practical experience rather than theoretical
motivation.
• To get best fitting spherical semivariogram...

model.1 <- fit.variogram(logzinc.var,

vgm(psill=1,model="Sph",
range=300,nugget=1))
plot(logzinc.var, model=model.1)
title("Best Fit Spherical")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.11/27
 Best spherical

Best Fit Spherical

0.6

0.4
semivariance

0.2

0.0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.12/27
 Step 2: Fit Theoretical
Semivariograms

• To get best fitting exponential semivariogram...

model.1 <- fit.variogram(logzinc.var,

vgm(psill=1,model="Exp",
range=300,nugget=1))
plot(logzinc.var, model=model.1)
title("Best Fit Exponential")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.13/27
 Best Exponential

Best Fit Exponential

0.6

0.4
semivariance

0.2

0.0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.14/27
 Step 2: Fit Theoretical
Semivariograms

• To get best fitting Gaussian semivariogram...

model.1 <- fit.variogram(logzinc.var,

vgm(psill=1,model="Gau",
range=300,nugget=1))
plot(logzinc.var, model=model.1)
title("Best Fit Gaussian")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.15/27
 Best Gaussian
Best Fit Gaussian

0.6

0.4
semivariance

0.2

0.0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.16/27
 Step 2: Fit Theoretical
Semivariograms

• To get best fitting Matern semivariogram...

model.1 <- fit.variogram(logzinc.var,

vgm(psill=1,model="Mat",
range=300,nugget=1))
plot(logzinc.var, model=model.1)
title("Best Fit Matern (kappa=0.25)")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.17/27
 Best Matern (κ = 0.25)

Best Fit Matern (kappa=0.5)

0.6

0.4
semivariance

0.2

0.0
0 500 1000 1500
distance

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.18/27
 What’s best?

• We’ll use the spherical theoretical semivariogram.

• (Skipping formal AIC or other comparison).
• What are the best fit values?

> model.1
model psill range
1 Nug 0.0506555 0.0000
2 Sph 0.5906009 896.9702
• Nugget: 0.051, Partial sill: 0.591, Range: 896.97

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.19/27
 Kriging

# define the theoretical semivariogram

m <- vgm(psill=0.591,"Sph",
range=896.97,nugget=0.051)
# ordinary kriging:
logzinc.krige <- krige(log(zinc)˜1,
˜x+y, model = m, data = meuse,
newd = meuse.grid)
# make ’levelplot’
levelplot(var1.pred˜x+y,
data=logzinc.krige,
aspect = mapasp(logzinc.krige),
main = "ordinary kriging predictions")

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.20/27
 Prediction levels
ordinary kriging predictions

−7.5

333000
−7

−6.5

332000
y

−6

331000
−5.5

−5

330000

178500 179000 179500 180000 180500 181000 181500

x

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.21/27
 Kriging Variance
ordinary kriging variance

−0.5

−0.45
333000

−0.4

−0.35

332000

−0.3
y

−0.25

331000
−0.2

−0.15

330000 −0.1

178500 179000 179500 180000 180500 181000 181500

x

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.22/27
 Prediction levels
Prediction (log zinc)

333000
332000
Pre
d
icti
on

331000
y

330000

178500 179500 180500 181500

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.23/27
 Prediction levels (zinc)

Prediction (zinc)

333000
332000
Pre
dic
onti

331000
y

330000

178500 179500 180500 181500

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.24/27
 Kriging variance

Kriging variance surface

333000
332000
Kr.
var

331000
y

330000

178500 179500 180500 181500

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.25/27
 Notes

• MSPE highest where fewest observations.

• Back-transforming to original zinc levels gives idea
of prediction, but is no longer BLUP (we lose “L”
and “U”).
• gstat does much more complicated analyses too.

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.26/27
 Summary of geostatistical methods

• We’ve really only seen the whirlwind tour of the

basic elements.
• Semivariogram.
• Ordinary kriging.
• Many other options too.
• Questions?

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.27/27
 HW 4

• Repeat this analysis for the Smoky Mountain pH

data.
• Should be able to use the data break code to do
this.
• Due April 28, graphs, output.

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.28/27
 Project 2

• Write up results of HW 4 as a report with

Introduction, Methods, Results, and Discussion
sections.
• Due May 3, 5pm.

Data Break 8: Kriging the Meuse RiverBIOS 737 Spring 2004 – p.29/27