Unidade I – Modelos de Variação Espacial
Temas: A variação espacial e os fatores que a determinam. Introdução à modelagem de dados espaciais. O modelo discreto de explicação da variação espacial. As visões determinística e estocástica do modelo discreto de variação espacial. O modelo contínuo de variação espacial. A visão determinística do modelo contínuo de variação espacial.
1 Introdução à modelagem espacial
Natureza da variação espacial
2 Modelo discreto de variação espacial
2.1 Modelo discreto: visão determinística
# Pacotes
library(magrittr)
library(dplyr)
library(glue)
library(lattice)
# 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)\[Solo = f(cl, o, r, p, t, \ldots)\]
Vejamos um exemplo de mapa pedológico produzido usando o modelo discreto de variação espacial em sua visão determinística. Para isso será preciso 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)
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 slotO próximo passo consiste em produzir uma figura de pedologia25. Para isso usamos a função sp::spplot, criada especificamente para a produção de figuras de objetos espaciais. Antes disso, vamos criar uma rampa de cores específica para pedologia25, a qual será constituída de tantas cores quantas foram as unidades de mapeamento. Vejamos o resultado na Figura 2.1.
col_soil_um <- terrain.colors(nlevels(pedologia25$um))
sp::spplot(pedologia25, scales = list(draw = TRUE), col.regions = col_soil_um)
Figure 2.1: Mapa pedológico preparado na escala cartográfica de 1:25000 e suas nove unidades de mapeamento.
Descarregar dados do conjunto de dados ctbo0030:
perfis25_o <- febr::observations('ctb0030', which.cols = 'all')##
|
| | 0%
|
|=================================================================| 100%perfis25_l <- febr::layers('ctb0030', which.cols = 'all', missing.data = 'keep')##
|
| | 0%
|
|=================================================================| 100%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)Criar objeto espacial com os perfis e gerar figura.
sp::coordinates(perfis25) <- c('coord_x', 'coord_y')
sp::proj4string(perfis25) <- sirgas2000
perfis25 <- sp::spTransform(perfis25, wgs84utm22s)
sp::spplot(
pedologia25, scales = list(draw = TRUE),
xlim = extendrange(c(perfis25@bbox[1, ], pedologia25@bbox[1, ])),
ylim = extendrange(c(perfis25@bbox[2, ], pedologia25@bbox[2, ])),
col.regions = terrain.colors(nlevels(pedologia25$um)),
main = "Localização dos perfis modais") +
lattice::xyplot(coord_y ~ coord_x, data = as.data.frame(perfis25@coords),
pch = 20, col = 'red', lwd = 2, cex = 2) %>%
latticeExtra::as.layer()
## observacao_id taxon_sibcs_2009 ca_kcl_aas
## 1 PERFIL-01 Argissolo Bruno-Acinzentado Alítico abrúptico 0.00
## 2 PERFIL-02 Neossolo Litólico Distro-úmbrico típico 4.30
## 3 PERFIL-03 Argissolo Vermelho Alítico típico 1.00
## 4 PERFIL-04 Cambissolo Háplico Ta Eutrófico típico 5.30
## 5 PERFIL-05 Neossolo Litólico Eutro-úmbrico típico 3.80
## 6 PERFIL-06 Planossolo Háplico Alítico típico 4.40
## 7 PERFIL-07 Neossolo Flúvico Tb Eutrófico fragmentário 7.60
## 8 PERFIL-08 Neossolo Regolítico Distro-Úmbrico típico 1.44
## 9 PERFIL-09 Neossolo Regolítico Distro-Úmbrico típico 2.05
## 10 PERFIL-10 Neossolo Regolítico Distro-Úmbrico típico 2.71
## 11 PERFIL-11 Neossolo Litólico Distro-Úmbrico típico 2.35
## areia_ argila_
## 1 823 82
## 2 266 303
## 3 199 384
## 4 843 84
## 5 847 73
## 6 60 393
## 7 750 100
## 8 763 107
## 9 669 127
## 10 635 107
## 11 534 2562.2 Modelo discreto: visão estocástica
pontos400_o <- febr::observations('ctb0003', which.cols = 'all')##
|
| | 0%
|
|=================================================================| 100%pontos400_l <- febr::layers('ctb0003', which.cols = 'all', missing.data = 'keep')##
|
| | 0%
|
|=================================================================| 100%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)
sp::spplot(
pedologia25, scales = list(draw = TRUE),
col.regions = terrain.colors(nlevels(pedologia25$um)),
main = "Localização das 400 observações") +
lattice::xyplot(coord_y ~ coord_x,
data = pontos400@coords %>% as.data.frame(), pch = 21, col = 'red') %>%
latticeExtra::as.layer()
Nós podemos calcular a média e desvio padrão da variável de interesse para cada categoria a partir das observações do solo que estão dentro de cada categoria. Mas primeiro temos que identificar a categoria dentro da qual cada observação do solo se encontra. Para isso usamos a função sp::over e armazenamos o resultado em uma nova coluna de pontos400 chamada um. Assumindo que nossa amostra é grande o suficiente, que as observações são independetes, e que a variável de interesse possui distribuição nomal, calculamos o erro padrão (EP) da média para obter os limites inferior e superior do intervalo de confiança de 95%.
Calcular estatísticas de uma variável do solo para cada unidade de mapeamento.
pontos400$um <- sp::over(x = pontos400, y = pedologia25) %>% unlist()
pontos400@data %>%
filter(!is.na(um)) %>%
group_by(um) %>%
summarise(
media = mean(areia_),
ic95 = 1.96 * sd(areia_) / sqrt(n()),
n = n()) %>%
mutate(
inf = round(media - ic95),
media = round(media),
sup = round(media + ic95)) %>%
select(um, inf, media, sup, n)## # A tibble: 9 x 5
## um inf media sup n
## <fctr> <dbl> <dbl> <dbl> <int>
## 1 Argissolo Bruno-Acinzentado 728 749 771 80
## 2 Argissolo Vermelho 221 238 254 38
## 3 Cambissolo - Neossolo 262 304 346 49
## 4 Neossolo Flúvico 768 864 959 2
## 5 Neossolo Litólico 513 552 591 169
## 6 Neossolo Litólico - Neossolo Regolítico 227 261 295 13
## 7 Neossolo Regolítico 531 600 670 8
## 8 Neossolo Regolítico 813 859 905 4
## 9 Planossolo Háplico 320 465 609 13Vamos comparar o resultado acima com os valores obtidos diretamente dos perfis modais. Será que os dados produzidos a partir dos perfis modais encontram-se dentro do intervalo de confiança de 95% da média estimada acima?
perfis25@data %>%
select(taxon_sibcs_2009, areia_) %>%
arrange(taxon_sibcs_2009)## taxon_sibcs_2009 areia_
## 1 Argissolo Bruno-Acinzentado Alítico abrúptico 823
## 2 Argissolo Vermelho Alítico típico 199
## 3 Cambissolo Háplico Ta Eutrófico típico 843
## 4 Neossolo Flúvico Tb Eutrófico fragmentário 750
## 5 Neossolo Litólico Distro-úmbrico típico 266
## 6 Neossolo Litólico Distro-Úmbrico típico 534
## 7 Neossolo Litólico Eutro-úmbrico típico 847
## 8 Neossolo Regolítico Distro-Úmbrico típico 763
## 9 Neossolo Regolítico Distro-Úmbrico típico 669
## 10 Neossolo Regolítico Distro-Úmbrico típico 635
## 11 Planossolo Háplico Alítico típico 603 Modelo contínuo de variação espacial
3.1 Modelo contínuo: visão determinística
Modelos digitais de elevação: http://www.scielo.br/img/revistas/rbcs/v40//0100-0683-rbcs-18069657rbcs20150022-gf02.jpg
Perfil topográfico de modelos digitais de elevação: http://www.scielo.br/img/revistas/rbcs/v40//0100-0683-rbcs-18069657rbcs20150022-gf03.jpg
Interpolação determinística usando o inverso da distância como peso da influência das observações vizinhas.
grid <- sp::spsample(pedologia25, 10000, type = 'regular')
modelo_idw <- gstat::gstat(
id = 'argila', formula = argila_ ~ 1,
data = pontos400, nmax = 8, set = list(idp = 0.5))
mapa_idw <- predict(modelo_idw, grid)## [inverse distance weighted interpolation]sp::gridded(mapa_idw) <- TRUE
sp::spplot(mapa_idw, 'argila.pred', col.regions = col_soil_var,
main = "Interpolação determinística") +
latticeExtra::as.layer(lattice::xyplot(
coord_y ~ coord_x, data = pontos400@coords %>% as.data.frame(),
pch = 21, col = 'red', cex = 0.5))