Analise Espacial - Processos continuos
Análise Geoespacial da Base de Dados
Este script realiza uma análise geoespacial exploratória de uma base de dados contendo medições realizadas em diferentes profundidades. Utilizaremos pacotes específicos para manipulação de dados geoespaciais e visualização. A base de dados foi importada do arquivo castro.txt e contém coordenadas geográficas, além de valores registrados para diferentes faixas de profundidade.
Pacotes Utilizados
## --------------------------------------------------------------
## Analysis of Geostatistical Data
## For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
## geoR version 1.9-4 (built on 2024-02-14) is now loaded
## --------------------------------------------------------------Importação da Base de Dados
solo <- read.table("C:/Users/Maria/Desktop/processos continuos/castro.txt", header = TRUE)
head(solo)## Lat Lon X0.10 X10.15 X15.20 X20.25 X25.30 X30.35 X35.40 Xcd
## 1 608039.6 7250137 2.06 2.55 1.83 1.54 1.12 0.84 0.80 -49.93059
## 2 608063.1 7250146 2.03 1.58 2.21 2.36 2.61 2.44 1.76 -49.93036
## 3 608086.6 7250154 1.67 1.38 1.61 1.79 2.13 1.46 1.36 -49.93013
## 4 608007.6 7250152 1.23 1.02 0.96 1.19 1.39 1.45 0.76 -49.93091
## 5 608031.1 7250161 0.57 0.75 1.22 1.69 2.12 2.16 2.34 -49.93068
## 6 608054.6 7250169 1.67 1.60 1.20 1.75 1.70 1.58 1.21 -49.93045
## Ycd
## 1 -24.85994
## 2 -24.85986
## 3 -24.85978
## 4 -24.85980
## 5 -24.85973
## 6 -24.85965Estrutura e Resumo dos Dados
Nesta seção, exploramos a estrutura do dataset para compreender as variáveis e realizar uma análise inicial.
## [1] 358 11## X0.10 X10.15 X15.20 X20.25
## Min. :0.560 Min. :0.540 Min. :0.360 Min. :0.500
## 1st Qu.:1.690 1st Qu.:1.500 1st Qu.:1.383 1st Qu.:1.320
## Median :1.950 Median :1.740 Median :1.605 Median :1.600
## Mean :2.091 Mean :1.822 Mean :1.710 Mean :1.713
## 3rd Qu.:2.360 3rd Qu.:2.055 3rd Qu.:1.938 3rd Qu.:2.020
## Max. :4.920 Max. :3.900 Max. :3.980 Max. :3.860
## X25.30 X30.35 X35.40
## Min. :0.490 Min. :0.460 Min. :0.400
## 1st Qu.:1.363 1st Qu.:1.462 1st Qu.:1.330
## Median :1.695 Median :1.885 Median :1.890
## Mean :1.786 Mean :1.910 Mean :1.924
## 3rd Qu.:2.090 3rd Qu.:2.280 3rd Qu.:2.388
## Max. :3.990 Max. :3.930 Max. :4.790Análise Descritiva: Profundidade 35-40 cm
Para a análise descritiva, focaremos na profundidade de 35-40 cm, convertendo a base para um objeto geoespacial e gerando uma visualização inicial.
## Number of data points: 358
##
## Coordinates summary
## Lat Lon
## min 607546.1 7250137
## max 608125.0 7250863
##
## Distance summary
## min max
## 24.98484 728.01244
##
## Data summary
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.400000 1.330000 1.890000 1.923827 2.387500 4.790000
Mapa de Pontos (Superior Esquerdo): Mostra a distribuição dos pontos de amostragem na área de estudo, com diferentes símbolos indicando grupos de medições. Os pontos estão bem distribuídos, cobrindo de forma uniforme a área.
Dispersão Vertical e Horizontal (Superior Direito e Inferior Esquerdo): Mostra como a variável 35-40 cm varia conforme a localização. Há variabilidade, sugerindo uma possível tendência espacial.
Histograma com Curva de Densidade (Inferior Direito): A curva de densidade sugere uma distribuição próxima do normal, mas com uma leve assimetria à direita e alguns possíveis outliers, indicando variabilidade nas medições.
Verificação dos Pressupostos
Nesta etapa, faremos a verificação dos pressupostos de normalidade e homogeneidade da variância, utilizando a transformação de Box-Cox para ajustar possíveis assimetrias e identificar a melhor transformação para a variável analisada. Em seguida, exploraremos a visualização espacial dos dados para verificar padrões e outliers.
O gráfico de log-verossimilhança sugere o valor ótimo de lambda para a transformação dos dados, próximo de zero. Isso indica que uma transformação logarítmica pode ser adequada para estabilizar a variância e aproximar os dados de uma distribuição normal.
Ajustando com Transformação Logarítmica
solo$log_35_40 <- log(solo[, 9] + 1)
geosolo <- as.geodata(solo, coords.col = c(1, 2), data.col = "log_35_40")
boxcox(geosolo)
Box-Cox Plot: Após a transformação, o gráfico sugere um lambda ótimo próximo de 1, indicando que os dados agora estão próximos da normalidade. O intervalo de confiança inclui 1, o que dispensa transformações adicionais. Podemos prosseguir com a análise utilizando os dados ajustados.
Gráfica Geosolo: A transformação estabilizou a variabilidade e melhorou a normalidade dos dados, permitindo prosseguir com as análises espaciais.
Visualização Espacial (Quintis): A distribuição espacial dos pontos mostra clusters nas extremidades, refletindo variabilidade nas medições. As bolhas maiores e vermelhas indicam valores mais altos, sugerindo áreas de maior concentração ou possíveis anomalias.
Ajuste do Variograma e Definição de Modelos
Nesta etapa, iremos ajustar visualmente o variograma experimental e testar diferentes distâncias máximas. O objetivo é identificar a distância de influência e ajustar modelos teóricos ao variograma, facilitando a estimativa de parâmetros para predição espacial.
Definindo as Bordas e Plotando o Variograma
bord <- read.csv("C:/Users/Maria/Desktop/ESPACIAL/borda.csv")
geosolo$borders <- with(bord, cbind(V1, V2))
points(geosolo, pt.div = "quintile", xlab = "leste", ylab = "norte")
## variog: computing omnidirectional variogram
Crescimento Inicial: A semivariança aumenta rapidamente para distâncias menores, indicando que há uma forte autocorrelação espacial entre pontos próximos.
Platô (Sill): A partir de aproximadamente 300 metros, o variograma começa a estabilizar, formando um platô. Isso indica que, além dessa distância, os pontos não estão mais correlacionados.
Alcance (Range): O alcance é a distância em que o variograma atinge o platô. Aqui, ele ocorre em torno de 300 metros, representando a extensão da dependência espacial.
Ajuste Visual do Variograma (Eyefit)
Nesta seção, ajustamos visualmente o modelo teórico ao variograma experimental usando a função eyefit().
## cov.model sigmasq phi tausq kappa kappa2 practicalRange
## 1 exponential 0.07 144.23 0.03 <NA> <NA> 432.07446580548O ajuste visual confirma a presença de dependência espacial até cerca de 350 metros. O modelo teórico se mostra adequado para predição espacial, validando os parâmetros estimados para a próxima etapa de krigagem.
Testando Outras Distâncias Máximas
Para avaliar a robustez do ajuste, testamos outras distâncias máximas (500 e 550).
## variog: computing omnidirectional variogram
## cov.model sigmasq phi tausq kappa kappa2 practicalRange
## 1 exponential 0.07 160.26 0.03 <NA> <NA> 480.096054149432## variog: computing omnidirectional variogram
## cov.model sigmasq phi tausq kappa kappa2 practicalRange
## 1 exponential 0.07 176.28 0.03 <NA> <NA> 528.087685170648Os ajustes para distâncias maiores mantêm consistência nos parâmetros, indicando que a dependência espacial é capturada até cerca de 350 metros, com patamar próximo de 0,08 e efeito pepita visível. Estes resultados são sólidos para prosseguir com a modelagem e krigagem, dado que os parâmetros permanecem estáveis mesmo com diferentes distâncias máximas analisadas.
Análise Direcional: Nuvens de Pontos do Variograma
Utilizamos variog4() para explorar a anisotropia, ou seja, a variação da dependência espacial em diferentes direções. Esta análise identifica padrões na autocorrelação espacial, fundamentais para a modelagem e a krigagem.
## variog: computing variogram for direction = 0 degrees (0 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
## tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
Diferenças Entre Direções: O gráfico mostra os variogramas nas direções 0º, 45º, 90º e 135º, revelando variações significativas e confirmando a presença de anisotropia.
Direções Preferenciais: A direção 90º (leste-oeste) apresenta maior semivariança, indicando menor correlação espacial. Já a direção 0º (norte-sul) possui menor semivariança, sugerindo maior continuidade.
Comprovada a anisotropia, a modelagem deve considerar essas variações para melhorar a precisão das predições.
Estimação de Parâmetros e Modelos de Krigagem
Nesta etapa, ajustaremos diferentes modelos de variograma utilizando a função likfit(). O objetivo é comparar esses modelos (circular, exponencial, cauchy e gaussiano) e identificar o que melhor se adequa aos dados, com base nos critérios de AIC e BIC. Em seguida, aplicaremos a krigagem para fazer predições espaciais.
Ajuste dos Modelos de Variograma
## kappa not used for the circular correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.## kappa not used for the gaussian correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
## arguments for the maximisation function.
## For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
## times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.## [1] 5.165721 3.664446 4.755153 8.229058## [1] 20.68785 19.18658 20.27728 23.75119O modelo exponencial apresentou os menores valores tanto para AIC quanto para BIC, indicando ser o modelo mais adequado para representar a dependência espacial dos dados.
Os modelos gaussiano e circular tiveram os piores ajustes, com maiores valores de AIC e BIC.
Visualizaçao do melhor modelo:
## Summary of the parameter estimation
## -----------------------------------
## Estimation method: maximum likelihood
##
## Parameters of the mean component (trend):
## beta
## 0.0734
##
## Parameters of the spatial component:
## correlation function: exponential
## (estimated) variance parameter sigmasq (partial sill) = 0.0426
## (estimated) cor. fct. parameter phi (range parameter) = 104.5
## anisotropy parameters:
## (fixed) anisotropy angle = 0 ( 0 degrees )
## (fixed) anisotropy ratio = 1
##
## Parameter of the error component:
## (estimated) nugget = 0.0395
##
## Transformation parameter:
## (fixed) Box-Cox parameter = 0 (log-transformation)
##
## Practical Range with cor=0.05 for asymptotic range: 313.0794
##
## Maximised Likelihood:
## log.L n.params AIC BIC
## "2.168" "4" "3.664" "19.19"
##
## non spatial model:
## log.L n.params AIC BIC
## "-46.78" "2" "97.55" "105.3"
##
## Call:
## likfit(geodata = geosolo, ini.cov.pars = c(0.05, 150), nugget = 0.05,
## lambda = 0, cov.model = "exponential")Visualização Gráfica dos Modelos
## variog: computing omnidirectional variogramlines.variomodel(ml1, col = "red") # Modelo Circular
lines.variomodel(ml2, col = "blue") # Modelo Exponencial
lines.variomodel(ml3, col = "green") # Modelo Cauchy
lines.variomodel(ml4, col = "black") # Modelo Gaussiano
legend("bottom", legend = c("Circular", "Exponencial", "Cauchy", "Gaussiano"),
col = c("red", "blue", "green", "black"), lty = 1, cex = 0.8, horiz = TRUE, bty = "n")
O gráfico mostra o ajuste dos modelos de variograma (circular, exponencial, cauchy e gaussiano) aos dados experimentais:
O modelo exponencial (azul) apresentou o melhor ajuste, seguindo de perto os pontos observados.
O modelo circular (vermelho) superestima a semivariança em maiores distâncias.
O modelo cauchy (verde) subestima a semivariança em curtas distâncias.
O modelo gaussiano (preto) oferece um ajuste intermediário, mas não acompanha bem os dados extremos.
O modelo exponencial é o mais adequado, apresentando o melhor equilíbrio entre ajuste e simplicidade.
Definindo o Grid de Predição
Criamos um grid de predição para aplicar a krigagem, cobrindo toda a área de estudo.
GR <- pred_grid(geosolo$borders, by = 10)
# Visualizando o grid de predição
points(geosolo)
points(GR, col = 2, pch = 19, cex = 0.3)
# Refinando o grid dentro das bordas
gr0 <- GR[.geoR_inout(GR, geosolo$borders),]
points(gr0, col = 4, pch = 19, cex = 0.3)
## [1] 2247 2Krigagem e Predição Espacial
Nesta etapa, aplicamos a krigagem para estimar os valores espaciais, utilizando os modelos ajustados.
Krigagem com o Modelo Circular
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodimage(geosolo.kr1, col = terrain.colors(15), x.leg = c(607450, 607880), y.leg = c(7250170, 7250200))
Krigagem com o Modelo Exponencial
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodimage(geosolo.kr2, col = terrain.colors(15), x.leg = c(607500, 607700), y.leg = c(7250150, 7250200))
Krigagem com o Modelo Cauchy
## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodimage(geosolo.kr3, col = terrain.colors(15), x.leg = c(607500, 607700), y.leg = c(7250150, 7250200))
Todos os modelos captaram bem a continuidade espacial, refletida pela transição suave das cores nos mapas. O modelo exponencial destacou-se como o mais adequado, oferecendo maior precisão nas estimativas e identificando padrões claros de variação espacial. Os modelos circular e cauchy também mostraram boas estimativas, embora menos detalhadas, indicando consistência nas tendências espaciais.
Predição Espacial com Krigagem
Nesta etapa, aplicaremos a krigagem utilizando o modelo mais adequado, o exponencial (ml2), para gerar um grid de predição e realizar a interpolação dos dados. Além disso, serão apresentados mapas de probabilidades condicionais e quantis para melhor interpretação dos resultados.
Krigagem com o Modelo Exponencial
Utilizaremos a krigagem ordinária com o controle definido para o modelo exponencial ajustado (ml2). A seguir, aplicamos a interpolação espacial para estimar os valores e criamos o mapa resultante.
# Configurações da krigagem com o modelo exponencial
KC <- krige.control(type = "OK", obj.model = ml2, lam = 0)
# Controle de saída com limiar para probabilidade condicional
OC <- output.control(thres = 1.890, quantile = c(0.1, 0.9))
# Aplicando a krigagem
geosolo.kc <- krige.conv(geosolo, locations = gr0, krige = KC, output = OC)## krige.conv: results will be returned only for prediction locations inside the borders
## krige.conv: model with constant mean
## krige.conv: performing the Box-Cox data transformation## Warning in krige.conv(geosolo, locations = gr0, krige = KC, output = OC):
## .Random.seed not initialised. Creating it with runif(1)## krige.conv: sampling from the predictive distribution (conditional simulations)
## krige.conv: back-transforming the simulated values
## krige.conv: back-transforming the predicted mean and variance
## krige.conv: Kriging performed using global neighbourhoodimage(geosolo.kc, loc = GR, bor = geosolo$borders,
col = terrain.colors(15), x.leg = c(607500, 607770), y.leg = c(7250150, 7250200))
Análise de Probabilidade Condicional
A seguir, avaliamos a probabilidade condicional para um limiar específico, gerando um mapa de probabilidades com base nos valores estimados.
image(geosolo.kc, loc = GR, bor = geosolo$borders,
col = terrain.colors(21), val = 1 - geosolo.kc$prob,
x.leg = c(4950, 5450), y.leg = c(4770, 4830))
O mapa apresentado mostra as probabilidades condicionais da krigagem para o modelo exponencial. O valor plotado é 1 - probabilidade, o que indica a chance de que a variável de interesse seja menor do que o limiar definido (no caso, 1.890).
As áreas verdes escuras representam regiões com alta probabilidade de que os valores estimados sejam menores do que o limiar, sugerindo maior certeza nas predições para estas áreas.
As regiões amarelas e rosadas indicam menor probabilidade, ou seja, maior incerteza e possibilidade de que os valores observados excedam o limiar definido.
O mapa reflete boa continuidade espacial, com transição suave entre áreas de maior e menor probabilidade, o que indica um ajuste eficiente da krigagem.
O mapa captura bem as áreas de maior confiabilidade nas estimativas. A predominância da cor verde sugere que a maior parte da área estudada apresenta alta probabilidade de valores abaixo do limiar, reforçando a adequação do modelo exponencial ajustado.
Mapa de Quantis
Os mapas de quantis permitem identificar as áreas com maior e menor confiabilidade nas estimativas. Apresentamos os mapas para o 10º e 90º percentis, representando “baixos confiáveis” e “altos confiáveis”, respectivamente.
image(geosolo.kc, loc = GR, bor = geosolo$borders,
col = terrain.colors(21), val = geosolo.kc$quantile$q90,
x.leg = c(4950, 5450), y.leg = c(4770, 4830),
main = "Mapa do 90º Percentil - Baixos Confiáveis")
image(geosolo.kc, loc = GR, bor = geosolo$borders,
col = terrain.colors(21), val = geosolo.kc$quantile$q10,
x.leg = c(4950, 5450), y.leg = c(4770, 4830),
main = "Mapa do 10º Percentil - Altos Confiáveis")
Mapa do 10º Percentil (“Altos Confiáveis”): Indica as áreas com maior probabilidade de altos valores, destacadas em amarelo e laranja. As regiões verdes mostram concentrações mais baixas.
Mapa do 90º Percentil (“Baixos Confiáveis”): Exibe as zonas com alta probabilidade de valores baixos, em verde. As áreas em amarelo indicam menor intensidade, mostrando uma transição suave entre regiões.
Conclusão
A análise geoespacial indicou forte dependência espacial, com o modelo exponencial apresentando o melhor ajuste. Os variogramas e a krigagem representaram bem a estrutura dos dados.
Os mapas de krigagem evidenciaram áreas com diferentes intensidades, e os percentis identificaram zonas de maior e menor risco. Esses resultados oferecem uma boa base para interpretação e validação dos padrões espaciais observados.