#getwd() #setwd("/home/paulojus... ## 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) ## agora considerando tendências plot(s100, low=T, trend="1st") plot(s100, low=T, trend="2nd") plot(s100, low=T, trend=~coords[,1]) ## tendência? dependência espacial? ou ambos? ## análise descritiva de padrão espacial - variograma v <- variog(s100, max.dist=1) plot(v) v.env <- variog.mc.env(s100, obj.v=v) plot(v, env=v.env) args(variog) plot(variog(s100, max.dist=1, opt="cloud")) plot(variog(s100, max.dist=0.6)) plot(variog(s100, max.dist=1, trend="1st")) plot(variog(s100, max.dist=1, trend="1st")) ## 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))