library(geoR)

data(s100)  

summary(s100) 

plot(s100)


cloud1 <- variog(s100, option = "cloud", max.dist=1)

bin1 <- variog(s100, uvec=seq(0,1,l=11))

	  
par(mfrow=c(2,2))

plot(cloud1, main = "Variogram cloud")
 
plot(bin1, main = "Binned cloud")
 

bin2 <- variog(s100,uvec = seq(0,1,l=11), bin.cloud = T)



plot(bin2, bin.cloud = T, main = "Cloud boxplot")

#Check for significant spatial correlation

 env.mc <- variog.mc.env(s100,obj.var=bin1)

 plot(bin1,envelope = env.mc)

par(mfrow=c(1,1))


plot(bin1)

lines.variomodel(cov.model = "exp", cov.pars = c(1,0.3), nugget = 0, max.dist = 1,  lwd = 3)

smooth <- variog(s100, option = "smooth", max.dist = 1, n.points = 100, kernel = "normal", band = 0.2)

lines(smooth, type ="l", lty = 2)


vario60 <- variog(s100, max.dist = 1, direction=pi/3) 

plot(vario60)

title(main = expression(paste("directional, angle = ", 60 * degree)))
	 
vario.4 <- variog4(s100, max.dist = 1)

plot(vario.4, lwd=2)

####### fit empircal variogram to a parametric form

# Fitting models with nugget fixed to zero
ml <- likfit(s100, ini = c(1,0.5), fix.nugget = T)

wls <- variofit(bin1, ini = c(1,0.5), fix.nugget = T)
	 
# Fitting models with a fixed value for the nugget
ml.fn <- likfit(s100, ini = c(1,0.5), fix.nugget = T, nugget = 0.15)

wls.fn <- variofit(bin1, ini = c(1,0.5), fix.nugget = T, nugget = 0.15)

# Fitting models estimated nugget
ml.n <- likfit(s100, ini = c(1,0.5), nug = 0.5)

wls.n <- variofit(bin1, ini = c(1,0.5), nugget=0.5)

# Now, plotting fitted models against empirical variogram

par(mfrow = c(1,3))

plot(bin1, main = expression(paste("fixed ", tau^2 == 0)))

lines(ml, max.dist = 1)

lines(wls, lty = 2, lwd = 2, max.dist = 1)


plot(bin1, main = expression(paste("fixed ", tau^2 == 0.15)))
lines(ml.fn, max.dist = 1)

lines(wls.fn, lty = 2, lwd = 2, max.dist = 1)

plot(bin1, main = expression(paste("estimated  ", tau^2)))
lines(ml.n, max.dist = 1)

lines(wls.n, lty =2, lwd = 2, max.dist = 1)



summary(ml.n) 


#As a first example consider the prediction at four locations labeled 1, 2, 3, 4 and
#indicated in the figure below.

plot(s100$coords, xlim=c(0,1.2), ylim=c(0,1.2), xlab="Coord X", ylab="Coord Y")

loci <- matrix(c(0.2, 0.6, 0.2, 1.1, 0.2, 0.3, 1.0, 1.1), ncol=2)

text(loci, as.character(1:4), col="red")

polygon(x=c(0,1,1,0), y=c(0,0,1,1), lty=2)
	  
#The command to perform ordinary kriging using the parameters estimated by
#weighted least squares with nugget fixed to zero would be:

kc4 <- krige.conv(s100, locations = loci, krige = krige.control(obj.m = wls))

#Consider now a second example. The goal is to perform prediction on a grid covering
#the area and to display the results. Again, we use ordinary kriging. The commands are:


pred.grid <-  expand.grid(seq(0,1, l=51), seq(0,1, l=51))


kc <- krige.conv(s100, loc = pred.grid, krige = krige.control(obj.m = ml))


image(kc, loc = pred.grid, col=gray(seq(1,0.1,l=30)), xlab="Coord X", ylab="Coord Y")

##Crossvalidation
##Leave-one-out

xv.ml <- xvalid(s100, model = ml) 

par(mfcol = c(5, 2), mar = c(3, 3, 1, 0.5), mgp = c(1.5,0.7, 0)) 

plot(xv.ml)

##Bayesian kriging

bsp4 <- krige.bayes(s100, loc = loci, prior = prior.control(phi.discrete = seq(0,
5, l = 101), phi.prior = "rec"), output = output.control(n.post = 5000))

par(mfrow = c(1, 3), mar = c(3, 3, 1, 0.5), mgp = c(2,1, 0))
hist(bsp4$posterior$sample$beta, main = "", xlab = expression(beta),prob = T)
hist(bsp4$posterior$sample$sigmasq, main = "",xlab = expression(sigma^2), prob = T)
hist(bsp4$posterior$sample$phi, main = "", xlab = expression(phi),prob = T)

##Predictive distributions

par(mfrow = c(2, 2))

for (i in 1:4) {
       kpx <- seq(kc4$pred[i] - 3 * sqrt(kc4$krige.var[i]),
           kc4$pred[i] + 3 * sqrt(kc4$krige.var[i]),
           l = 100)
       kpy <- dnorm(kpx, mean = kc4$pred[i], sd = sqrt(kc4$krige.var[i]))
       bp <- density(bsp4$predic$simul[i, ])
       rx <- range(c(kpx, bp$x))
       ry <- range(c(kpy, bp$y))
       plot(cbind(rx, ry), type = "n", xlab = paste("Location",
           i), ylab = "density", xlim = c(-4, 4),
           ylim = c(0, 1.1))
       lines(kpx, kpy, lty = 2)
       lines(bp)
   }