#
# Atelier SFdS: Introduction aux méthodes spatiales et spatio-temporelles
#                 23 et 24 juin 2016, IHP (Paris)
#
# Author: Denis Allard, 24th of July 2016
#
# To be used and circulated freely
################################################################


# -----------------
#
# 1ere partie:  Utilisation de RGeostats
#
# -----------------

library(RGeostats)

# Lecture des données
jura = read.table("http://ciam.inra.fr/biosp/sites/ciam.inra.fr.biosp/files/jura_tout.txt",header=T)
Ni   = jura$Ni
Co   = jura$Co
rt = jura$rt
lu = jura$lu

# Visualisation
par(mfrow=c(1,2))
plot(jura$x,jura$y,main="Ni",xlab="X",ylab="Y",
     pch=19,cex=Ni/20,col="blue")
plot(jura$x,jura$y,main="Co",xlab="X",ylab="Y",
     pch=19,cex=Co/10,col="red")

par(mfrow=c(1,2))
plot(Ni,Co,main = "Co vs. Ni",pch=19)
plot(log(Ni),log(Co),main = "log(Co) vs. log(Ni)",pch=19)


jura.db = db.create(x1=jura$x,x2=jura$y,z1=jura$Co,z2=jura$Ni,z3=jura$Zn)
jura.db = db.rename(jura.db,name=2,newname="X")
jura.db = db.rename(jura.db,name=3,newname="Y")
jura.db = db.rename(jura.db,name=4,newname="Co")
jura.db = db.rename(jura.db,name=5,newname="Ni")
jura.db = db.rename(jura.db,name=6,newname="Zn")
jura.db

# Covariance omindirectionnelle
jura.cov = vario.calc(jura.db,lag=0.2,nlag=10,calcul="cov")
plot(jura.cov,npairdw=TRUE,title="Covariance (Co,Ni,Zn)")

# Variogramme omnidirectionnel
jura.vario = vario.calc(jura.db,lag=0.2,nlag=10,calcul="vg")
plot(jura.vario,npairdw=TRUE,title="Variogramme (Co,Ni,Zn)")

# Variogramme 4 directions
jura.vario4d = vario.calc(jura.db,lag=0.2,nlag=10,calcul="vg",dir=c(0,45,90,135))
plot(jura.vario4d,npairdw=TRUE,title="Variogramme (Co,Ni,Zn)")

# Ajustements variogrammes
melem.name()
jura.model1 = model.auto(jura.vario,struct=melem.name(c(2)),title="Exp",verbose=1) # exp
jura.model2 = model.auto(jura.vario,struct=melem.name(c(8)),title="Matern",verbose=1) # exp

jura.model3 = model.auto(jura.vario,struct=melem.name(c(1,2,2)),title="Nug + Exp + Exp",verbose=1) # nugget + exp + exp
jura.model4 = model.auto(jura.vario,struct=melem.name(c(1,3,3)),title="Nug + Sph + Sph",verbose=1) # nugget + sph + sph
jura.model5 = model.auto(jura.vario,struct=melem.name(c(1,8,8)),title="Nug + Matern + Matern",verbose=1) # nugget + Matern + Matern                       



# Krigeage
grille = db.create(flag.grid=T,nx=c(108,108),x0=c(0.4,0.4),dx=c(0.05,0.05))
unique.neigh=neigh.init(type=0,ndim=2)

# Krigeage homotopique
resultm1 = kriging(jura.db,grille,model=jura.model1,neigh=unique.neigh,uc=c("1"))
plot(resultm1,name=4,pos.legend=5,col=topo.colors(200))
plot(resultm1,name=7,pos.legend=5)
plot(resultm1,name=5,pos.legend=5,col=topo.colors(200))
plot(resultm1,name=8,pos.legend=5)
plot(resultm1,name=6,pos.legend=5,col=topo.colors(200))
plot(resultm1,name=9,pos.legend=5)


resultm3 = kriging(jura.db,grille,model=jura.model3,neigh=unique.neigh)
plot(resultm3,name=4,pos.legend=5,col=topo.colors(200))
plot(resultm3,name=7,pos.legend=5)
plot(resultm3,name=5,pos.legend=5,col=topo.colors(200))
plot(resultm3,name=8,pos.legend=5)
plot(resultm3,name=6,pos.legend=5,col=topo.colors(200))
plot(resultm3,name=9,pos.legend=5)
xvalid(jura.db,model=jura.model3,neigh=unique.neigh,uc=c("1"))

# Essai de validation

idx.data = sample(dim(jura)[1])[1:159]
plot(jura$x[idx.data],jura$y[idx.data],main="Ni",xlab="X",ylab="Y",
     pch=19,cex=Ni/20,col="blue")
points(jura$x[-idx.data],jura$y[-idx.data],main="Ni",xlab="X",ylab="Y",
     pch=19,cex=Ni/20,col="red")
jura.data  = jura[idx.data,]
jura.valid = jura[-idx.data,]
jura.data.db  = db.create(x1=jura.data$x,x2=jura.data$y,z1=jura.data$Co,z2=jura.data$Ni,z3=jura.data$Zn)
jura.valid.db = db.create(x1=jura.valid$x,x2=jura.valid$y,z1=jura.valid$Co,z2=jura.valid$Ni,z3=jura.valid$Zn)

# Validation de co-krigeage
grille.valid <- db.create(x1=jura.valid$x,x2=jura.valid$y)
krigeage = kriging(jura.data.db,grille.valid,model=jura.model3,neigh=unique.neigh,uc=c("1"))
rmse.Co = sqrt(mean((krigeage[,4]-jura.valid.db[,4])^2))
rmsse.Co = sqrt(mean((krigeage[,4]-jura.valid.db[,4])^2/krigeage[,7]^2))
rmse.Ni = sqrt(mean((krigeage[,5]-jura.valid.db[,5])^2))
rmsse.Ni = sqrt(mean((krigeage[,5]-jura.valid.db[,5])^2/krigeage[,8]^2))
rmse.Zn = sqrt(mean((krigeage[,6]-jura.valid.db[,6])^2))
rmsse.Zn = sqrt(mean((krigeage[,6]-jura.valid.db[,6])^2/krigeage[,9]^2))


# -----------------
#
# Utilisation de gstat
#
# -----------------


library(gstat)
library(sp)



# Visualisation
jura.sp = SpatialPointsDataFrame(jura[,1:2],jura[,c(6,9,11)])

par(mfrow=c(1,2))
bubble(jura.sp,"Co")
bubble(jura.sp,"Ni",col="red")

vario.Ni = variogram(Ni ~ 1 ,jura.sp)
plot(vario.Ni,xlim=c(-0.05,2))

vario.Co = variogram(Co ~1 ,jura.sp)
plot(vario.Co,xlim=c(-0.05,2))

vario.Zn = variogram(Zn ~1 ,jura.sp)
plot(vario.Zn,xlim=c(-0.05,2))

par(mfrow=c(1,2))
hist(Zn);hist(log(Zn))
vario.Zn = variogram(log(Zn) ~1 ,jura.sp)
plot(vario.Zn,xlim=c(-0.05,2))

boxplot(Ni ~ rt)
vario.Nirt = variogram(Ni ~ 1+rt,jura.sp)
plot(vario.Nirt)

v.fit1 = fit.variogram(vario.Ni, model = vgm(model="Exp", psill= 80, range= 0.5, nugget=1))
v.fit1
plot(vario.Ni,model=v.fit1)
attr(v.fit1,"SSErr")
v.fit2 = fit.variogram(vario.Ni, model = vgm(model="Mat", psill= 80, range= 0.5,kappa=0.35,nugget=1))
v.fit2
plot(vario.Ni,model=v.fit2)
attr(v.fit2,"SSErr")                      
v.fit3 = fit.variogram(vario.Ni, model = vgm(model="Mat", psill= 80, range= 0.5,kappa=0.35,nugget=1),
                                            fit.kappa=TRUE) # Essayer d'autres valeurs de kappa

v.fit3 = fit.variogram(vario.Ni, model = vgm(model="Mat", psill= 74, range= 0.73,kappa=0.5,nugget=12),
                       fit.kappa=TRUE) 
v.fit3
plot(vario.Ni,model=v.fit3)
attr(v.fit3,"SSErr") 



plot(vario.Ni,model=v.fit1)
plot(vario.Ni,model=v.fit3)


jura.multi = gstat(NULL, "Ni", Ni ~ 1, jura.sp)
jura.multi = gstat(jura.multi, "Co", Co ~ 1, jura.sp)
jura.multi = gstat(jura.multi, "log(Zn)", log(Zn) ~ 1, jura.sp)
jura.multi
vario.multi = variogram(jura.multi)
plot(vario.multi)


g = gstat(jura.multi, model = vgm(model="Exp", psill= 80, range= 0.7,nugget=12), fill.all=TRUE)
v.fit5 = fit.lmc(vario.multi,g)
plot(vario.multi,model=v.fit5)
v.fit5

out = gstat.cv(v.fit5)
summary(out)
mean(out$residual)
mean(out$residual^2)
v.fit5$model$Ni$psill[1]/(v.fit5$model$Ni$psill[1] + v.fit5$model$Ni$psill[2])
mean(out$zscore^2)
plot(Ni,out$Ni.pred)
abline(0,1)  # On observe une sous-estimation des valeurs les plus élevées
plot(Ni,out$residual)
abline(h=0)  # Cette observation est encore plus visible ici
plot(out$Ni.pred,out$residual)
abline(h=0)  # Pas de corrélation entre valeurs prédites et erreur de prédiction