#===APRESENTAÇÃO DE PROGRAMA GEOESTATÍSTICO===# # #==Pacote geoestatístico utilizado: gstat==# # #==Aluno: Tito Nunes de Castro==# ## importando dados dado <- read.table("dadosintro.txt", head=T) ## instalando os pacotes necessaŕios para a análise estatística e suas dependências require(maptools) require(sp) require(gstat) #convertendo o data frame em um SpatialPointsDataFrame -> pacote sp dados = na.omit(dado) #eliminando os NAs do dataframe, pois o pacote gstat não trabalha com NAs coordinates(dados) = ~coords.x+coords.y dados #Visualizando os dados spplot(dados) bubble(dados, zcol='dados',fill = FALSE, col = "black", do.sqrt=FALSE, maxsize=2, identify = FALSE) #equivalente ao points() no pacote geoR #criando o grid de predição para a interpolação dos dados #1º passo: definir o grid x.range = as.integer(range(dados@coords[,1])) y.range = as.integer(range(dados@coords[,2])) #2º passo: criar o grid gr = expand.grid(x=seq(0,100, len=51), y=seq(0,100, len=51)) #3º passo: converter para a classe SpatialPixel coordinates(gr) = ~x+y gridded(gr) = TRUE #4ºpasso: visualização do grid em um gráfico plot(gr, col = "red") points(dados, pch=20, col="blue", cex=0.9) #verificação de anisotropia vari = variogram(g, alpha=c(0 ,90 ,180 ,360)) fitvar <- fit.variogram(vari, model=vgm(psill=300,"Exp",range=4,nugget=200,kappa=1.5)) plot(vari, model = fitvar, as.table=TRUE) #o ajuste se "encaixa" em todas as direções #Ajustando o variograma #1ºpasso: transformando o conjunto de dados em um objeto gstat #2ºpasso: plotar o variograma plot(variogram(dados ~ 1, dados)) #2º passo: criar o modelo do variograma mod = vgm(psill=300, "Exp", range=10, nugget=200, kappa=1.5) #3º passo: ajustar o modelo fitmodel = fit.variogram.reml(dados~1, dados, model=mod) # ajuste pela máxima verossimilhança restrita plot(variogram(dados~1,dados) ,fitmodel) #Validação cruzada da krigagem cross = krige.cv(dados~1, dados, model=fitmodel) #Plotando os valores observado pelo predito plot(dados$dados, cross$var1.pred, asp=1, xlab="Observado", ylab="Predito") abline(0,1,col="red",cex=0.5) #criação do mapa de krigagem map<-krige(dados~1,dados,model=fitmodel,newdata=gr) image(map) spplot(map,"var1.pred",col.regions=terrain.colors(50)) #==COMPARAÇÃO COM O geoR==# require(geoR) dG <- as.geodata(dado, coords.col=c(1,2), data.col=3) lf <- likfit(gdata, ini.cov.pars=c(300,10), nugget=200, cov.model="exponential") grid <- expand.grid(seq(0,100, len=51), seq(0,100, len=51)) krig <- krige.conv(dG, loc=grid, krige=krige.control(obj.m=lf)) image(krig); title('geoR') X11() image(map); title('gstat')