#-------------------------------------------------------------------------------------------------- # Programa para Análise Geoestatística canônica comparativa para os pacotes geoR, gstat e geospt # # Carvalho, Samuel de Padua Chaves e # Departamento de Ciencias Florestais # Universidade de São Paulo - ESALQ # Piracicaba, Brasil # spcarvalho@usp.br # # Revisado em: 18-10-2012 #-------------------------------------------------------------------------------------------------- ## Instalaçao dos pacotes necessários ## Necessita versao >= 2.15.0 para o geospt #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("maptools", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) #install.packages("sp", repos="http://brieger.esalq.usp.br/CRAN", dep=TRUE) ## Carregando os pacotes require(geospt) require(gstat) require(maptools) ## funções para importação/exportação e manipulação de mapas e dados geográficos gpclibPermit() ## função para habilitar licença de uso. A resposta deve ser TRUE require(sp) ## representação de dados espaciais no R require(geoR) require(lattice) ## Importando os dados setwd("E:/Doc_USP/LCE 5700 - Geoestatistica/Trabalho/apresentacao_programa") 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") ## transforma os dados para a classe do tipo SpatialPointsDataFrame ## Analise exploratoria dos dados ## Analises sugeridas no livro "Applied Spatial Data Analysis with R" pagina 192 X11() bubble(dados_gstat, "dados", key.space = "bottom") ## funcao do pacote sp X11() points(dG) ## Estimando a correlacao espacial: Variograma ## geoR v <- variog(dG, max.dist=summary(dG)[[3]][2]) X11() plot(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))) ## Mostra 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) ## Sill=533,1; Range=18,19; Nugget=16,16 vf_geor <- variofit(v, ini=ef) vf_geor; ## gstat vgm(ef[[1]][[2]][1],'Exp',ef[[1]][[2]][2], ef[[1]][3][[1]]) vf_gstat <- fit.variogram.reml(dados~1, ~coords.x+coords.y, dados, model=vgm(ef[[1]][[2]][1],'Exp',ef[[1]][[2]][2], ef[[1]][3][[1]])) vf_gstat ## Usando a likfit da geoR ml <- likfit(dG, ini=c(ef[[1]][[2]][1], ef[[1]][[2]][2])) ml #mesmos resultados do anterior porem diferente do variofit. Beta é o parametro para media? Por que este valor é diferente de ## Resultados diferentes entre variofit e likfit mean(dados$dados);ml$beta; ## Predicao Espacial ## Definindo o grid de predicao gr <- expand.grid(seq(0,summary(dG)[[3]][2], by=0.5), seq(0,summary(dG)[[3]][2], by=2)) names(gr) <- c("coords.x","coords.y") ## geoR krig.geor <- krige.conv(dG, loc=gr, krige=krige.control(obj.m=ml)) ## 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() spplot(krig.gstat["var1.pred"], main = "ordinary kriging predictions\npackage-gstat", col.regions=terrain.colors(21), scale=list(draw=T)) X11() spplot(krig.geor$predict, main = "ordinary kriging predictions\npackage-geoR", col.regions=terrain.colors(21), scale=list(draw=T)) ## nao rodou ## Usando o pacote geospt summary(dados_gstat) names(dados_gstat) m <- vgm(ef[[1]][[2]][1],'Exp',ef[[1]][[2]][2], ef[[1]][3][[1]]) # leave-one-out cross validation: out <- krige.cv(dados~1, dados_gstat, m, nmax = 100) # Saidas: MPE, ASEPE, RMSPE, MSPE, RMSSPE, R2 c.cv <- criterio.cv(out) #cross validation c.cv args(rbf) pred <- rbf(eta=c.cv$RMSPE, z=dados$dados, coordinates=dados[,1:2], newdata=gr, n.neigh=10, func="EXPON") #pred é um arquivo do tipo SpatialPixelsDataFrame names(pred) coordinates(pred) = c("x", "y") gridded(pred) <- TRUE # show prediction map X11() spplot(pred["var1.pred"], col.regions=terrain.colors(21),main = "Predictions package-geospt", key.space=list(space="right", cex=0.8)) X11() image(pred["var1.pred"],col=terrain.colors(21)) ## Graficos de preditos plot((krig.gstat)[[1]],krig.geor$predict) plot((krig.gstat)[[1]],pred[[1]]) plot(pred[[1]], krig.geor$predict) ## 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. ## Curiosidade: Analises de estacionaridade: Nao ok para o dataframe dados !!! data(coalash) # Probabilidade de estacionaridade direcao norte-sul pocket.plot(coalash, "PPR", coalash$x, coalash$y, coalash$coalash, FALSE) #pocket.plot(dados, "PPR", as.integer(dados$coords.x), as.integer(dados$coords.y), as.numeric(dados$dados), FALSE) #rodou com erro # Analise da estacionaridade local da variancia na direcao norte-sul pocket.plot(coalash, "PVR", coalash$x, coalash$y, coalash$coalash, FALSE) #pocket.plot(dados, "PVR", as.integer(dados$coords.x), as.integer(dados$coords.y), as.numeric(dados$dados), FALSE) #rodou com erro # Probabilidade de estacionaridade direcao leste-oeste pocket.plot(coalash, "PPC", coalash$x, coalash$y, coalash$coalash, FALSE) #pocket.plot(dados, "PPC", as.integer(dados$coords.x), as.integer(dados$coords.y), as.numeric(dados$dados), FALSE) #nao rodou # Analise da estacionaridade local da variancia na direcao leste-oeste pocket.plot(coalash, "PVC", coalash$x, coalash$y, coalash$coalash, FALSE) #pocket.plot(dados, "PVC", as.integer(dados$coords.x), as.integer(dados$coords.y), as.numeric(dados$dados), FALSE) #nao rodou