##### AULA 10 - GEOSTATISTICA

### AULA 10.4 - PREDICAO

#getwd()
#rm(list=ls())

### INSTAL E CARREG DOS PCTES E IMPORT DOS ARQUIVOS

install.packages("gstat") # pacote q faz 'idw'

require(gstat)
require(geoR)
require(maptools)     ## funCOes para importacao/exportacao e manipulacao 
#gpclibPermit()        ## funcao para habilitar licenca de uso
require(sf)           ## (classes) para representacao de dados espaciais 
require(spdep)

# 1. colocar arq 'deng.nv' e 'den.likfit' na pasta da workspace

source("./resol_aulas_12_1a4/deng.nv.R")
plot(deng.nv,low=T)

source("./resol_aulas_12_1a4/den.likfit.R")
den.likfit

### 2. DEFINICAOO DE UM GRID DE PREDICAO

summary(deng.nv) #para ver os minimos e max das coord para uso para fazer o grid
# pegar os min e max da borda

# $borders.summary
#        [,1]    [,2]
# min 662788.3 7691802
# max 673269.6 7705258

# criacao de um grid de predicao
# com base no retangulo envolvente acima, contendo 100 col e 100 linhas

gr=expand.grid(seq(662788,673270,len=100),seq(7691802,7705258,len=100))

# funcao 'points'
# 'pch' - simbolo a ser usado para plotagem - 'pch=19' - circulo solido

# 'cex' - quantidade pela qual os caracteres de plotagem e os simbolos serao escalonados em relacao ao valor default (cex=1.0)
# no caso do comando usado acima, ele define a dimensao de cada ponto do grid
# 'col' - codigo da cor

points(deng.nv)
points(gr,pch=19,cex=0.25,col=5)
head(gr)

### 3. CRIACAO DE UM DATA FRAME COM OS DADOS A SEREM UTILIZADOS

##criando dada.frame com as incid, as 3 covariaveis e as coordenadas

CoordX=deng.nv$coords[,1]
head(CoordX)
CoordY=deng.nv$coords[,2]
head(CoordY)
dados.inc=deng.nv$data
head(dados.inc)
head(deng.nv$covariates)
dados.ch.3inst=deng.nv$covariates[,7]
head(dados.ch.3inst)
dados.ch.2sm=deng.nv$covariates[,12]
head(dados.ch.2sm)
dados.prop.cas=deng.nv$covariates[,24]
head(dados.prop.cas)

dados.tdos=list()##dizendo que eh uma lista

dados.tdos=data.frame(coordx=CoordX, coordy=CoordY,inc=dados.inc,ch3inst=dados.ch.3inst,propcas=dados.prop.cas,ch2sm=dados.ch.2sm)##criando o data.frame

names(dados.tdos)##para ver os nomes das colunas
head(dados.tdos)
as.data.frame(dados.tdos)##comando para transform lista em data.frame; so para ter certeza q eh um d.f
class(dados.tdos)##para ver a classe

## para deixar os nomes da coord do grid igual ao do "dados.tdos"
head(gr)
grnovo=list()
grnovo=data.frame(coordx=gr$Var1,coordy=gr$Var2)
head(grnovo)
class(grnovo)

### 4. PREDICAO DOS VALORES DAS VAR DE TEND PARA OS PONTOS DO GRID

# variavies: CH_3INST, PROP_CASAS e CH_ATE2SM
# predicao nao estocastica, isto eh, sem usar geoestatistica
# usar um algoritmo, por exemplo, idw (inverso da distancia)
# os valores dessas tres variaveis devem ser criados para cada ponto do grid
# fazer uma variavel de cada vez

# predicaoo de 'ch3inst'
IDWch3inst <- idw(ch3inst~1,locations=~coordx+coordy,data=dados.tdos, newdata=grnovo)
head(IDWch3inst)

#predicaoo dos valores da var "propcas"
IDWpropcas <- idw(propcas~1,locations=~coordx+coordy,data=dados.tdos, newdata=grnovo)
head(IDWpropcas)

#predicaoo dos valores da var "ch2sm"
IDWch2sm <- idw(ch2sm~1,locations=~coordx+coordy,data=dados.tdos, newdata=grnovo)
head(IDWch2sm)

# juntando os tres bancos em um unico data frame
# procedimento interessante para verificar se os valores foram preditos nos mesmos ptos XY

##criacao do data.frame com os valores preditos das covariaveis

cov.pred=list()
cov.pred=data.frame(IDWch3inst,IDWpropcas, IDWch2sm)
class(cov.pred)
head(cov.pred)

##refazendo com os mesmos nomes do deng.nv e grnovo

head(grnovo)
str(deng.nv)
head(cov.pred)
cov.pred2=list()
cov.pred2=data.frame(coordx=cov.pred[,1],coordy=cov.pred[,2],CH_3INST=cov.pred[,3],PROP_CASAS=cov.pred[,7], CH_ATE2SM=cov.pred[,11])

# tirando as coord XY de cov.pred2; foram mantidas inicialmente para garantir q os dados estejam alinhados com as coords

head(cov.pred2)
cov.pred3=cov.pred2[,-c(1,2)]
head(cov.pred3)
#cov.pred4=cov.pred3[,-1]
#head(cov.pred4)

# agora tenho o arq com os valores preditos das tres variaeis para todo o grid - vou usar isso para predizer os valores de inc em td o grid

## 5. PREDICAO DOS VALORES DA INCIDENCIA EM TODOS OS PTOS DO GRID

# usar o bd 'deng.nv'
# usar o grid 'grnovo'
# usar os valores do variograma ajustados por 'likfit' (obj 'den.likfit')
# usar as informacao sobre a tendencia
# usar 'lam=0.5' - trab com a raiz da incid (mas os val pred serao de incid e nao da sua raiz)

kc=krige.conv(deng.nv,loc=grnovo,krige=krige.control(type.krige="OK",
              obj.m=den.likfit,trend.d=~CH_3INST+PROP_CASAS+CH_ATE2SM,
              trend.l=~as.matrix(cov.pred3),lam=0.5))

##mapa da predicao por krigagem

##x.leg(a,b): a = coord x de inicio da posicion da leg e b = coord de fim
##y.leg(a,b): a = coord y de inicio e b = coord de fim

image(kc,col=terrain.colors(21))
points(deng.nv,add=T)##adicionando os pontos

## mapa da variancia da krigagem

image(kc, col=terrain.colors(21), val=kc$krige.var, x.leg=c(670000,675000), y.leg=c(7692000,7692500))
points(deng.nv, add=T)

## mapa de erro padrao de predicao/krigagem (E[S|y])

image(kc, col=terrain.colors(21), val=sqrt(kc$krige.var), x.leg=c(670000,675000),	y.leg=c(7692000,7692500))
points(deng.nv, add=T)

## dividindo o mapa de krigagem em algums poucas (4) classes 

names(kc)#para ver o nomes dos objetod dentro de kc
summary(kc$predict)##para ver os valores minimos e maximo dos valores preditos
image(kc, breaks=c(435.5, 1703.0, 2752.0,3633.0,7376.0),x.leg=c(670000,675000),	y.leg=c(7692000,7692500),col=c("blue","green","yellow","red"))

summary(kc$krige.var) # para ver os val mim e max da variancia