date()
## [1] "Wed Apr 24 13:58:38 2013"
Krige the Parana rainfall.
suppressMessages(install.packages("knitr"))
## Error: trying to use CRAN without setting a mirror
suppressMessages(require(knitr))
suppressMessages(install.packages("geoR"))
## Error: trying to use CRAN without setting a mirror
suppressMessages(require(geoR))
suppressMessages(install.packages("maps"))
## Error: trying to use CRAN without setting a mirror
suppressMessages(require(maps))
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")
plot(parana, trend = "2nd")
The maximum pairwise distance is 619.5 degrees 2. Compute empirical variograms.
plot(variog4(parana, trend = "2nd", max.dist = 309.75), omni = TRUE)
## 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
iv = c(600, 100)
summary(likfit(parana, ini = iv, cov.model = "exp", fix.nug = TRUE, message = FALSE))$likelihood$AIC
## [1] 1394
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "1st", fix.nug = TRUE,
message = FALSE))$likelihood$AIC
## [1] 1365
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = TRUE,
message = FALSE))$likelihood$AIC
## [1] 1347
summary(likfit(parana, ini = iv, cov.model = "exp", trend = "2nd", fix.nug = FALSE,
message = FALSE))$likelihood$AIC
## [1] 1338
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = TRUE,
message = FALSE))$likelihood$AIC
## [1] 1350
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE,
message = FALSE))$likelihood$AIC
## [1] 1337
summary(likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE))
## 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.
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta0 beta1 beta2 beta3 beta4 beta5
## 424.945 0.075 -0.637 0.000 0.000 0.001
##
## Parameters of the spatial component:
## correlation function: spherical
## (estimated) variance parameter sigmasq (partial sill) = 329
## (estimated) cor. fct. parameter phi (range parameter) = 159
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 403
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 1 (no transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 159.4
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "-660" "9" "1337" "1364"
##
## non spatial model:
## log.L n.params AIC BIC
## "-671" "7" "1356" "1376"
##
## Call:
## likfit(geodata = parana, trend = "2nd", ini.cov.pars = iv, fix.nugget = FALSE,
## cov.model = "sph")
plot(variog(parana, trend = "2nd", uvec = seq(0, 309.75, l = 29)))
## variog: computing omnidirectional variogram
lines.variomodel(cov.model = "sph", cov.pars = c(328.9, 159.4), nug = 403, max.dist = 309.75,
col = "green")
t.modelS = likfit(parana, ini = iv, cov.model = "sph", trend = "2nd", fix.nug = FALSE,
messages = TRUE)
## 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.
loc.grd = expand.grid(seq(150.122, 768.5087, l = 200), seq(70.36, 461.9581,
l = 200))
kcE = krige.conv(parana, loc = loc.grd, krige = krige.control(trend.d = "2nd",
trend.l = "2nd", obj.m = t.modelS))
## 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
kcS = krige.conv(parana, loc = loc.grd, krige = krige.control(trend.d = "2nd",
trend.l = "2nd", obj.m = t.modelS))
## 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(kcE, col = rev(heat.colors(20)), xlab = "East", ylab = "North")
title(main = "Spherical Model w/ Nugget = 381.2")
map("state", add = TRUE)
contour(kcE, add = TRUE, nlevel = 20)
xv = xvalid(parana, model = t.modelS)
## 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 Rainfall", ylab = "Predicted Rainfall",
asp = 1, col = "red", pch = 20)
title(main = "Spherical Model")
abline(0, 1)
abline(lm(xv$predicted ~ xv$data), col = "red")
mean(xv$error^2)
## [1] 534.3
mean(abs(xv$error))
## [1] 17.89