Temas: Observações dependentes e preferenciais. Resíduos espacialmente correlacionados. A visão estocástica do modelo contínuo de variação espacial. O processo gerador aleatório e espacialmente autocorrelacionado. Semivariância e modelos matemáticos do semivariograma empírico.

Configuração inicial.

# Pacotes
library(magrittr)
library(dplyr)
library(glue)
library(lattice)
library(latticeExtra)

# Sistemas de referência de coordenadas (Fonte: http://spatialreference.org/ref/epsg/)
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sirgas2000 <- sp::CRS('+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')

# Rampas de cores
col_soil_var <- topo.colors(100)

Vamos descarregar os arquivos do mapa pedológico para um diretório local. O nome desse diretório é definido abaixo pelo objeto data_folder. Altere o caminho para esse diretório conforme necessário. Caso você já tenha descarregado os arquivos do mapa pedológico, então data_folder deve ser o caminho para o diretório onde esses arquivos se encontram.

data_folder <- '../data/'
ext <- c('dbf', 'prj', 'shp', 'shx')
files <- glue('pedologia25.{ext}')
download <- !all(files %in% list.files(data_folder)) 
if (download) {
  url <- 'https://github.com/samuel-rosa/UFSM-SOL-843/tree/master/data/'
  url <- glue('{url}{files}')
  destfile <- glue('{data_folder}{files}')
  for (i in 1:length(files)) {
    download.file(url = url[i], destfile = destfile[i])
  }
}

Agora você já pode carregar o mapa pedológico para o seu ambiente de trabalho. Para isso, use a função raster::shapefile. Note que se o sistema operacional do seu computador for o MS Windows, pode ser necessário incluir o o argumento encoding = 'UTF-8' na função raster::shapefile para garantir que os carácteres especiais usados nos nomes das unidades de mapeamento seja reconhecidos. Verifique se a estrutura de pedologia25 está conforme o esperado para um objeto do tipo SpatialPolygonsDataFrame. Note que a função sp::spTransform é usada para projetar as coordenadas original no plano cartesiano (UTM).

pedologia25 <- 
  glue('{data_folder}pedologia25.shp') %>% 
  raster::shapefile(stringsAsFactors = TRUE) %>% 
  sp::spTransform(wgs84utm22s)
col_soil_um <- terrain.colors(nlevels(pedologia25$um))
str(pedologia25, 2)
## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  25 obs. of  1 variable:
##   ..@ polygons   :List of 25
##   ..@ plotOrder  : int [1:25] 12 8 11 5 6 18 24 10 22 23 ...
##   ..@ bbox       : num [1:2, 1:2] 226655 6715750 231890 6721725
##   .. ..- attr(*, "dimnames")=List of 2
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot

1 Introdução à amostragem espacial

A amostragem é um dos maiores contribuintes para os custos da modelagem do ambiente. Assim, o adequado planejamento amostral é essencial para reduzir a necessidade de recursos (financeiros, psicológicos, humanos, entre outros) e maximizar o número de observações possíveis.

A situação ideal é aquela em que os recursos disponíveis não impõe quaisquer limitações à amostragem. Nesse caso seria possível fazer observações e coletar amostras em diversas etapas. Por exemplo, iniciaríamos com um levantamento exploratório para identificar a estrutura da variação espacial da variável de interesse. De posse desse novo conhecimento, executaríamos uma nova etapa amostral a fim de atender a algum critério como, por exemplo, obter uma cobertura espacial aproximadamente uniforme da área sendo modelada. Caso os resultados ainda sejam insatisfatórios, uma nova etapa amostral poderia ser executada para, por exemplo, fazer observações em locais específicos cujas condições ambientais estejam sub-representadas na base de dados. Por fim, calibrado o modelo espacial do solo e feitas as predições espaciais, coletaríamos amostras para validação do modelo preditivo em número suficiente para garantir um nível de confiança pré-determinado.

Mas a situação ideal está longe de ser o que acontece na prática. Em geral temos que fazer todas as observações em uma única fase, coletando todo o material possível, inclusive para as amostras de validação. Isso requer que o tipo de amostragem mais apropriado seja utilizado a fim de otimizar o uso dos recursos disponíveis e obter o melhor modelo preditivo possível. E qual seria o tipo de amostragem mais apropriado? Uma resposta universal à essa pergunda continua desconhecida. A melhor estratégia costuma ser avaliar os diferentes tipos de amostragem frente (1) os objetivos do projeto, (2) os recursos disponíveis, e (3) as dificuldades operacionais encontradas na área sendo modelada.

Podemos dizer que existem dois tipos fundamentais de amostragem espacial:

1.1 Amostragem probabilística

A característica fundamental da amostragem probabilística é que a chance de um determinado local ser amostrado (probabilidade de inclusão) é conhecida e maior do que zero. Em outras palavras, todo e qualquer local possui alguma chance de ser amostrado, mesmo que alguns tenham maior chance do que outros. Um local que não pode ser amostrado tem probabilidade de inclusão igual a zero.

A amostragem probabilística é muito utilizada em experimentos controlados, como aqueles desenvolvidos em campos experimentais, casas de vegetação e laboratórios. No caso da modelagem espacial, a amostragem probabilística costuma ser usada para a validação das predições espaciais. Entretanto, ela também pode ser usada para obter observações para a calibração dos modelos preditivos. São exemplo a amostragem aleatória simples e a amostragem aleatória estratificada simples (Figura 1.1).

pts_random <- list()
for (i in 1:2) {
  pts_random[[i]] <- sp::spsample(x = pedologia25, n = 20, type = "random")
}
sp::spplot(
  pedologia25, col.regions = col_soil_um, alpha.regions = 0.3, colorkey = FALSE) +
  as.layer(xyplot(y ~ x, data = as.data.frame(pts_random[[1]]@coords), col = 'red', pch = 17, cex = 1.5)) +
  as.layer(xyplot(y ~ x, data = as.data.frame(pts_random[[2]]@coords), col = 'blue', pch = 17, cex = 1.5))
Dois conjuntos de pontos de observação espacial com localização selecionada de maneira aleatória.

Figure 1.1: Dois conjuntos de pontos de observação espacial com localização selecionada de maneira aleatória. 

prob <- rep(0.1, 10)
x <- rnorm(10)
mean(x)

prob <- rnorm(10, 100, 2) 
prob <- prob / sum(prob)
sum(prob)

sum(x * prob)

1.2 Amostragem não-probabilística

Na amostragem não-probabilística, como o próprio nome já diz, não são considerados os valores de probabilidade de inclusão para a seleção dos locais de amostragem. A escolha dos locais de amostragem depende da definição de um critério a ser atendido. A amostragem não-probabilística costuma ser dividida em três categorias:

  • casual,
  • conveniente, e
  • intencional.

Na amostragem casual os locais de amostragem são escolhidos, fundamentalmente, em função da subjetividade da pessoa conduzindo a amostragem. Não existe um critério claro a ser atendido. Outros locais amostrais podem ser escolhidos caso outra pessoa conduza a amostragem, mesmo que não haja justificativa plausível para isso.

A amostragem conveniente está diretamente relacionada à otimização do uso dos recursos disponíveis. Ela consiste em evitar realizar observações em locais de difícil acesso como áreas densamente florestadas, distantes de estradas, terrenos íngremes, ou áreas que apresentem risco para a saúde ou à vida devido à presença de, por exemplo, animais peçonhentos. Assim sendo, o critério usado para a definição dos locais de observação é a soma dos custos financeiro e operacional. Quanto menor forem os custos financeiro e operacional, maior será o número de observações.

A amostragem intencional é semelhante à amostragem conveniente no sentido de que em ambas os locais amostrais são definidos a fim de otimizar um critério pré-determinado. A diferença fundamental entre as duas é a natureza desse critério. Enquanto na amostragem conveniente o critério tem origem puramente econômica, na amostragem intencional o critério tem origem pedológica e/ou estatística, podendo-se agregar critérios de origem econômica.

Você já deve ter percebido que, no contexto espacial, a amostragem não-probabilística é mais comum do que a amostragem probabilística. Os dois principais entraves à análise estatística clássica dos dados obtivos via amostragem não-probabilística são a falta de independência entre as observações e a seleção preferencial dos locais de amostragem.

1.2.1 Observações dependentes

Malhas amostrais espaciais podem ser preparadas usando a função sp::spsample. Para preparar uma malha regular, que contenha pontos de observação separados por distâncias idênticas, tanto no eixo x, como no eixo y, usamos o argumento type = 'regular'. Se quisermos que a malha seja exatamente centrada em relação à região amostral, então usamos o argumento offset = c(0.5, 0.5). A Figura 1.2 mostra uma malha regular gerada usando sp::spsample.

pts_regular <- sp::spsample(pedologia25, n = 50, type = "regular", offset = c(0.5, 0.5))
sp::spplot(
  pedologia25, col.regions = col_soil_um, alpha.regions = 0.3, colorkey = FALSE) +
  as.layer(xyplot(x2 ~ x1, data = as.data.frame(pts_regular@coords), col = 'red', pch = 17))
Conjunto de pontos de observação posicionados em malha centrada, de espaçamento fixo e regular em ambas as direções horizontal e vertical.

Figure 1.2: Conjunto de pontos de observação posicionados em malha centrada, de espaçamento fixo e regular em ambas as direções horizontal e vertical.

1.2.2 Observações preferenciais

Como exemplo de amostragem preferencial, vamos analizar novamente os perfis modais do mapa pedológico preparado na escalada cartográfica de 1:25000. Para isso seguimos os mesmos três passos da unidade anterior:

  1. Descarregamos os dados do Repositório Brasileiro de Dados de Ferro do Solo usando febr:observations e febr::layers;
  2. Agregamos os dados das observações e das camadas num data.frame único;
  3. Transformamos esse data.frame para um SpatialPointsDataFrame.
perfis25_o <- febr::observations('ctb0030', which.cols = 'all', progress = FALSE)
perfis25_l <- febr::layers('ctb0030', which.cols = 'all', missing.data = 'keep', progress = FALSE)
id <- c('dataset_id', 'observacao_id')
perfis25 <- 
  merge(x = perfis25_o, y = perfis25_l, by.x = id, by.y = id) %>% 
  select(observacao_id, coord_x, coord_y, taxon_sibcs_2009, ca_kcl_aas, areia_, argila_)
rm(perfis25_l, perfis25_o)
sp::coordinates(perfis25) <- c('coord_x', 'coord_y')
sp::proj4string(perfis25) <- sirgas2000
perfis25 <- sp::spTransform(perfis25, wgs84utm22s)

A Figura 1.3 mostra a localização de 10 dos 11 perfis modais (lembre que um dos perfis modais está localizado para além dos limites da área mapeada).

sp::spplot(
  pedologia25, col.regions = col_soil_um, alpha.regions = 0.3, colorkey = FALSE) +
  as.layer(xyplot(coord_y ~ coord_x, data = as.data.frame(perfis25@coords), col = 'red', pch = 17))
Localização de 10 dos perfis modais usados para caracterização das unidades do mapa pedológico.

Figure 1.3: Localização de 10 dos perfis modais usados para caracterização das unidades do mapa pedológico.

1.2.3 Resíduos correlacionados

Pela 1ª Lei da Geografia sabe-se que objetos mais próximos entre si são mais similares do que objetos que estão distantes uns dos outros. Disso resulta que os resíduos dos modelos são, também, mais parecidos entre si quanto mais próximas entre si estiverem as observações. A existência de correlação espacial faz com que o modelo discreto seja subótimo, portanto menos realista. O mesmo se aplica ao modelo contínuo determinístico, baseado em interpoladores locais, pois os mesmos geralmente ignoram a extensão da correlação espacial e/ou usam critérios heurísticos.

pontos400_o <- febr::observations('ctb0003', which.cols = 'all', progress = FALSE)
pontos400_l <- febr::layers('ctb0003', which.cols = 'all', missing.data = 'keep', progress = FALSE)
id <- c('dataset_id', 'observacao_id')
pontos400 <- 
  merge(pontos400_o, pontos400_l, by.x = id, by.y = id) %>% 
  select(observacao_id, coord_x, coord_y, taxon_sibcs_2009, ca_kcl_aas, argila_, areia_)
rm(pontos400_l, pontos400_o)
sp::coordinates(pontos400) <- ~ coord_x + coord_y
sp::proj4string(pontos400) <- sirgas2000
pontos400 <- sp::spTransform(pontos400, wgs84utm22s)
pontos400$um <- sp::over(x = pontos400, y = pedologia25) %>% unlist()
pontos400in <- pontos400[!is.na(pontos400$um), ]

Figura 1.4.

sp::spplot(
  pedologia25, col.regions = col_soil_um, alpha.regions = 0.3, colorkey = FALSE) +
  latticeExtra::as.layer(
    lattice::xyplot(coord_y ~ coord_x, data = as.data.frame(pontos400in@coords), col = 'red', pch = 17))
Distribuição espacial das 376 observações localizadas dentro dos limites da área de estudo.

Figure 1.4: Distribuição espacial das 376 observações localizadas dentro dos limites da área de estudo.

Como no encontro anterior, vamos estimar a média da variável de interesse para cada categoria (unidade de mapeamento) a partir das observações do solo que estão dentro de cada uma delas. Usamos a função sp::over para identificar a categoria dentro da qual cada observação do solo se encontra. O resultado é armazenado na coluna um. Para facilitar nosso trabalho, vamor criar um novo objeto espacial chamado pontos400in para armazenar apenas as observações que estão no interior (in) dos limites da área de estudo. Depois de estimada a média de cada categoria, podemos calcular o resíduo para cada observação:

\[e(\boldsymbol{s}_{ik}) = y(\boldsymbol{s}_{ik}) - \bar{y}(\boldsymbol{s}_k)\]

A equação acima mostra que o resíduo, \(e(\boldsymbol{s}_{ik})\), de uma observação espacial qualquer, \(y(\boldsymbol{s}_{ik})\), localizada na categoria \(k\), é dado pela diferença entre o valor naquela observação e a média da categoria onde ela se encontra, \(\bar{y}(\boldsymbol{s}_k)\).

pontos400in@data <- 
  pontos400in@data %>% 
  group_by(um) %>%
  mutate(
    residuo = areia_ - mean(areia_))

Calculado os resíduos, vamos agora verificar sua distribuição espacial. Para isso usamos um gráfico de bolhas, o qual pode ser produzido usando a função sp::bubble. Um gráfico de bolha é uma forma de representação gráfica dos valores da variável de interesse em cada observação, onde o tamanho e cor do símbolo são proporcionais ao valor da variável de interesse. Na Figura 1.5, as observações com resíduo negativo são indicadas em magenta, enquanto as observações em verde representam aquelas cujo resíduo é positivo. Quanto maior o diâmetro do círculo, maior é o resíduo em termos absolutos.

sp::bubble(pontos400in, 'residuo', do.sqrt = TRUE, main = "")
Gráfico de bolha dos resíduos da variável de interesse, calculados em função da média de cada categoria do mapa pedológico.

Figure 1.5: Gráfico de bolha dos resíduos da variável de interesse, calculados em função da média de cada categoria do mapa pedológico.

2 Modelo contínuo de variação espacial: visão estocástica

2.1 Processo aleatório espacialmente autocorrelacionado

\[\LaTeX\]

\[Y(\boldsymbol{s})\]

\[y(\boldsymbol{s}_i)\]

2.2 Medida da dependência espacial

Como quantificar a dependência espacial? Podemos usar a variância entre os pares de observações.

\[\sigma^2 = \left\{ y(\boldsymbol{s}_i) - \bar{y}(\boldsymbol{s}) \right\}^2 + \left\{y(\boldsymbol{s}_j) - \bar{y}(\boldsymbol{s})\right\}^2\]

Isso é exatamente o mesmo que:

\[\sigma^2 = \frac{1}{2} \left\{ y(\boldsymbol{s}_i) - y(\boldsymbol{s}_j) \right\}^2\]

Para usar a variância \(\sigma^2\), precisamos assumir que ela depende apenas da distância geográfica que separa as observações \(\boldsymbol{s}_i\) e \(\boldsymbol{s}_j\). Em outras palavras, assumir que a variância não depende da posição absoluta das observações no espaço geográfico.

Também precisamos assumir que a diferença esperada – a média – seja igual à zero, pelo menos dentro de pequenos intervalos de distância entre as observações. Em outras palavras, assumir que, quando avaliada localmente, a média não varia.

Vamos lembrar como é calculada a variância dos dados no contexto não-espacial.

A Figura 2.1 mostra a nuvem semivariográfica dos resíduos da variável de interesse calculados acima.

variograma <- gstat::variogram(areia_ ~ 1, pontos400in, cloud = TRUE, cutoff = Inf)
# mean(variograma$gamma); var(pontos400in$areia_)
plot(variograma, ylab = "Semivariância", xlab = "Distância (m)", ylim = c(0, 400000)) +
  layer(panel.abline(h = res_var, col = 'red', lty = 'dashed')) +
  layer(panel.text(y = res_var, x = 1500, label = "Variância amostral", pos = 3))

# (n * n - 1) / 2
# hist(pontos400in$residuo)
variograma <- gstat::variogram(residuo ~ 1, pontos400in, cloud = TRUE, cutoff = Inf)
res_var <- var(pontos400in$residuo)
plot(variograma, ylab = "Semivariância", xlab = "Distância (m)", ylim = c(0, 400000)) +
  layer(panel.abline(h = res_var, col = 'red', lty = 'dashed')) +
  layer(panel.text(y = res_var, x = 1500, label = "Variância amostral", pos = 3))
Nuvem semivariográfica.

Figure 2.1: Nuvem semivariográfica.

Semivariograma empírico suavizado (Figura ??).

Parâmetros do semivariograma:

  • Variância total (patamar)
  • Variância residual (pepita)
  • Alcance
limites <- seq(0, 2500, length.out = 15)
variograma <- gstat::variogram(residuo ~ 1, pontos400in, boundaries = limites)
plot(variograma, ylab = "Semivariância", xlab = "Distância (m)", pch = 20, cex = 1.5) +
  layer(panel.abline(
    h = c(res_var, min(variograma$gamma)), v = 1500, col = 'red', lty = 'dashed')) +
  layer(panel.text(
    x = c(600, 600, 1500), 
    y = c(res_var, min(variograma$gamma), max(variograma$gamma) / 2), cex = 1.5,
    label = c("Patamar (Variância total)", "Pepita (Variância residual)", 'Alcance'), 
    pos = c(3, 3, 1), srt = c(0, 0, 90))) +
  
  layer(panel.arrows(
    x0 = c(2250, 2250, 0), x1 = c(2250, 2250, 1500), y0 = c(min(variograma$gamma), 0, 18000), 
    y1 = c(res_var, min(variograma$gamma), 18000), length = 0.1, ends = 'both')) +
  
  layer(panel.text(
    x = c(2270, 2270, 750, 0), cex = 1.5,
    y = c(min(variograma$gamma) + diff(range(variograma$gamma)) / 2, 
          min(variograma$gamma) - min(variograma$gamma) / 2, 18000, 2000), 
    label = c('Variância\nexplicada', 'Variância\ninexplicada', 'Correlação\nEspacial', 
              'Erro de medida\nCorrelação de curto alcance'), srt = c(90, 90, 0, 0), 
    pos = c(1, 1, 3, 4)))

Vamos conhecer três funções autorizadas clássicas que podem ser usadas para modelar o semivariograma empírico:

  • Puro ruído (sem estrutura)
  • Exponencial (com patamar e alcance)
  • Gaussiana (com patamar e alcance)
m_exp <- gstat::vgm(psill = 37000 - 8000, model = 'Exp', range = 300, nugget = 8000)
m_gau <- gstat::vgm(psill = 37000 - 12000, model = 'Gau', range = 400, nugget = 12000)
m_nug <- gstat::vgm(psill = 37000, model = 'Nug', range = 0, nugget = 0)
plot(variograma, m_exp, ylab = "Semivariância", xlab = "Distância (m)", pch = 20, cex = 1.5) +
  layer(panel.lines(
    gstat::variogramLine(
      m_exp, dist_vector = seq(0, 2500, length.out = 1500)), col = "darkgreen", lwd = 2)) +
  layer(panel.lines(
    gstat::variogramLine(m_gau, dist_vector = seq(0, 2500, length.out = 1500)), col = "orange", lwd = 2)) +
  layer(panel.lines(
    gstat::variogramLine(m_nug, dist_vector = seq(0, 2500, length.out = 1500)), col = "purple", lwd = 2)) +
  layer(panel.text(
    x = rep(1000, 3), y = seq(10000, 15000, length.out = 3), c('Exponencial', 'Gaussiana','Puro ruído'),
    col = c('darkgreen', 'orange', 'purple'), pos = 4))
Três funções autorizadas para modelar o semivariograma empírico.

Figure 2.2: Três funções autorizadas para modelar o semivariograma empírico.

Vejamos o que cada uma dessas funções representa no mundo real.

grid <- sp::spsample(pedologia25, 10000, type = 'regular')
v <- gstat::gstat(formula = z ~ 1, dummy = TRUE, beta = 0, model = m_gau, nmax = 10)
simula <- predict(v, grid, nsim = 1, debug.level = 0)
names(simula) <- 'Gau'
v <- gstat::gstat(formula = z ~ 1, dummy = TRUE, beta = 0, model = m_exp, nmax = 10)
simula$Exp <- predict(v, grid, nsim = 1, debug.level = 0)$sim1
v <- gstat::gstat(formula = z ~ 1, dummy = TRUE, beta = 0, model = m_nug, nmax = 10)
simula$Nug <- predict(v, grid, nsim = 1, debug.level = 0)$sim1
sp::gridded(simula) <- TRUE
sp::spplot(simula)
Realização de três campos aleatórios definidos por funções gaussiana, exponencial e puro ruído.

Figure 2.3: Realização de três campos aleatórios definidos por funções gaussiana, exponencial e puro ruído.

Vejamos outras funções existentes.

gstat::show.vgms(main = "Modelos do semivariograma", xlab = "Distância", ylab = "Semivariância")

Ajuste do modelo ao variograma empírico

Figura ??

m_exp
##   model psill range
## 1   Nug  8000     0
## 2   Exp 29000   300
m_exp <- gstat::fit.variogram(variograma, m_exp)
plot(variograma, m_exp, col.line = 'red', xlab = 'Distância (m)', ylab = 'Semivariância')

Predições espaciais

mapa_exp <- gstat::krige(residuo ~ 1, pontos400in, grid, m_exp, debug.level = 1)
## [using ordinary kriging]
sp::gridded(mapa_exp) <- TRUE
sp::spplot(mapa_exp, 1, col.regions = col_soil_var)

mapa_exp$var1.var <- sqrt(mapa_exp$var1.var)
sp::spplot(mapa_exp, 2, col.regions = col_soil_var)