setwd("C:/Users/mpavlic/Desktop/TEMP")
library(sqldf)
library(gstat)
library(lattice)

DF<-read.csv("temp1.csv", sep=";", dec=".", head=T)
attach(DF)

coordinates(DF) <-~Y+X
spplot(DF, zcol="Globina", col.regions=c("red", "green", "blue", "orange", "black"), key.space="right", cex=sqrt(Globina)/3, pch=20, cuts=seq(0,50, by=10))


# variogram globin vrtin ---------------
### OK without ANIS
v.ok<-variogram(Globina~1, data=DF, cutoff=3600)
plot(v.ok, type="b", plot.numbers=T)

vm.ok<-vgm(40, "Pen", 1300, 15)
(vmf.ok <- fit.variogram(v.ok, vm.ok ))
print(plot(v.ok, plot.numbers = F, pch = 20, col = "darkblue", model = vmf.ok))



##### CREATE GRID --------------------------
### set rages of GRID
x.range <- as.integer(range(DF@coords[,1]))
y.range <- as.integer(range(DF@coords[,2]))

### set resolution of grid and create grid
GridRes<-100
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

### sync names of grid and points
dimnames(coordinates(DF))
dimnames(coordinates(grd))
dimnames(grd@coords) <-list(NULL,c("X", "Y"))
dimnames(grd@bbox) <-list(c("X", "Y"))

### check the plots of grid and point
plot(grd, col="lightgrey")
points(DF, col="red", pch=14, cex=0.33)


##### PREDICTION --------------------------
### OK without anisotropy
z.ok<-krige(Globina~1, DF, grd, vmf.ok)
print(spplot(z.ok, zcol="var1.pred", col.regions=terrain.colors(64),contour=T, pretty=T, cuts=15, key.space="right"), position=c(0,0,.33, 1), more=T)
print(spplot(z.ok,  zcol="var1.var", col.regions=terrain.colors(64), contour=T, pretty=T, cuts=15, key.space="right"))


##### UNIVERSAL KRIGING
### trend surface 2.order
ts2<-krige(Globina~I(Y^2)+I(X^2)+I(Y)+I(X), loc=DF, newdata=grd, model=NULL)
print(spplot(ts1, zcol="var1.pred", contour=T, pretty=T, col.regions=terrain.colors(15)))
summary(ts2)
#linear model ts2
(lm.ts2<-lm(Globina~I(Y^2)+I(X^2)+I(Y)+I(X), data=DF))
summary(lm.ts2)





### variogram of residuals
v.u<-variogram(Globina~I(Y^2)+I(X^2)+I(X)+I(Y), data=DF, cutoff=3600)
plot(v.u2, type="b", plot.numbers=T)
vm.u<-vgm(35, "Cir", 1000, 15)
(vf.u<-fit.variogram(v.u, vm.u))
print(plot(v.u, plot.numbers = F, pch = 20, col = "darkblue", model = vf.u))

### variogram and variogram of residuals

compare.variogram<-data.frame(np = v.ok$np, separation = v.ok$dist, gamma.ok = v.ok$gamma, gamma.uk = v.u$gamma, gamma.dif = v.ok$gamma -v.u$gamma)
print(compare.variogram)

### compare variogram and variogram of residuals
plot(compare.variogram$gamma.ok ~ compare.variogram$separation, pch = 20, 
     col = "blue", type = "b", xlab = "Separation", ylab = "Semivariance",main = "Variogram, Jura Cobalt (ppm)", sub = "OK: blue, UK: green")
points(compare.variogram$gamma.uk ~ compare.variogram$separation, pch = 20, 
       col = "green", type = "b")


### fit variogram of residuals
vr.u <- variogram(Globina ~ I(Y^2)+I(X^2)+I(X)+I(Y), loc = DF, cutoff = 3600)
vrm.u<-vgm(35, "Cir", 1300, 15)
(vrmf.u<-fit.variogram(vr.u, vrm.u))
print(plot(vr.u, plot.numbers = F, pch = 20, col = "darkblue", model = vrmf.u))


### difference OK - UK
(r.range<-vmf.ok$range-vrmf.u$range)
(r.sill<-vmf.ok$sill-vrmf.u$sill)

z.uk<-krige(Globina~I(Y^2)+I(X^2)+I(X)+I(Y), DF, grd, vrmf.u)
print(spplot(z.uk, zcol="var1.pred", col.regions=terrain.colors(64),contour=T, pretty=T, cuts=15, key.space="right"), position=c(0.33, 0, .66, 1), more=T)
print(spplot(z.uk,  zcol="var1.var", col.regions=terrain.colors(64), contour=T, pretty=T, cuts=15, key.space="right"))
summary(z.uk$var1.pred)
summary(DF$Globina)


### plot differences
diff.ok.uk <- data.frame(pred = z.ok$var1.pred - z.uk$var1.pred, 
                         var = z.ok$var1.var - z.uk$var1.var)

str(diff.ok.uk)
coordinates(diff.ok.uk) <- coordinates(z.uk)
print(spplot(diff.ok.uk, zcol = "pred", pretty = T, cuts = 8, col.regions = terrain.colors(64),key.space = "right"))