## 1. Baixar os arquivos da página do curso em um diretório (pasta) de trabalho (workspace) ## apontar o R para este "workspace" #getwd() #setwd(file.choose()) # visializar os dados #file.show("dadosIntro.txt") dGeo <- read.table("dadosIntro.txt", header=T) ## se necessário apotar para outro diretório pode-se usar o mecanismo para procura do arquivo #dGeo <- read.table(file.choose(), header=T) ## Informação básica sobre os dados head(dGeo) dim(dGeo) str(dGeo) summary(dGeo) ## visualizando as coordenadas dos dados plot(dGeo[,1:2], asp=1) # ver janela gráfica ## definindo uma malha de pontos para predição gridGeo <- expand.grid(coords.x=0:100, coords.y=0:100) points(gridGeo, pch=19, col=2, cex=0.25) # ver janela gráfica ## Alguns métodos de interpolação ## Vamos considerar: ## ## 1. IDW Innverse distance weighting require(gstat) ## transformando os dados para o formato do pacotes "sp" ## (conveniente para trabalhar com funcoes do pacote gstat) ## vamos duplicar os dados (para preservar o original) ... dSP <- dGeo class(dSP) ## ... e converter para classe sp coordinates(dSP) <- ~coords.x + coords.y #coordinates(dSP) <- c("coords.x", "coords.y") class(dSP) plot(dSP) spplot(dSP, "dados") summary(dSP) str(dSP) ## da mesma forma vamos duplicar e converter a malha de predição gridSP <- gridGeo class(gridSP) coordinates(gridSP) <- ~coords.x + coords.y class(gridSP) gridded(gridSP) <- TRUE class(gridSP) summary(gridSP) ## IDW predIDW <- idw(dados ~1, dSP, gridSP) spplot(predIDW) ## O idw pode usar diferentes potências das distâncias predIDW <- idw(dados ~ 1, dSP, gridSP, idp=1) spplot(predIDW) class(predIDW) str(predIDW) ## armezanando no objeto da malha gridSP <- SpatialPixelsDataFrame(gridSP, predIDW@data) ## interpolando diferentes potências da distância gridSP$idp0.5 <- idw(dados ~ 1, dSP, gridSP, idp = 0.5)$var1.pred gridSP$idp02 <- idw(dados ~ 1, dSP, gridSP, idp = 2)$var1.pred gridSP$idp05 <- idw(dados ~ 1, dSP, gridSP, idp = 5)$var1.pred gridSP$idp10 <- idw(dados ~ 1, dSP, gridSP, idp = 10)$var1.pred ## comparando as interpolacoes names(gridSP) spplot(gridSP, c("idp0.5", "idp02", "idp05", "idp10"), main = "IDW") #################################################################### # "Local trend"Interpolacao por polinomios locais gridLT = krige(dados ~ coords.x+coords.y, dSP, gridSP, nmax = 10) # visualizacao nao fica bo pois mistura escalas da predicao e veriancia spplot(gridLT) ## melhor visualizar cada uma individualmente spplot(gridLT["var1.pred"]) spplot(gridLT["var1.var"]) gridLT3 = krige(dados ~ 1, dSP, gridSP, nmax = 10, degree=3) spplot(gridLT3) spplot(gridLT3["var1.pred"]) spplot(gridLT3["var1.var"]) ######################################################## # 3. Interpolação por GAM (generalised additiva models) require(mgcv) mGAM <- gam(dados ~ s(coords.x, coords.y), data=dGeo) mGAM gridSP$GAM <- predict(mGAM, newdata=gridGeo) spplot(gridSP, c("GAM"), main = "GAM") # poderiam ter sido obtidas as veriancias (erros pardao) de predicao) ################################################# # 4. Interpolação geoestatística - krigagem ## primeiro usa-se os dados para inferir estrutura de dependencia dVario <- variogram(dados ~ 1, data=dSP) plot(dVario) ## ajusta-se visualmente algum modelo para o variograma dVM <- vgm(400, "Sph", range=20, nugget=100) plot(dVario) lines(variogramLine(dVM, maxdist=100)) # o ajuste pode ser refinado por algum metosdo dVariofit = fit.variogram(dVario, vgm(400, "Sph", range=20, nug=100)) lines(variogramLine(dVariofit, 100)) ## e com isto pode-se obter a interpolacao por krigagem OK <- krige(dados~1, dSP, gridSP, model = dVariofit) spplot(OK["var1.pred"], main = "Krigagem Ordinária") spplot(OK["var1.var"], main = "Variancias de Krigagem") spplot(OK, main = "Krigagem Ordinária")