date()
## [1] "Sun Apr 21 18:54:44 2013"
Krige the Parana rainfall.
1.Examine the data for spatial trends and normality.
suppressWarnings(require(geoR, quietly = TRUE))
## -------------------------------------------------------------- Analysis of
## geostatistical data For an Introduction to geoR go to
## http://www.leg.ufpr.br/geoR geoR version 1.7-4 (built on 2012-06-29) is
## now loaded --------------------------------------------------------------
data(parana)
class(parana)
## [1] "geodata"
summary(parana)
## $n
## [1] 143
##
## $coords.summary
## east north
## min 150.1 70.36
## max 768.5 461.97
##
## $distances.summary
## min max
## 1.0 619.5
##
## $data.summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 163 234 270 274 318 414
##
## $borders.summary
## east north
## min 138.0 46.77
## max 798.6 507.93
##
## $others
## [1] "loci.paper"
##
## attr(,"class")
## [1] "summary.geodata"
plot(parana)
plot(parana, trend = "1st")
With 1st trend removed the variation of values appears to be roughly symmetric.
2.Compute empirical variograms.
par(mfrow = c(2, 1))
plot(variog(parana, trend = "1st", max.dist = 309.7), main = "Classical")
## variog: computing omnidirectional variogram
plot(variog(parana, trend = "1st", max.dist = 309.7, estimator.type = "modulus"),
main = "Modulus")
## variog: computing omnidirectional variogram
3.Fit a variogram model to empirical variogram.
par(mfrow = c(1, 1))
plot(variog(parana, trend = "1st", max.dist = 309.7))
## variog: computing omnidirectional variogram
p.vm = likfit(parana, trend = "1st", ini = c(700, 150), nug = 350, cov.model = "exponential")
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
lines(p.vm, col = "blue")
# p.vm1 = likfit(parana, trend = '1st', ini = c(700, 150), nug = 350,
# cov.model = 'spherical') lines(p.vm, col = 'red')
4.Construct an interpolated surface using a 200 by 200 grid.
plot(parana$coords, pch = 10, col = "red", xlab = "east", ylab = "north")
pred.grid = expand.grid(seq(150, 769, l = 200), seq(70, 462, l = 200))
points(pred.grid, pch = 15, cex = 0.18)
p.ks = krige.conv(parana, loc = pred.grid, krige = krige.control(trend.d = "1st",
trend.l = "1st", obj.m = p.vm))
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with mean given by a 1st order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
str(p.ks)
## List of 6
## $ predict : num [1:30414] 362 361 360 359 359 ...
## $ krige.var : num [1:30414] 530 523 513 501 497 ...
## $ beta.est : Named num [1:3] 419.452 -0.144 -0.396
## ..- attr(*, "names")= chr [1:3] "beta0" "beta1" "beta2"
## $ distribution: chr "normal"
## $ message : chr "krige.conv: Kriging performed using global neighbourhood"
## $ call : language krige.conv(geodata = parana, locations = pred.grid, krige = krige.control(trend.d = "1st", trend.l = "1st", obj.m = p.vm))
## - attr(*, "sp.dim")= chr "2d"
## - attr(*, "prediction.locations")= symbol pred.grid
## - attr(*, "parent.env")=<environment: R_GlobalEnv>
## - attr(*, "data.locations")= language parana$coords
## - attr(*, "class")= chr "kriging"
image(p.ks, col = terrain.colors(20), xlab = "east", ylab = "north")
text(parana$coords, labels = parana$data, cex = 0.5)
contour(p.ks, add = TRUE, nlev = 20)
# image(p.ks, val = sqrt(p.ks$krige.var), col = heat.colors(16), xlab =
# 'east', ylab = 'north') points(parana$coords, pch = 20) contour(p.ks,
# val = sqrt(p.ks$krige.var), add = TRUE, nlev = 16)
5.Examine the prediction errors.
xv = xvalid(parana, model = p.vm)
## xvalid: number of data locations = 143
## xvalid: number of validation locations = 143
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143,
## xvalid: end of cross-validation
plot(xv$data, xv$predicted, xlab = "Observed Parana rainfall", ylab = "Predicted Parana rainfall at observed location",
asp = 1, col = "red", pch = 20)
abline(0, 1)
abline(lm(xv$predicted ~ xv$data), col = "red")
mean(xv$error^2)
## [1] 527.2
mean(abs(xv$error))
## [1] 17.74