#-------------------------------------------------------------------------------------------------- # Programa para Análise Geoestatística canônica: Krigagem convencional krigagem e em blocos para os pacotes geoR e gstat # # Carvalho, Samuel de Padua Chaves e # Departamento de Ciencias Florestais # Universidade de São Paulo - ESALQ # Piracicaba, Brasil # spcarvalho@usp.br # # Revisado em: 2012-10-26 #-------------------------------------------------------------------------------------------------- ## Instalaçao dos pacotes necessários #install.packages("geoR", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) #install.packages("geospt", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) #install.packages("sp", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) #install.packages("lattice", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) ## Carregando os pacotes require(gstat) require(geoR) ## Importando os dados setwd("Direcionar a pasta de trabalho") dados <- read.table('dadosintro.txt', h=T) head(dados) ## Cria o objeto dG(geodata) para analise com a geoR dG <- as.geodata(dados, coords.col=c(1,2), data.col=3); ## funcao do pacote geoR dados_gstat <- dados coordinates(dados_gstat) <- c("coords.x","coords.y") ## Analise exploratoria dos dados X11() bubble(dados_gstat, "dados", key.space = "bottom") ## funcao do pacote sp plot(dG, low=T) ## Estimando a correlacao espacial: Variograma Empírico ## geoR v <- variog(dG, max.dist=summary(dG)[[3]][2]) X11() plot(v) lines(variofit(v)) ## gstat X11() plot(variogram(dados ~ 1, dados_gstat)) X11() ## em todas direções plot(variogram(dados ~ 1, dados_gstat, alpha = c(0, 45, 90, 135))) ## Mostrando os modelos de variograma ##show.vgms() # pacote gstat ##show.vgms(model = "Mat", kappa.range = c(0.1, 0.2, 0.5, 1, 2, 5, 10), max = 10) ## Modelando o variograma ## geoR X11() plot(v) ef <- eyefit(v) vf_geor <- variofit(v, ini=ef) vf_geor; ## gstat vgm(ef[[1]][[2]][1],'Exp',ef[[1]][[2]][2], ef[[1]][3][[1]]) dados_g <- gstat(id = "dados", formula=dados~1, locations =~coords.x+coords.y, data = dados) vf_gstat <- fit.variogram(variogram(dados_g), model=vgm(psill=ef[[1]][[2]][1],'Exp',range=ef[[1]][[2]][2], nugget=ef[[1]][3][[1]])) vf_gstat X11() plot(variogram(dados ~ 1, dados_gstat), vf_gstat) #resultado estranho ## Usando a likfit da geoR ml <- likfit(dG, ini=ef) ml #mesmos resultados do anterior porem diferente do variofit. Beta é o parametro para media? Por que este valor é diferente de mean(dados$dados);ml$beta; # Valores diferentes, indicando correlaçao espacial nos dados ## Predicao Espacial ## Definindo o grid de predicao gr=expand.grid(y=seq(1,100,1), x=seq(1,100,1)) names(gr) <- c("coords.x","coords.y") # necessario para usar no krige ## geoR krig.geor <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ml)) points(dG) points(gr, pch='+', cex=0.25, col=2) ## gstat gr1 <- gr gridded(gr1) <- ~coords.x+coords.y ## transforma em SpatialPixels krig.gstat <- krige(dados ~ 1, dados_gstat, gr1, vf_gstat) ## Visualizando os mapas X11() image(krig.geor, col=terrain.colors(21)) points(dG, add=T, pch=19) contour(krig.geor, add=T, nlev=21) X11() image(krig.gstat, col=terrain.colors(21)) points(dG, add=T, pch=19) contour(krig.gstat, add=T, nlev=21) ## Outra forma de visualizacao X11() class(krig.gstat["var1.pred"]) spplot(krig.gstat["var1.pred"], main = "ordinary kriging predictions\npackage-gstat", col.regions=terrain.colors(21), scale=list(draw=T)) ## Graficos de preditos X11() par(mfrow=c(1,3)) plot((krig.gstat)[[1]],krig.geor$predict, xlab='Preditos - gstat', ylab='Preditos - geoR'); abline(0,1, col='red') ## Krigagem em bloco ## geoR grid <- gr X11() plot(grid$coords.x,grid$coords.y, pch='+', cex=.25) axis(1,at=seq(0,100,by=5)) abline(v=seq(25,100,25));abline(h=seq(25,100,25)); abline(v=0);abline(h=0); ## Definindo o tamanho do bloco len_block_x <- 25 len_block_y <- 25 ## Calculando o numero de blocos regulares cabíveis na area grid$nblock=max(grid$coords.x)/len_block_x; grid$nblock=max(grid$coords.y)/len_block_y; ## Atribuindo os valores preditos em cada bloco grid$xlab=ifelse((grid$coords.x <= 25),'A', ifelse((grid$coords.x > 25 & grid$coords.x <= 50),'B', ifelse((grid$coords.x > 50 & grid$coords.x <= 75),'C','D'))) grid$ylab=ifelse((grid$coords.y <= 25),'A', ifelse((grid$coords.y > 25 & grid$coords.y <= 50),'B', ifelse((grid$coords.y > 50 & grid$coords.y <= 75),'C','D'))) grid$bloco=paste(grid$xlab,grid$ylab) grid$xlab<-NULL; grid$ylab<-NULL; grid$nblock<-NULL; grid$preditos <- krig.geor$predict ## Numero de observaçoes por bloco tapply(grid$coords.x, grid$bloco, length) ## Calculo da media por bloco media_pred <- aggregate(list(ybarra=grid$preditos), list(bloco=grid$bloco), mean) ## Associando ao dataframe grid <- merge(grid, media_pred) ## Krigagem em bloco pela gstat ## Janela regular de 25 por 25 krig.gstat.block <- krige(dados ~ 1, dados_gstat, gr1, vf_gstat, block=c(25,25)) plot((krig.gstat)[[1]],krig.geor$predict) plot((krig.gstat)[[1]],(krig.gstat.block)[[1]],) plot((krig.gstat.block)[[1]],grid$ybarra) ## Visualizando os mapas da krigagem convencional comparado a krigagem em bloco pela gstat X11() image(krig.gstat, col=terrain.colors(21)) points(dG, add=T, pch=19) contour(krig.gstat, add=T, nlev=21) X11() image(krig.gstat.block, col=terrain.colors(21)) points(dG, add=T, pch=19) contour(krig.gstat.block, add=T, nlev=21) ## End run #========================================================================================================================### ## Material consultado # Programas desenvolvidos durante a disciplina LCE 5700 - Geoestatística # BISVAND, R. S. et al. Applied Spatial Data Analysis with R. Springer. 2008. New York. 390p.