# 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,192.5),pwidth=75,pheight=75)
# block prediction by constrained kriging on the log scale
CK=CKrige(formula=log(zinc)~sqrt(dist),data=meuse,
locations=~x+y,object=preCK,ex.out=TRUE)
# backtransformation to the original scale for the CK prediction
beta=CK$parameter$beta.coef
M=meuse.blocks@data$M
c1 <- 0.2
c2 <- beta[2]^2 * meuse.blocks@data$M
CK$object@data$Zn=exp(CK$object@data$prediction
+ 0.5*(0.2+beta[2]^2*M-unlist(preCK@covmat)))
# U: upper limits of the relative 95# U multiplied by the predictions CK$object@data$Zn gives
# the upper limits of the 95
CK$object@data$U=exp(CK$object@data$prediction
+1.96*CK$object@data$prediction.se) /CK$object@data$Zn
# plots with spplot, generic function in the sp package
# the plot shows the constrained kriging predictions on
# the orginal scale
# function ck.colors(n) create a rainbow-like color vector
breaks <- seq(0, 1850, by = 185)
spplot(CK$object,zcol="Zn",at=breaks,col.regions=ck.colors(10),
colorkey=list(labels=list(at=breaks,labels=breaks)))
# plot of the factor to get the upper bound of the
97.5
breaks=seq(1,3.2,by=0.2)
spplot(CK$object,zcol="U",at=breaks,col.regions=ck.colors(11),
colorkey=list(labels=list(at=breaks,labels=breaks)))
Run the code above in your browser using DataLab