## ## Análise geoestatística utilizando o pacote geoR do programa R ## Exemplo introdutório ## ## carregando o pacote require(geoR) ## 1. OBTENDO VISUALIZADO/RESUMINDO DADOS ## carregando um conjunto de dados já disponível no pacote: data(s100) ## obtendo informações sobre o dado já disponível no pacote ## verificando o "tipo" de objeto ## neste caso o objeto já é da classe "geodata" class(s100) ## verificando a estrutura do objeto str(s100) ## obtendo um resumo dos dados: summary(s100) ## visualizando os dados (da classe geodata) points(s100) ## ver gráfico ## 2. ESTIMAÇÃO DE PARÃMETROS # 2.1 ESTIMAÇÃO POR MÁXIMA VEROSSIMILHANÇA ## estimando os parâmetros do modelo pelo métodos da máxima verossimilhança ## a estimação exige valores ("chutes") iniciais da variância e alcance e vamos utilizar os seguintes: INI <- c(var(s100$data), max(dist(s100$coords))/10) fitML <- likfit(s100, ini=INI) ## verificando o resultado do ajuste (estimativas dos parâmetros) fitML ## resultados mais detalhados do ajuste summary(fitML) ## 3. INTERPOLAÇÃO ESPACIAL ## definindo coordenada de um ponto para predição LOC <- t(c(0.5, 0.5)) ## visualizando no mapa points(s100) points(LOC, pch=4, cex=1.5, col=2) ## ver gráfico ## predição (interpolação) por krigagem pred1 <- krige.conv(s100, loc=LOC, krige=krige.control(obj.model=fitML)) ## vendo o valor predito e a variância de predição pred1[1:2] ## definindo a predição em outros pontos LOC <- cbind(c(0.5, 0.1, 0.9, 0.6, 0.6, 0.05), c(0.5, 0.9, 0.9, 0.7, 0.3, 0.35)) LOC ## visualizando no mapa points(s100) text(LOC[,1], LOC[,2], 1:6, cex=1.5, col=2) ## ver gráfico ## predição (interpolação) por krigagem pred1 <- krige.conv(s100, loc=LOC, krige=krige.control(obj.model=fitML)) ## vendo o valor predito e a variância de predição ## analisar e discutir os resultados em função da localização dos pontos a serem preditos pred1[1:2] ## definindo agora uma malha de pontos para predição gr <- expand.grid(seq(0,1, length=50), seq(0,1, length=50)) ## visualizando os pontos de predição no mapa points(s100) points(gr, cex=0.3, col=2) pred.gr <- krige.conv(s100, loc=gr, krige=krige.control(obj.model=fitML)) ## mapa de valores preditos image(pred.gr) ## ver gráfico ## sobrepondo os dados ao mapa points(s100, add=T) ## mudando a palheta de cores image(pred.gr, col=terrain.colors(20)) ## ver gráfico points(s100, add=T) ## mapa das variâncias de krigagem image(pred.gr, val=krige.var, col=terrain.colors(20)) ## ver gráfico points(s100, add=T) ## mapa dos erros padrão de krigagem image(pred.gr, val=sqrt(krige.var), col=terrain.colors(20)) ## ver gráfico points(s100, add=T) ## NOTE QUE: ## - AS PREDIÇÕES DEPENDEM DOS VALORES VIZINHOS ## - AS VARIÃNCIAS DE PREDIÇÃO DEPENDEM APENAS DA CONFIGURAÇÃO DOS VIZINHOS (E NÃO DOS VALORES!!) ## - ISTO É UMA CARACTERÍSTICA DO MODELO GAUSSIANO ## 2.2 ESTIMAÇÃO POR VARIOGRAMAS ## Variogramas possibilitam: ## - fazer uma análise exploratória/descritiva da estrutura de dependência dos dados ## - estimar os parâmetros por um método alternativo ao de máxima verossimillhança (** com ressalvas...) ## inicialmente vamos obtr a "nuvem variográfica" vc <- variog(s100, option="cloud") plot(vc) ## ver gráfico ## e, alternativamente, os pontos podem ser agrupados por classes para uma vesialização melhor v <- variog(s100) plot(v) ## ver gráfico ## note entretanto que não é recomendado utilizar todas as distâncias estre pontos vc <- variog(s100, option="cloud", max.dist=0.6) plot(vc) ## ver gráfico ## v <- variog(s100, max.dist=0.6) plot(v) ## ver gráfico ## mas neste caso ainda não verificamos a estabilizaçõa do variograma, o que pode ser devido a um alcance longo. ## vamos então permitir uma distância um pouco maior, mas que forneça um número razoável de pontos vc <- variog(s100, option="cloud", max.dist=1) plot(vc) ## ver gráfico ## v <- variog(s100, max.dist=1) plot(v) ## ver gráfico ## vendo resultados: v[1:3] ## neste caso ainda nao estabilizou claramente mas vamos manter este ultimo ## estimando parâmetros pelo variograma fitVG <- variofit(v) ## vendo os valores ajustados (compare cpm os de max. verossimilhança obtidos anteriormente) fitVG # adicionando a fnção de variograma ajustada ao gráfico lines(fitVG) ## (experimentar com diferentes opções para esta função) ## refazendo a krigagem com o ajuste de variograma predVG.gr <- krige.conv(s100, loc=gr, krige=krige.control(obj.model=fitVG)) ## mapa de valores preditos image(predVG.gr, col=terrain.colors(20)) ## ver gráfico ## sobrepondo os dados ao mapa points(s100, add=T) ## colocando as predições obtidas pelos dos métodos de estimação lado a lado para comparação par(mfrow=c(1,2)) image(pred.gr, col=terrain.colors(20)) image(predVG.gr, col=terrain.colors(20)) ## ver gráfico ## melhorando a comparação forçando uma escala comum para os dois gráficos ## e tb comparando os erros padrão de predição par(mfrow=c(2,3)) ZL <- range(c(pred.gr$pred, predVG.gr$pred)) image(pred.gr, col=terrain.colors(20), zlim=ZL) image(predVG.gr, col=terrain.colors(20), zlim=ZL) ## ver gráfico plot(pred.gr$pred, predVG.gr$pred, asp=1, xlim=ZL, ylim=ZL) ; abline(0,1) ZL1 <- range(sqrt(c(pred.gr$krige.var, predVG.gr$krige.var))) image(pred.gr, val=sqrt(krige.var), col=terrain.colors(20), zlim=ZL1) image(predVG.gr, val=sqrt(krige.var), col=terrain.colors(20), zlim=ZL1) ## ver gráfico plot(sqrt(pred.gr$krige.var), sqrt(predVG.gr$krige.var), asp=1, xlim=ZL1, ylim=ZL1) ; abline(0,1)