library(lattice)
library(rgdal)
library(gstat)

#####CREATE DATASET
n<-100
DF <- data.frame(X = rnorm(n, mean=100, sd=200), Y = rnorm(n, mean=500, sd=500), 
                 SoilFactor = c("nic", "ena","dva","tri", "stiri"))


print(xyplot(X~Y, data=DF, panel=function(x, y,...){
  ltext(x, y, labels=DF[,"SoilFactor"], cex=0.7, col="red")}))

## create table
Soil.table<-table(DF [, "SoilFactor"])
Soil.table
tmp<-as.vector(DF[,"SoilFactor"])

DF[,"SoilFactor"] <-factor(tmp)
lvls<-levels(DF[,"SoilFactor"])
print(lvls)

### CREATE GSTAT OBJECT -01 classes
nmx<-10
dat.c <- gstat(id = lvls[1], formula = I(SoilFactor == lvls[1]) ~ 1, loc=~X+Y, data = DF, nmax = nmx, set = list(order = 2, zero = 0.01))
dat.c <- gstat(dat.c, lvls[2], I(SoilFactor == lvls[2]) ~ 1, ~X + Y, DF, nmax = nmx)
dat.c <- gstat(dat.c, lvls[3], I(SoilFactor == lvls[3]) ~ 1, ~X + Y, DF, nmax = nmx)
dat.c <- gstat(dat.c, lvls[4], I(SoilFactor == lvls[4]) ~ 1, ~X + Y, DF, nmax = nmx)
dat.c <- gstat(dat.c, lvls[5], I(SoilFactor == lvls[5]) ~ 1, ~X + Y, DF, nmax = nmx)


#fit variogram
plot(variogram(dat.c))
dat.c <- gstat(dat.c, model = vgm(.2, "Cir", 400, 0), fill.all = T)
x <- variogram(dat.c)

# fit lmc
dat.fit = fit.lmc(x, dat.c, fit.ranges=F, fit.sills=c(T,T,T))
plot(x, model = dat.fit)



####CREATE PREDICTION GRID
GridRes<-10
coordinates(DF)<-~X+Y
x.range <- as.integer(range(DF@coords[,1]))
y.range <- as.integer(range(DF@coords[,2]))
grd<-expand.grid(x=seq(from=x.range[1], to=x.range[2], by=GridRes),
y=seq(from=y.range[1], to=y.range[2], by=GridRes))
coordinates(grd)<-~y+x
gridded(grd)<-T
plot(grd)
#CHECK THAT GRID AND POINTS HAVE SAME COORDINATE NAMES
dimnames(coordinates(DF))
dimnames(coordinates(grd))
dimnames(grd@coords) <-list(NULL,c("X", "Y"))
dimnames(grd@bbox) <-list(c("X", "Y"))

#PREDICT
z<-predict(dat.fit, newdata=grd)
str(z)

print(spplot(z, z=seq(1,10, by=2), main="Prediction", col.regions=heat.colors(64)), split=c(1,1,1,2), more=T)
print(spplot(z, z=seq(2,19, by=2), main="Variance",col.regions=heat.colors(64)), split=c(1,2,1,2))