## Instalando o pacote geoR do programa R
install.packages("geoR", dep=T)
## carregando o pacote
require(geoR)
## visualizando um conjunto de dados disponível como exemplo do pacote
## Análises descritivas visuais
points(s100)
plot(s100)
plot(s100, low=T)
## análise descritiva de padrão espacial - variograma
plot(variog(s100, max.dist=1))
args(variog)
plot(variog(s100, max.dist=1, opt="cloud"))
plot(variog(s100, max.dist=1))
## estimando os parâmetros de um modelo geoestatístico simples
ml <- likfit(s100, ini=c(1, 0.2))
ml
## preparando um grid e fazendo um mapa de predição na área
args(pred_grid)
points(s100)
gr <- expand.grid(seq(0,1, leng=100), seq(0,1, leng=100))
points(gr, cex=0.3, col=2)
## algoritmo de interpolação (krigagem)
args(krige.conv)
kc <- krige.conv(s100, loc=gr, krige=krige.control(obj=ml))
## visualizando a predição
image(kc, col=terrain.colors(51))
points(s100, add=T)
## visualizando a incerteza das predições (erros padrão de krigagem)
image(kc, val=sqrt(kc$krige.var), col=terrain.colors(51))
points(s100, add=T)
## outras palhetas de cores e fatiamentos da imagem de predição
image(kc, col=gray(seq(1, 0, length=21)))
image(kc, col=gray(seq(1, 0, length=5)))
image(kc, col=gray(seq(1, 0, length=3)))
image(kc, col=gray(seq(1, 0, length=21)))
#################################################################
## Um outro conjunto de dados
points(parana)
str(s100)
str(parana)
plot(parana)
plot(parana, low=T)
## notar o efeito de uma "tendência" clara dos dados
plot(parana, tren="1st", low=T)
plot(variog(parana), max.dist=450)
plot(variog(parana, trend="1st"), max.dist=450)
## ajustando e comparando 2 modelosd: um sem e outro com tendência
ml.pr0 <- likfit(parana, ini=c(6000, 100))
ml.pr1 <- likfit(parana, ini=c(4000, 100), trend="1st")
ml.pr0
ml.pr1
AIC(ml.pr0)
AIC(ml.pr1)
## Interpolandop com o modelo com tendência
gr.pr <- pred_grid(parana$borders, by=10)
points(parana)
points(gr.pr, cex=0.25, col=2)
kc.pr <- krige.conv(parana, loc=gr.pr, krige=krige.control(obj=ml.pr1, trend.d="1st", trend.l="1st"),
borders=parana$borders)
image(kr.pr)
image(kc.pr, col=terrain.colors(51))