# load data
data(meuse,meuse.blocks)
# approximation of block variance
# pixel area = 75m x 75m
# exponential covariance function with measurement error = 0, nugget = 0.05,
# part. sill = 0.15 and range parameter = 192.5
preCK <- preCKrige(newdata = meuse.blocks,
model = covmodel("exponential", 0, 0.05, 0.15, scale = 192.5),
pwidth = 75, pheight = 75)
# block prediction by constrained kriging
ck <- CKrige(formula = log(zinc) ~ sqrt(dist),
data = meuse, locations = ~ x + y, object = preCK)
# transforming back to the original scale can be done for CK directly as
# the variance of the ck predictor matches the true block variance
ck@data$Zn <- exp(ck@data$prediction + 0.5 * unlist(preCK@covmat))
# factor: the factor multiplied with the corresponding block mean
# prediction gives the upper bound of the 97.5% confidence interval
# of the prediction
ck@data$factor <- exp(ck@data$prediction + 1.96*ck@data$prediction.se)/ck@data$Zn
# plots with spplot, generic function in the sp package
# plot of the constrained block kriging prediction
# function ck.colors(n) create a vector of n colors
ticks <- seq(0, 1850, by = 185)
spplot(ck, zcol ="Zn" , col.regions = ck.colors(10), at = ticks,
colorkey= list(labels = list(at = ticks, labels = ticks)))
# plot of the factor to get the upper bound of the 97.5% confidence interval
spplot(ck, zcol = "factor", col.regions = ck.colors(10), at = seq(1,3, by = 0.2))
Run the code above in your browser using DataLab