setwd("C:/Users/mpavlic/Documents/2-STUDIJ/9-Analiza/2-v2/Kategoricni kriging")
remove(list=ls())

library(lattice)
library(rgdal)
library(sqldf)
library(gstat)
library(sp)


##### LOAD DATA
DF<-read.csv("zerodist.csv", dec=".", sep=";")
DF$SoilID<-as.factor(DF$SoilID)
str(DF$SoilID)
coordinates(DF)<-~X+Y
zerodist(DF)
DF<-remove.duplicates(DF, zero=0, remove.second = TRUE)
zerodist(DF)
write.table(DF, "zerodist1.csv", sep=";", dec=".")
DF<-read.csv("zerodist1.csv", dec=".", sep=";")
DF$SoilID<-as.factor(DF$SoilID)
str(DF)

##### create table of SoilID
Soil.table<-table(DF [, "SoilID"])
Soil.table

lvls<-levels(DF[,"SoilID"])
print(lvls)


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


##### fit variogram
dat.c <- gstat(dat.c, model = vgm(1, "Cir", 1200, 1), 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<-100
coordinates(DF)<-~X+Y
x.range <- as.integer(range(DF@coords[,1]))
y.range <- as.integer(range(DF@coords[,2]))
str(x.range)

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)<-~x+y
gridded(grd)<-T


##### MAKE SURE THE 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)

##### PLOT THE PREDICTION ---> "Error in parse(text = x) : <text>:1:3: unexpected symbol 1: 0.pred"

print(spplot(z, zcol="0.pred"), col.regions=terrain.colors(64), contour=T, pretty=T, cuts=8, key.space="right")