date()
## [1] "Wed Apr 23 02:08:38 2014"
The objective of this problem set is to have you demonstrate skill at using functions in the geoR (or gstat) package for creating an interpolated surface using the methods of kriging.
require(geoR)
## Loading required package: geoR
## Loading required package: sp
## Loading required package: MASS
## --------------------------------------------------------------
## 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
## --------------------------------------------------------------
require(ggplot2)
## Loading required package: ggplot2
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, bor = borders)
(1.) Examine the data for spatial trends and normality.
plot(parana, trend = "1st")
plot(parana, tren = "2nd")
(2.) Compute empirical variograms on the residuals.
plot(variog(parana, trend = "2nd", max.dist = 309.75))
## variog: computing omnidirectional variogram
plot(variog4(parana, trend = "2nd", max.dist = 309.75))
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
(3.) Fit a variogram model.
plot(variog(parana, trend = "2nd", uvec = seq(0, 309.7462, l = 29)))
## variog: computing omnidirectional variogram
lines.variomodel(cov.model = "sph", cov.pars = c(360.8187, 214.5301), nug = 419.9149,
max.dist = 309.7462, col = "blue")
p.vm = likfit(parana, trend = "2nd", ini = c(360.8187, 214.5301), nug = 419.9149,
cov.model = "spherical")
## kappa not used for the spherical 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.
p.vm
## likfit: estimated model parameters:
## beta0 beta1 beta2 beta3 beta4 beta5
## "435.4928" " 0.0477" " -0.7009" " -0.0004" " 0.0001" " 0.0006"
## tausq sigmasq phi
## "419.9155" "360.8170" "214.5301"
## Practical Range with cor=0.05 for asymptotic range: 214.5
##
## likfit: maximised log-likelihood = -660
(4.) Construct an interpolated surface using a 200 by 200 grid.
plot(parana$coords, pch = 10, col = "blue", 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 = "2nd",
trend.l = "2ndst", 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 2nd order polynomial on the coordinates
## krige.conv: Kriging performed using global neighbourhood
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)
(5.) Examine the prediction errors.
xv.wk = 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
df = data.frame(x = xv.wk$data, y = xv.wk$predicted)
ggplot(df, aes(x, y)) + geom_point() + geom_abline(intercept = 0, slope = 1) +
geom_smooth(method = lm, color = "green")
mean(xv.wk$error^2)
## [1] 536.6
mean(abs(xv.wk$error))
## [1] 17.87