#### AUAL 12 - GEOESTATiSTICA

### AULA 12.3 - VARIOGRAMA

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

### INSTAL E CARREG DOS PCTES

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

# 1. IMPORTA??O DO 'geodata' 'deng.nv

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

### 2. MODELAGEM DO VARIOGRAMA (considerando a tendencia)

summary(deng.nv)# para ver distancia maxima entre os pontos = dmax = 12730m

# fazer o variograma para metade da distancia - d = 6400m
# distancias grandes - poucos pares de pontos - gde variabil do variograma

# funcao 'variog'

variog1 <- variog(deng.nv,max.dist=6400)
plot(variog1) # plotagem do variograma

# variograma considerando transf ('lam=0.5)

variog2 <- variog(deng.nv,max.dist=6400, lam=0.5)
plot(variog2) # plotagem do variograma

# variograma considerando a transformacao ('lam=0.5')
# e a tendencia ('trend=~CH_3INST+ PROP_CASAS+CH_ATE2SM')

variog3 <- variog(deng.nv,max.dist=6400,lam=0.5,trend=~CH_3INST+ PROP_CASAS+CH_ATE2SM)
plot(variog3)# plotagem do variograma

# envelopes do variograma (monte carlo - 99 permutacoes)

var.env <- variog.mc.env(deng.nv,obj.variog=variog3)
plot(variog3,env=var.env) # variograma com envelope

# ajuste do variograma experim a um variograma teorico
# fazer e ajuste e depois salvar

ef.variog <- eyefit(variog3) 

# usar funcao exponencial e valores em torno de sigma: 628.4; phi: 3659; nug=409.8; ou digitar

ef.variog # mostra os valores dos parametros escolhidos

## PARES DE PONTOS NO VARIOGRAMA

# opcao 'cloud' - variograma em nuvem

variog3.nuven <- variog(deng.nv,max.dist=6400,lam=0.5,trend=~CH_3INST+ PROP_CASAS+CH_ATE2SM,option="cloud")

### OBS.: sobre uso ou nao de aspas
# uso aspas para informar alguma caracteristica
# nao uso aspas qdo me refiro a objeto

plot(variog3.nuven) # posso avaliar a existencia de outliers, princ., em peq distancias

##isotropia e anisotropia

# funcao 'variog4' - faz variograma em quatro direcoes 
# medidas em relacao ao eixo y (0, 45, 90 e 135 graus)

variog3.4dir <- variog4(deng.nv,max.dist=6400,lam=0.5,trend=~CH_3INST+ PROP_CASAS+CH_ATE2SM)
plot(variog3.4dir)

## Observacao: se nao tiver um bom motivo para concluir que existe anisotropia
## melhor usar modelo isotropico
## Muito do que eh considerado anisotropia, na realidade eh apenas variacao aleatoria
## no geoR - angulo de anisotropia - em relacao ao eixo 'y'

### 3. SIGNIFICADO DOS PARAMETROS DO VARIOGRAMA

ef.variog # mostra os valores dos paremetros escolhidos

# modelo exponencial - kappa = 0.5
# outros tipos de modelo exigem a especificacao de kappa e phi (Matern, por exemplo)

ef.variog2 <- eyefit(variog3) 

# Sigamsq - sill (ou contribuicao)
# tausq - nugett ou efeito pepita
# alcance pratico - funcao de kappa e phi

# funcao 'praticalRange' informa alcance pratico em funcao do modelo, do kappa e do phi

practicalRange("exp",phi=3659,kappa=0.5)

### 4. ESTIMACAO DOS PARAMETROS DO VARIOGRAMA POR MAX VEROSS - FUNCAO 'likfit'

# Variograma - da uma ideia dos valores de phi, tausk, sigma2 e dos betas
# uso funcao 'likfit' para estimar os valores dos parametros por max verossimilhanca

ef.variog

den.likfit <- likfit(deng.nv,lam=0.5,ini=c(628.4,3659),nug=409.8,trend=~CH_3INST+ PROP_CASAS+CH_ATE2SM)
den.likfit # mostra os valores dos parametros estimados com base na max veross
summary(den.likfit)
ls(den.likfit)# componentes do objeto 'den.likfit'
den.likfit$beta # mostra o conteudo de 'beta'
den.likfit$model.components # mostra as componentes de tend, espacial e residuos

# valores estimados por max verossim
den.likfit

## likfit: estimated model parameters:
##      beta0       beta1       beta2       beta3       tausq     sigmasq         phi 
## "  56.4888" "   0.4661" "   0.2934" "   0.3006" " 402.2121" "1231.1728" "3658.9993" 
## Practical Range with cor=0.05 for asymptotic range: 10961.38
# likfit: maximised log-likelihood = -3642

##AIC = (-2)x(-3642) + 2x5 = 7274

### 5. EXPORTACAO DO OBJETO OBTIDO COM O 'likfit'

dump("den.likfit","./resol_aulas_12_1a4/den.likfit.R")