##
## Fundamentos de geoestatística com o pacote geoR
##

## importando dados
dados <- read.table("http://www.leg.ufpr.br/~paulojus/CE097/dados/dadosIntro.txt", head=T)
dados
head(dados)
str(dados)

## instalando o pacote geoR (se necessário) e suas dependencias
#install.packages("geoR", dep=TRUE)

## carregando o pacote geoR
require(geoR)
## convertendo para objeto da classe geodata 
dG <- as.geodata(dados, coords.col=c(1,2), data.col=3)
class(dG)

## nomes de conjuntos de dados disponíveis no pacote
data(package="geoR")

## visualizando os dados
points(dG)
args(points.geodata)
points.geodata(dG, pt.div="quint")
points.geodata(dG, pt.div="quint", cex.min=1, cex.max=1)

## Um resumo dos dados
summary(dG)
summary(dG, lambda=0) ## trsnformação log (Box-Cox)  
dG	
## uma outra visualização dos dados
plot(dG)
args(plot.geodata)
plot(dG, lambda=0)
plot(dG, lowess=TRUE)

## outras visualizações com points.geodata()
points(dG)
points(dG, lambda=0)
points(dG, pt.div="quint", cex.min=1, cex.max=1)
?points.geodata
points(dG, pt.div="quint", cex.min=1, cex.max=1,permute=TRUE)

## obtendo um variograma (empírico)
v <- variog(dG)	
plot(v)		

## outras opções de variograma
vC <- variog(dG, option="cloud")	
vC <- variog(dG, max.dist=100, option="cloud")	
plot(vC)		

par(mfrow=c(2,2))
v <- variog(dG, max.dist=100)	
plot(v)		
v <- variog(dG, max.dist=80)	
plot(v)		
v <- variog(dG, uvec=seq(0, 80, by=10))	
plot(v)		
v <- variog(dG, uvec=seq(0, 80, by=5))	
plot(v)		
par(mfrow=c(1,1))

## escolhendo um deles
v <- variog(dG, uvec=seq(0, 80, by=10))	
plot(v)		

## envelope do variograma obtido sob permutação dos dados
v.mc <- variog.mc.env(dG, obj.variog=v)
plot(v, env=v.mc)

## "ajustando" função de variograma (visualmente)
v <- variog(dG, max.dist=80)	
plot(v)		
ef <- eyefit(v)
ef

## definindo um "grid" de predição
gr <- expand.grid(seq(0,100, len=51), seq(0,100, len=51))
## visualizando a malha na qual será feita a predição
points(dG)
points(gr, pch=19, cex=0.25, col=2)	

kc.ef <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef))

## visualizando os valores preditos na forma de um mapa (com diferentes opções de cores)
x11()
par(mfrow=c(2,2), mar=c(2,2,1,1))
image(kc.ef)
image(kc.ef, col=terrain.colors(21))
image(kc.ef, col=gray(seq(1,0, l=21)))
image(kc.ef, col=gray(seq(1,0, l=5)))
par(mfrow=c(1,1))

## mais sobre palhetas de cores
hcl.pals()
image(kc.ef, col=hcl.colors(20, "viridis"))
image(kc.ef, col=hcl.colors(20, "viridis", rev = TRUE))

## para ver e escolher...
par(mfcol=c(2,3), mar=c(2,2,3,1))
for(i in hcl.pals()){
    image(kc.ef, col=hcl.colors(20, i), main = i)
    image(kc.ef, col=hcl.colors(20, i, rev = TRUE), main = paste(i, "rev"))
    Sys.sleep(3)
}
par(mfrow=c(1,1))

## Comentário: paletas "robustas"
par(mfrow=c(2,3), mar=c(2,2,3,1))
image(kc.ef, col=hcl.colors(20, "Blues 3", rev = TRUE))
image(kc.ef, col=hcl.colors(20, "YlGnBu", rev = TRUE))
image(kc.ef, col=hcl.colors(20, "Viridis", rev = TRUE))
##
image(kc.ef, col=hcl.colors(20, "Purple-Green", rev = TRUE))
image(kc.ef, col=hcl.colors(20, "Blue-Red 3"))
par(mfrow=c(1,1))


summary(kc.ef$predict)
image(kc.ef, col=terrain.colors(4), breaks=c(0, 30, 80, 100, 150))
image(kc.ef, val = sqrt(kc.ef$krige.var), coords.data=dG$coords)
##
## Fim da análise inicial/básica
##

##
## Podemos revisitar passos e usar outros métodos 
##
plot(v)		
ef <- eyefit(v)  ## ajustar e salvar 3 escolhas diferentes de modelos
ef


vf1 <- variofit(v, ini=ef[[1]])
vf1
vf2 <- variofit(v, ini=ef[[2]])
vf2
vf3 <- variofit(v, ini=ef[[3]])
vf3
## ?variofit  para mais opções e sobre a função

par(mfrow=c(1,3))
plot(v)
lines(ef[[1]], col=2)
lines(vf1, col=4)

plot(v)
lines(ef[[2]], col=2)
lines(vf2, col=4)

plot(v)
lines(ef[[3]], col=2)
lines(vf3, col=4)
par(mfrow=c(1,1))

## outro envelope: de model ajustado
vf1.env <- variog.model.env(dG, obj.variog = v, model.pars = vf1)
plot(v)
lines(vf1)
lines(vf1.env)

## estimando a media
ONES <- rep(1, nrow(dados))
D <- as.matrix(dist(dG$coords))
names(vf1)
vf1$cov.pars
## supondo que vá usar o modelo exponencial
SIGMA <- with(vf1, nugget * diag(nrow(dados)) + cov.pars[1] * exp(-D/cov.pars[2]))
dim(SIGMA)
## $\hat{\mu} = (1'\Sigma^{-1}1)^{-1} 1'\Sigma^{-1}y$
iS.ONES <- solve(SIGMA, ONES)
(mu.gls <- solve(crossprod(iS.ONES, ONES), crossprod(iS.ONES, dG$data)))
## ou ...
(mu.gls <- sum(iS.ONES * dG$data)/sum(iS.ONES))
mean(dG$data)

rm(ONES,D,SIGMA,iS.ONES)

## Fazendo Função para calcular média a partir de um modelo ajustado
names(vf1)
args(cov.spatial)
args(varcov.spatial)

mu.fun <- function(fit, geodata){
    ONES <- rep(1, nrow(geodata$coords))
    SIG <- with(fit,
                  varcov.spatial(coords = geodata$coords,
                                 cov.model = cov.model,
                                 nugget = nugget,
                                 cov.pars = cov.pars))$varcov
    iSIG.ONES <- solve(SIG, ONES)
    return(sum(iSIG.ONES * dG$data)/sum(iSIG.ONES))
}
mu.fun(vf1, dG)
mu.fun(vf2, dG)
mu.fun(vf3, dG)
mean(dG$data)

## estimativa por máxima verossimilhança
## usando a funcao da geoR
ml <- likfit(dG, ini=c(300, 20))
ml
mu.fun(ml, dG)

## calculando a verossimilhanca de um ajuste variofit
loglik.GRF(dG, obj.model = ml)
loglik.GRF(dG, obj.model = vf1)
loglik.GRF(dG, obj.model = vf2)
loglik.GRF(dG, obj.model = vf3)

## voltando a variogramas
plot(v)
lines(ml, col="blue")
lines(vf1, col="red")
lines(vf2, col="green")
lines(vf3, col="black")
legend("topleft", legend=c("ml","vf1","vf2","vf3"),
       col=c("blue","red","green","black"), lty=1)

## definindo um "grid" de predição
gr <- expand.grid(seq(0,100, len=51), seq(0,100, len=51))
points(dG)
points(gr, pch=19, cex=0.25, col=2)	
	

## fazendo a predição espacial (krigagem)

## com opções de ajusta "a olho"(eyefit)
kc.ef1 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[1]]))
kc.ef2 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[2]]))
kc.ef3 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ef[[3]]))

## com um variograma ajustado (variofit)
kc.vf1 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=vf1))
kc.vf2 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=vf2))
kc.vf3 <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=vf3))

## com parâmetros ajustados por verossimilhança (likfit)
kc.ml <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ml))
names(kc.ml)

## visualizando os valores preditos na forma de um mapa
par(mfrow=c(2,3), mar=c(2,2,1,1))
image(kc.ef1, col=terrain.colors(21))
image(kc.ef2, col=terrain.colors(21))
image(kc.ef3, col=terrain.colors(21))
image(kc.vf1, col=terrain.colors(21))
image(kc.vf2, col=terrain.colors(21))
image(kc.ml, col=terrain.colors(21))
par(mfrow=c(1,1), mar=c(2,2,1,1))

## escala comum
ZL <- range(
    range(kc.ef1$pred),
    range(kc.ef2$pred),
    range(kc.ef3$pred),
    range(kc.vf1$pred),
    range(kc.vf2$pred),
    range(kc.ml$pred)
)
par(mfrow=c(2,3), mar=c(2,2,1,1))
image(kc.ef1, col=terrain.colors(21), zlim=ZL)
image(kc.ef2, col=terrain.colors(21), zlim=ZL)
image(kc.ef3, col=terrain.colors(21), zlim=ZL)
image(kc.vf1, col=terrain.colors(21), zlim=ZL)
image(kc.vf2, col=terrain.colors(21), zlim=ZL)
image(kc.ml, col=terrain.colors(21), zlim=ZL)
par(mfrow=c(1,1))

## visualizando os erros padrão de predição (mapa de incerteza de predição)
ZL <- range(
    range(sqrt(kc.ef1$krige.var)),
    range(sqrt(kc.ef2$krige.var)),
    range(sqrt(kc.ef3$krige.var)),
    range(sqrt(kc.vf1$krige.var)),
    range(sqrt(kc.vf2$krige.var)),
    range(sqrt(kc.ml$krige.var))
)
par(mfrow=c(2,3), mar=c(2,2,1,1))
image(kc.ef1, val=sqrt(kc.ef1$krige.var), coords.data=dG$coords, zlim=ZL)
image(kc.ef2, val=sqrt(kc.ef2$krige.var), coords.data=dG$coords, zlim=ZL)
image(kc.ef3, val=sqrt(kc.ef3$krige.var), coords.data=dG$coords, zlim=ZL)
image(kc.vf1, val=sqrt(kc.vf1$krige.var), coords.data=dG$coords, zlim=ZL)
image(kc.vf2, val=sqrt(kc.vf2$krige.var), coords.data=dG$coords, zlim=ZL)
image(kc.ml, val=sqrt(kc.ml$krige.var), coords.data=dG$coords, zlim=ZL)
par(mfrow=c(1,1))

par(mfrow=c(1,3))
plot(kc.ml$pred, kc.vf1$pred, asp = 1); abline(0,1)
plot(kc.ml$pred, kc.vf2$pred, asp = 1); abline(0,1)
plot(kc.vf1$pred, kc.vf2$pred, asp = 1); abline(0,1)
par(mfrow=c(1,1))


## outras opções de visualização
par(mfrow=c(2,2))
contour(kc.ml)
image(kc.ml, col=terrain.colors(21))
points(dG, add=T)
contour(kc.ml, add=T, nlev=21)
persp(kc.ml)
persp(kc.ml, theta=20, phi=35)
par(mfrow=c(1,1))

##
## exportando predições
##
## exportando somente um vetor dos atributos (formato texto)
write(kc.ml$pred, file="kc1.txt", ncol=1)
file.show("kc.ml.txt")    ## tecle c para terminar a exibição do arquivo

## exportando somente os atributos como matriz (formato texto)
write(kc.ml$pred, ncol=51, file="kc2.txt")

## exportando atributos e locações de predição
names(kc.ml)
write.table(cbind(gr, data.frame(kc.ml[c(1,2)])),file="kc3.txt",
            row.names = FALSE)
## ou...
write.csv(cbind(gr, data.frame(kc.ml[c(1,2)])),file="kc4.csv",
            row.names = FALSE)
## ou...
write.csv2(cbind(gr, data.frame(kc.ml[c(1,2)])),file="kc5.csv",
            row.names = FALSE)

##
## simulacoes condicionais
##
args(krige.conv)
args(output.control)

kc.ml <- krige.conv(dG, loc=gr, 
                    krige=krige.control(obj.m=ml),
                    output = output.control(n.pred=1000))
names(kc.ml)
dim(kc.ml$simulations)

## a krigagem é aproximada por simulação condicional
predSim <- apply(kc.ml$simulations, 1, mean)
plot(kc.ml$pred, predSim); abline(0,1)
summary(kc.ml$pred - predSim)


## mapa de prob de estar acima de 65
limiar <- apply(kc.ml$simulations, 1, function(x) mean(x>65))
image(kc.ml, val=limiar)
contour(kc.ml, val=limiar)

## suponha agora que definimos como regiao de alto risco
## as quye possuiram prob > 0.7 de estar acima de 65
limiar.alto <- ifelse(limiar < 0.7, 0, 1)
image(kc.ml, val=limiar.alto, col=c(1,0))

## agora definindo 5 niveis de risco em fc da prob
image(kc.ml, val=limiar, breaks=seq(0, 1, by=0.2),
      col=c("blue","green","yellow","red","black"))

##
## definindo um polígono de interesse
points(dG)
bor <- locator(type="l")
polygon(bor, col=4)

## completar aqui a krigagem dentro do polygono


plot(dG)
v <- variog(dG, max.dist=80)
plot(v)

args(variog)
v0 <- variog(dG, max.dist=80, direction=0)
v45 <- variog(dG, max.dist=80, direction=pi/4)
v90 <- variog(dG, max.dist=80, direction=pi/2)
v135 <- variog(dG, max.dist=80, direction=3*pi/4)

lines(v0)
lines(v45, col=2)
lines(v90, col=3)
lines(v135,col=4)

v4 <- variog4(dG, max.dist=100)
plot(v4)
plot(v4, omn=T)


ml.iso <- likfit(dG, ini=c(500, 20))
ml.aniso <- likfit(dG, ini=c(500, 20), fix.psiR=F, fix.psiA=F)
ml.anisoR <- likfit(dG, ini=c(500, 20), fix.psiR=F, fix.psiA=T, psiA=0)

logLik(ml.iso)
logLik(ml.anisoR)
logLik(ml.aniso)


## Inf. Bayesiana
summary(dG)

args(krige.bayes)
args(model.control)
MC <- model.control()
args(prior.control)
PC <- prior.control(phi.discrete=seq(0, 40, len=41))
args(output.control)
OC <- output.control(n.post=1000)

kb <- krige.bayes(dG, model=MC, prior=PC, output=OC)
names(kb)
plot(kb)

## outra priori
PC <- prior.control(phi.discrete=seq(0, 40, len=41), 
                    phi.prior="rec")
kb <- krige.bayes(dG, model=MC, prior=PC, output=OC)
plot(kb)


## priori tb tem tau^2_Rel
PC <- prior.control(phi.discrete=seq(0, 50, len=26), 
                    phi.prior="rec",
                    tausq.rel.discrete=seq(0, 1.5, length=16),
                    tausq.rel.prior="rec")
kb <- krige.bayes(dG, model=MC, prior=PC, output=OC)
par(mfrow=c(1,2))
plot(kb)

##kb <- krige.bayes(dG, loc=gr, model=MC, prior=PC, output=OC)
##save(kb, file = "kb.RData")
download.file("http://www.leg.ufpr.br/~paulojus/dadosgeo/kb.RData", destfile="kb.RData")
load("kb.RData")

names(kb)
names(kb$post)
names(kb$pred)

image(kb)

p65bayes <- apply(kb$pred$simul, 1, function(x) mean(x>65))
image(kb, values = p65bayes)

## comparando as predicoes
image(kc.ml)
image(kb)

plot(kc.ml, kb$predictive$pred)