Introdução

O presente trabalho tem como objetivo principal a modelagem da distribuição espacial e a interpolação de dados de altitude a partir de 169 amostras coletadas em uma área de 1000 x 1000 metros. A metodologia geoestatística foi empregada para analisar a estrutura de dependência espacial dos dados e gerar um mapa de predição contínuo para toda a área de estudo, bem como um mapa de incerteza das predições.

Preparação de Ambiente e Carregamento dos Dados

Os dados são carregados a partir do arquivo Dados_altitude.txt.

library(sp)
library(gstat)
library(ggplot2)
library(dplyr)
library(readr)

dados <- read.table("Dados_altitude.txt", header = TRUE, row.names = 1)
dados <- dados %>% select(X, Y, altitude)
head(dados)
##           X Y  altitude
## 1   0.00000 0  94.99051
## 2  83.33333 0  98.88930
## 3 166.66667 0 103.57085
## 4 250.00000 0 106.03346
## 5 333.33333 0 104.57286
## 6 416.66667 0 101.79335

Análise Descritiva

A primeira etapa consiste em entender a distribuição estatística da variável de interesse, altitude.

summary(dados$altitude)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   94.99  110.89  120.33  121.09  129.20  152.45
sd(dados$altitude)
## [1] 13.61839
var(dados$altitude)
## [1] 185.4604

Os dados de altitude são consistentes, com uma tendência central bem definida em torno de 121 metros e uma variabilidade moderada. A distribuição aparenta ser simétrica.

Histograma

O histograma permite visualizar a forma da distribuição da variável.

ggplot(dados, aes(x = altitude)) +
  geom_histogram(aes(y = ..density..), binwidth = 5, fill = "lightblue", color = "black") +
  geom_density(alpha = .2, fill = "#FF6666") +
  labs(title = "Histograma de Frequência da Altitude", x = "Altitude (m)", y = "Densidade") +
  theme_minimal()

O histograma oferece uma confirmação visual de que a variável altitude é bem comportada. A distribuição é unimodal (possui um único pico) e simétrica

Análise Espacial e Estacionariedade

Nesta seção, investigamos como os dados se distribuem no espaço e avaliamos a premissa de estacionariedade de primeira ordem.

Mapa de Pontos

ggplot(dados, aes(x = X, y = Y)) +
  geom_point(aes(color = altitude, size = altitude)) +
  scale_color_viridis_c() +
  labs(title = "Distribuição Espacial das Amostras de Altitude",
       x = "Coordenada X (m)", y = "Coordenada Y (m)",
       color = "Altitude (m)", size = "Altitude (m)") +
  coord_fixed(ratio = 1) +
  theme_minimal()

A distribuição dos valores de altitude não é aleatória. Existe um padrão espacial muito bem definido. As altitudes mais baixas (cores escuras, roxo/azul) estão concentradas na porção sudoeste do mapa (canto inferior esquerdo), enquanto as altitudes mais altas (cores claras, amarelo) se concentram na porção nordeste (canto superior direito).

Estacionariedade de Primeira Ordem

Para confirmar a tendência, plotamos a altitude contra cada uma das coordenadas.

ggplot(dados, aes(x = X, y = altitude)) + geom_point() + geom_smooth(method = "lm", col = "red") + labs(title="Altitude vs. Coordenada X")

ggplot(dados, aes(x = Y, y = altitude)) + geom_point() + geom_smooth(method = "lm", col = "red") + labs(title="Altitude vs. Coordenada Y")

  • No gráfico Altitude vs. Coordenada X, a linha de tendência vermelha sobe da esquerda para a direita, indicando que, em média, a altitude aumenta à medida que nos movemos na direção Leste.

  • No gráfico Altitude vs. Coordenada Y, a linha de tendência também é ascendente, mostrando que, em média, a altitude aumenta à medida que nos movemos na direção Norte.

A estacionariedade de primeira ordem pressupõe que a média da variável é constante em toda a área de estudo. Como os gráficos demonstram que a média da altitude depende claramente das coordenadas X e Y, a premissa de estacionariedade de primeira ordem é violada.

Modelagem da Estrutura de Dependência Espacial

O semivariograma é a ferramenta utilizada para modelar a autocorrelação espacial dos dados após a remoção da tendência.

Análise de Anisotropia

Verificamos se a dependência espacial é a mesma em todas as direções.

coordinates(dados) <- ~X+Y
f_trend <- as.formula(altitude ~ X + Y)
v_aniso <- variogram(f_trend, data = dados, alpha = c(0, 45, 90, 135))
plot(v_aniso, main = "Semivariogramas Direcionais")

  • Direção 45° (NE-SO): A semivariância cresce de forma mais lenta. Os pontos demoram mais para atingir o patamar. Isso indica uma maior continuidade espacial nesta direção.

  • Direção 135° (SE-NO): A semivariância neste painel sobe muito mais rapidamente, atingindo valores mais altos em distâncias mais curtas. Isso indica uma menor continuidade espacial e um alcance menor nesta direção.

  • Direções 0° e 90°: Apresentam um comportamento intermediário entre as direções de 45° e 135°.

A análise revelou a presença de anisotropia, com a maior continuidade espacial ocorrendo na direção Nordeste-Sudoeste (45°), que coincide com a direção da tendência que identificamos anteriormente.

Semivariograma Experimental Isotrópico

v_iso <- variogram(f_trend, data = dados)
plot(v_iso, main = "Semivariograma Experimental Isotrópico (Resíduos)")

  • Efeito Pepita: Olhando para onde a curva de pontos interceptaria o eixo Y (em distância zero), notamos que ela não começa na origem (0,0). O primeiro ponto já tem uma semivariância de aproximadamente 10. Isso indica a presença de um efeito pepita.

  • Patamar: O patamar parece estar se formando em torno de um valor de semivariância de 35 a 40. ## Ajuste de Modelos

  • Alcance: O gráfico parece atingir esse patamar a uma distância de aproximadamente 350 a 400 metros. Na prática, significa que duas amostras separadas por mais de ~400 metros não estão espacialmente correlacionadas

Ajuste de Modelos Teóricos

Esférico

fit_sph <- fit.variogram(v_iso, model = vgm("Sph"))
fit_sph
##   model     psill  range
## 1   Nug  2.985885   0.00
## 2   Sph 37.150624 600.61
plot(v_iso, fit_sph, main = "Ajuste do Modelo Esférico")

  • Efeito Pepita (C_0): 2.99

  • Patamar Parcial (C): 37.15

  • Patamar Total (Sill, C_0+C): 2.99 + 37.15 = 40.14

  • Alcance (a): 600.61 metros

O alcance estimado de 601 metros parece um pouco longo. Observando os pontos experimentais, eles parecem se estabilizar (atingir o patamar) um pouco antes, talvez em torno de 400-450 metros. O modelo esférico está “esticando” a influência espacial um pouco além do que os dados brutos sugerem.

Exponencial

fit_exp <- fit.variogram(v_iso, model = vgm("Exp"))
fit_exp
##   model     psill    range
## 1   Nug  1.656171   0.0000
## 2   Exp 64.120364 556.8719
plot(v_iso, fit_exp, main = "Ajuste do Modelo Exponencial")

  • Efeito Pepita (C_0): 1.66

  • Patamar Parcial (C): 64.12

  • Patamar Total (Sill, C_0+C): 2.99 + 37.15 = 65.78

  • Alcance (a): 556.87 metros

O modelo estima um patamar de quase 66. Isso parece ser uma superestimação significativa quando olhamos para os pontos do semivariograma experimental, que não parecem passar muito de 40.

Gaussiano

fit_gau <- fit.variogram(v_iso, model = vgm("Gau"))
fit_gau
##   model     psill    range
## 1   Nug  7.911924   0.0000
## 2   Gau 28.484601 255.4246
plot(v_iso, fit_gau, main = "Ajuste do Modelo Gaussiano")

  • Efeito Pepita (C_0): 7.91

  • Patamar Parcial (C): 28.48

  • Patamar Total (Sill, C_0+C): 2.99 + 37.15 = 36.40

  • Alcance (a): 255.42 metros

O modelo estima um patamar e um alcance que parecem muito consistentes com os dados experimentais. Seu ponto fraco é o efeito pepita relativamente alto.

Modelo Final

O modelo Gaussiano se destaca por ter o patamar e o alcance mais realistas, mais alinhados com os dados brutos. Apesar do seu efeito pepita ser maior, ele parece capturar a distância e a variância totais de forma mais fidedigna.

Interpolação por Krigagem Universal

Com o modelo de variograma definido, realizamos a interpolação para gerar os mapas finais.

grade <- expand.grid(X = seq(0, 1000, 10), Y = seq(0, 1000, 10))
coordinates(grade) <- ~X+Y
gridded(grade) <- TRUE

krig_ord <- krige(f_trend, dados, grade, model = fit_gau)
## [using universal kriging]
krig_df <- as.data.frame(krig_ord)
names(krig_df) <- c("X", "Y", "predicao", "variancia")

ggplot(krig_df, aes(x = X, y = Y)) +
  geom_tile(aes(fill = predicao)) +
  scale_fill_viridis_c(name = "Valor Estimado") +
  coord_fixed(ratio = 1) +
  labs(title = "Mapa de Predição por Krigagem Universal", x = "Coordenada X (m)", y = "Coordenada Y (m)") +
  theme_minimal()

ggplot(krig_df, aes(x = X, y = Y)) +
  geom_tile(aes(fill = variancia)) +
  scale_fill_viridis_c(name = "Variância", option = "A") +
  coord_fixed(ratio = 1) +
  labs(title = "Mapa da Variância da Krigagem", x = "Coordenada X (m)", y = "Coordenada Y (m)") +
  theme_minimal()

O mapa de predição por Krigagem Universal recria com sucesso a tendência de aumento das altitudes do canto sudoeste (valores mais baixos, em roxo) para o canto nordeste (valores mais altos, em amarelo), que foi consistentemente identificada em todas as etapas anteriores da análise.

A variância é praticamente nula (cores escuras) exatamente sobre os 169 pontos de amostragem, pois a krigagem é um interpolador exato. A incerteza, e portanto a variância, aumenta à medida que nos afastamos dos pontos com informação real.

Conclusão

O processo geoestatístico foi aplicado com sucesso para transformar 169 amostras pontuais de altitude em uma superfície contínua e informativa. A análise detalhada da tendência espacial e da estrutura de autocorrelação permitiu a aplicação da Krigagem Universal, resultando em um mapa de predição robusto e um mapa de variância que quantifica a confiabilidade das estimativas. Os objetivos do trabalho foram alcançados, gerando um modelo digital de elevação útil e com controle de qualidade para a área de estudo.