Atividade sobre Métricas de Paisagem
Luiz Gustavo
25/01/2023
Parte inicial
Primeiramente, deve-se carregar os pacotes que serão utilizados no processo.
#Carregando os pacotes
library(tidyverse)
library(sf)
library(raster)
library(fasterize)
library(landscapemetrics)
library(tmap)
# Definir a pasta de trabalho
setwd("C:/Users/luiz8/Desktop/Metricas")Carregando os dados necessários para a análise. Nota-se que o diretório para download estava fora do ar, então baixei os arquivos diretamente do github e os utilizei.
# Carregando função multifit
source("https://raw.githubusercontent.com/phuais/multifit/master/multifit.R")
# options
options(timeout = 600)
# importar shapefile
uso <- sf::st_read("03_data/SP_3543907_USO.shp", quiet = TRUE)
# tabela de atributos
sf::st_drop_geometry(uso)## GEOCODIGO MUNICIPIO UF CD_UF CLASSE_USO AREA_HA
## 1 3543907 RIO CLARO SP 35 água 357.027
## 2 3543907 RIO CLARO SP 35 área antropizada 37297.800
## 3 3543907 RIO CLARO SP 35 área edificada 5078.330
## 4 3543907 RIO CLARO SP 35 formação florestal 7017.990
## 5 3543907 RIO CLARO SP 35 silvicultura 138.173# plot do mapa de uso e ocupação do solo importado
tm_shape(uso) +
tm_fill(col = "CLASSE_USO", title = "Legenda",
pal = c("blue", "orange", "gray", "forestgreen", "green"))
Rasterizando
Transformando vetor em raster para operações
# criar uma coluna numerica para as classes de uso da terra
uso$classe_num <- factor(uso$CLASSE_USO)
# criar um raster vazio
ra <- fasterize::raster(uso, res = 30)
# rasterizar baseado no shapefile
uso_raster <- fasterize::fasterize(sf = uso, raster = ra, field = "classe_num")
# plotar mapa fasterize
tm_shape(uso_raster) +
tm_raster(style = "cat", title = "Legenda",
palette = c("blue", "orange", "gray", "forestgreen", "green")) +
tm_grid(lines = FALSE, labels.rot = c(0, 90), labels.size = .8) +
tm_compass(position = c(.73, .08)) +
tm_scale_bar(position = c(.63, 0), text.size = .65) +
tm_layout(legend.position = c("left", "bottom")) 
Adicionando os pontos de amostragem
# amostragens de campo
amost <- readr::read_csv("03_data/pontos_amostragem.csv")
# transformar informações do csv em pontos georeferenciados
amost_vetor <- sf::st_as_sf(x = amost, coords = c("x", "y"), crs = raster::crs(uso))
# buffers, informações do tipo POINT são transformados em POLYGON, tornando-se 2d, criando uma área ao redor do ponto
buffers <- sf::st_buffer(x = amost_vetor, dist = 1000)
# mapa
tm_shape(uso_raster) +
tm_raster(style = "cat", title = "Legenda",
palette = c("blue", "orange", "gray", "forestgreen", "green")) +
tm_shape(buffers) +
tm_borders(col = "red", lwd = 2) +
tm_shape(amost_vetor) +
tm_dots(size = .7, shape = 20, alpha = .7) +
tm_grid(lines = FALSE, labels.rot = c(0, 90), labels.size = .8) +
tm_compass(position = c(.73, .08)) +
tm_scale_bar(position = c(.63, 0), text.size = .65) +
tm_layout(legend.position = c("left", "bottom")) 
Formando mapas ajustados a cada local de amostragem do estudo
# crop e mask das paisagens
paisagens <- raster::mask(uso_raster, buffers)
tm_shape(paisagens) +
tm_raster(style = "cat", title = "Legenda",
palette = c("blue", "orange", "gray", "forestgreen", "green")) +
tm_shape(buffers) +
tm_borders(col = "red", lwd = 2) +
tm_shape(amost_vetor) +
tm_dots(size = .7, shape = 20, alpha = .7)
# crop e mask das paisagens individualmente
paisagens_list <- NULL
for(i in 1:10){
# informacao
print(paste0("Ajustando a paisagem ", i))
# filter
buffers_i <- buffers[i, ]
# crop e mask
paisagens_list[[i]] <- uso_raster %>%
raster::crop(buffers_i) %>%
raster::mask(buffers_i)
}## [1] "Ajustando a paisagem 1"
## [1] "Ajustando a paisagem 2"
## [1] "Ajustando a paisagem 3"
## [1] "Ajustando a paisagem 4"
## [1] "Ajustando a paisagem 5"
## [1] "Ajustando a paisagem 6"
## [1] "Ajustando a paisagem 7"
## [1] "Ajustando a paisagem 8"
## [1] "Ajustando a paisagem 9"
## [1] "Ajustando a paisagem 10"# mapas
cores <- data.frame(val = freq(uso_raster)[1:5, 1],
col = c("blue", "orange", "gray", "forestgreen", "green"))
for(i in 1:10){
# nomes
nome_paisagem <- paste0("map_pai", i)
nome_floresta <- paste0("map_for", i)
# mapas
assign(nome_paisagem,
tm_shape(paisagens_list[[i]]) +
tm_raster(style = "cat", legend.show = FALSE,
palette = cores[cores$val %in% freq(paisagens_list[[i]])[, 1], 2]) +
tm_shape(buffers[i, ]) +
tm_borders(col = "red", lwd = 2) +
tm_shape(amost_vetor) +
tm_dots(size = .7, shape = 20, alpha = .7) +
tm_layout(main.title = names(paisagens_list)[i])
)
assign(nome_floresta,
tm_shape(raster::reclassify(x = paisagens_list[[i]], rcl = c(0,3,NA, 3,4,1, 4,6,NA))) +
tm_raster(legend.show = FALSE,
pal = "forestgreen") +
tm_shape(buffers[i, ]) +
tm_borders(col = "red", lwd = 2) +
tm_shape(amost_vetor) +
tm_dots(size = .7, shape = 20, alpha = .7) +
tm_layout(main.title = names(paisagens_list)[i])
)
}
# todos os mapas
tmap::tmap_arrange(map_pai1, map_pai2, map_pai3, map_pai4, map_pai5, map_pai6,
map_pai7, map_pai8, map_pai9, map_pai10)
Mostar apenas as áreas de floresta dos mapas
tmap::tmap_arrange(map_for1, map_for2, map_for3, map_for4, map_for5, map_for6,
map_for7, map_for8, map_for9, map_for10)
Tratando os dados para a geração dos mapas
# checar o raster
landscapemetrics::check_landscape(paisagens_list)## layer crs units class n_classes OK
## 1 1 projected m integer 3 v
## 2 2 projected m integer 3 v
## 3 3 projected m integer 3 v
## 4 4 projected m integer 3 v
## 5 5 projected m integer 3 v
## 6 6 projected m integer 3 v
## 7 7 projected m integer 4 v
## 8 8 projected m integer 3 v
## 9 9 projected m integer 4 v
## 10 10 projected m integer 3 v#' prerequisitos do raster
#' 1. sistema de referencias de coordenadas e projetada (crs)
#' 2. unidade esta em metros (units)
#' 3. classes como valores inteiros (class)
#' 4. numero de classes (n_class)
# metricas
all_metrics <- landscapemetrics::list_lsm()
# patch metrics
patch_metrics <- landscapemetrics::list_lsm() %>%
dplyr::filter(level == "patch") %>%
dplyr::arrange(type)
patch_metrics %>%
dplyr::group_by(type) %>%
dplyr::summarise(n = n())## # A tibble: 4 x 2
## type n
## <chr> <int>
## 1 aggregation metric 1
## 2 area and edge metric 3
## 3 core area metric 3
## 4 shape metric 5# class metrics
class_metrics <- landscapemetrics::list_lsm() %>%
dplyr::filter(level == "class") %>%
dplyr::arrange(type)
class_metrics_type <- class_metrics %>%
dplyr::group_by(type) %>%
dplyr::summarise(n = n())
# landscape metrics
landscape_metrics <- landscapemetrics::list_lsm() %>%
dplyr::filter(level == "landscape") %>%
dplyr::arrange(type)
landscape_metrics_type <- landscape_metrics %>%
dplyr::group_by(type) %>%
dplyr::summarise(n = n())Gerando alguns mapas de floresta, contendo informações de quantidade de manchas presentes e área core dessas manchas.
landscapemetrics::show_patches(landscape = paisagens_list[[1]],
class = 4, directions = 8)## $layer_1
landscapemetrics::show_patches(landscape = paisagens_list[[1]],
class = 4, directions = 4)## $layer_1
Nos mapas acima vemos que, com a mudança da propriedade “directions” em “show_patches”, ocorre uma mudança da quantidade das manchas; com o valor “8” a conexão entre áreas para definir como parte de uma mesma mancha precisa ser por todo o perímetro, incluindo as quinas, já o valor “4” implica uma conexão encima, abaixo e de ambos os lados, sem incluir quinas.
landscapemetrics::show_cores(landscape = paisagens_list[[1]],
class = 4, edge_depth = 1)## $layer_1
landscapemetrics::show_cores(landscape = paisagens_list[[1]],
class = 4, edge_depth = 2)## $layer_1
Já com a mudança no “edge_dept” em “show_core” o tamanho da região core, a qual é a mais preservada e conservada das características originais da floresta mudou. Assim, vemos 2 análises, uma mais otimista e a outra menos otimista sobre o impacto do efeito de borda.
landscapemetrics::show_lsm(landscape = paisagens_list[[1]],
what = "lsm_p_area", class = 4)## $layer_1
No mapa acima é calculada a área das manchas de floresta.
Vamos utilizar essa ideia para comparar duas outras áreas de amostragem:
landscapemetrics::show_patches(landscape = paisagens_list[[8]],
class = 4, directions = 8)## $layer_1
landscapemetrics::show_patches(landscape = paisagens_list[[9]],
class = 4, directions = 4)## $layer_1
landscapemetrics::show_cores(landscape = paisagens_list[[8]],
class = 4, edge_depth = 2)## $layer_1
landscapemetrics::show_cores(landscape = paisagens_list[[9]],
class = 4, edge_depth = 2)## $layer_1
Vemos que estamos comparando duas paisagens bem diferentes, a paisagem 8 é formada por uma grande mancha de floresta com uma grande área core em comparação à paisagem 9, bastante fragmentada e com uma área core bastante pequena mesmo considerando a soma de todos os frangmentos.
Cálculo das métricas
Área
Serão geradas tabelas com a área de cada “bloco” da paisagem, cada conjunto de células similares, de acordo com cada nível.
# area no nivel de mancha (patch - p)
paisagens_area_p <- landscapemetrics::lsm_p_area(landscape = paisagens_list[[1]])
paisagens_area_p_all <- landscapemetrics::lsm_p_area(landscape = paisagens_list)
# area no nivel de classe (class - c)
paisagens_area_c <- landscapemetrics::lsm_c_area_mn(landscape = paisagens_list[[1]])
paisagens_area_c_all <- landscapemetrics::lsm_c_area_mn(landscape = paisagens_list)
# area no nivel de paisagem (landscape - l)
paisagens_area_l <- landscapemetrics::lsm_l_area_mn(landscape = paisagens_list[[1]])
paisagens_area_l_all <- landscapemetrics::lsm_l_area_mn(landscape = paisagens_list)Vamos analisar individualmente a paisagem 1 e o que significa a tabela para cada nível:
paisagens_area_p## # A tibble: 32 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 patch 1 1 area 2.07
## 2 1 patch 1 2 area 0.09
## 3 1 patch 1 3 area 0.18
## 4 1 patch 1 4 area 2.07
## 5 1 patch 1 5 area 0.09
## 6 1 patch 2 6 area 247.
## 7 1 patch 2 7 area 0.09
## 8 1 patch 2 8 area 0.27
## 9 1 patch 2 9 area 0.09
## 10 1 patch 2 10 area 0.09
## # ... with 22 more rowsMostra todos os blocos de célula individualmente, a classe diz respeito ao uso e ocupação: 1 - água, 2 - área antropizada, 3 - área edificada, 4 - formação florestal e 5 - silvicultura.
paisagens_area_c## # A tibble: 3 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 class 1 NA area_mn 0.9
## 2 1 class 2 NA area_mn 25.1
## 3 1 class 4 NA area_mn 3.47Mostra a soma de todos os valores dentro de uma classe.
paisagens_area_l## # A tibble: 1 x 6
## layer level class id metric value
## <int> <chr> <int> <int> <chr> <dbl>
## 1 1 landscape NA NA area_mn 9.82Mostra a média de todos os agrupamentos de células.
Calcular todas as métricas, por nível
Calcula não apenas área, mas também: patch area; core area index; related circumscribing circle; contiguity index; core area; euclidean nearest neighbor distance; fractal dimension index; radius of gyration; number of core areas; perimeter-area ratio; patch perimeter; shape index; entre outras.
Primeiramente, vemos o calculo a nível de macha, para cada mancha na paisagem:
# patch level
lsm_patch <- landscapemetrics::calculate_lsm(landscape = paisagens_list[[4]],
level = "patch",
edge_depth = 1, # celulas
neighbourhood = 8, # oito celulas nas vizinhancas
full_name = TRUE,
verbose = TRUE,
progress = TRUE)##
> Progress nlayers: 1 / 1Então o nível de classe, que resume todas as manchas de uma mesma classe:
# class level
lsm_class <- landscapemetrics::calculate_lsm(landscape = paisagens_list[[4]],
level = "class",
edge_depth = 1, # celulas
neighbourhood = 8, # oito celulas nas vizinhancas
full_name = TRUE,
verbose = TRUE,
progress = TRUE)##
> Progress nlayers: 1 / 1E o nível paisagem que resume em um valor todas as classes:
# landscape level
lsm_landscape <- landscapemetrics::calculate_lsm(landscape = paisagens_list[[4]],
level = "landscape",
edge_depth = 1, # celulas
neighbourhood = 8, # oito celulas nas vizinhancas
full_name = TRUE,
verbose = TRUE,
progress = TRUE)##
> Progress nlayers: 1 / 1Espacializando os dados
# reclassificar
paisagen04_forest <- raster::reclassify(x = paisagens_list[[4]], rcl = c(0,3,NA, 3,4,1))
tm_shape(paisagen04_forest) +
tm_raster(pal = "forestgreen", legend.show = FALSE)
# calcular e espacializar
paisagen04_forest_patch <- landscapemetrics::spatialize_lsm(paisagen04_forest,
what = "patch",
progress = TRUE)##
> Progress nlayers: 1 / 1Mapas de métricas categóricas da paisagem 4
Mapa de área:
tm_shape(paisagen04_forest_patch$layer_1$lsm_p_area) +
tm_raster(pal = "viridis", title = "Área (ha)")
Mapa de forma:
tm_shape(paisagen04_forest_patch$layer_1$lsm_p_shape) +
tm_raster(pal = "viridis", title = "Formato")
Mapa de distância entre manchas:
tm_shape(paisagen04_forest_patch$layer_1$lsm_p_enn) +
tm_raster(pal = "viridis", title = "Distância (m)")
Múltiplas escalas
Aqui se consegue aumentar ou diminuir a escala do efeito que está sendo analizado na área desde o ponto de amostragem.
# tamanhos
tamanhos <- c(100, 200, 500, 1000)
# multiplos buffers
buffers_multi <- sf::st_buffer(x = sf::st_as_sf(purrr::map_dfr(amost_vetor, rep, 4)),
dist = rep(tamanhos, each = 10))
buffers_multi$buffer <- rep(tamanhos, each = 10)
# map
tm_shape(uso_raster) +
tm_raster(style = "cat", title = "Legenda", alpha = .5,
palette = c("blue", "orange", "gray", "forestgreen", "green")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") 
# metricas multiplas escalas
metricas_multiplas_escalas <- tamanhos %>%
set_names() %>%
map_dfr(~sample_lsm(landscape = uso_raster,
y = amost_vetor,
shape = "circle",
size = .,
what = c("lsm_c_np", "lsm_c_area_mn", "lsm_c_pland"),
all_classes = TRUE,
return_raster = TRUE,
verbose = TRUE,
progress = TRUE),
.id = "buffer")##
> Progress nlayers: 1 / 1
##
> Progress nlayers: 1 / 1
##
> Progress nlayers: 1 / 1
##
> Progress nlayers: 1 / 1# map
tm_shape(paisagens_list[[1]]) +
tm_raster(style = "cat", title = "Legenda", alpha = .7,
palette = c("blue", "orange", "forestgreen")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") +
tm_shape(amost_vetor[1, ]) +
tm_bubbles()
tm_shape(metricas_multiplas_escalas[metricas_multiplas_escalas$layer == 1 &
metricas_multiplas_escalas$buffer == 100, ]$raster_sample_plots[[1]],
bbox = buffers_multi[31,]) +
tm_raster(style = "cat", title = "Legenda", alpha = .7,
palette = c("blue", "orange", "forestgreen")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") +
tm_shape(amost_vetor[1, ]) +
tm_bubbles()
tm_shape(metricas_multiplas_escalas[metricas_multiplas_escalas$layer == 1 &
metricas_multiplas_escalas$buffer == 200, ]$raster_sample_plots[[1]],
bbox = buffers_multi[31,]) +
tm_raster(style = "cat", title = "Legenda", alpha = .7,
palette = c("blue", "orange", "forestgreen")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") +
tm_shape(amost_vetor[1, ]) +
tm_bubbles()
tm_shape(metricas_multiplas_escalas[metricas_multiplas_escalas$layer == 1 &
metricas_multiplas_escalas$buffer == 500, ]$raster_sample_plots[[1]],
bbox = buffers_multi[31,]) +
tm_raster(style = "cat", title = "Legenda", alpha = .7,
palette = c("blue", "orange", "forestgreen")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") +
tm_shape(amost_vetor[1, ]) +
tm_bubbles()
tm_shape(metricas_multiplas_escalas[metricas_multiplas_escalas$layer == 1 &
metricas_multiplas_escalas$buffer == 1000, ]$raster_sample_plots[[1]],
bbox = buffers_multi[31,]) +
tm_raster(style = "cat", title = "Legenda", alpha = .7,
palette = c("blue", "orange", "forestgreen")) +
tm_shape(buffers_multi) +
tm_borders(col = "red") +
tm_shape(amost_vetor[1, ]) +
tm_bubbles()
Podemos ver que os mapas gerados tem o mesmo ponto central, mas o buffer é diferente, permitindo uma escala maior ou menor da paisagem.
Multifit
A função multifit foi utilizada para seleção da melhor escala para se medir o efeito que alguns parâmetros têm sobre a amostra. O número de manchas, a área, a área core foram analizada.
# numero de especies por paisagem
n_sp <- tibble::tibble(id = 1:10, s = c(5, 3, 6, 5, 3, 2, 1, 10, 7, 4))
# preparar os dados numero de manchas
data_np <- metricas_multiplas_escalas %>%
dplyr::filter(metric == "np",
class == 4) %>%
dplyr::mutate(value = ifelse(is.na(value), 0, value)) %>%
tidyr::pivot_wider(id_cols = plot_id,
names_from = c(metric, buffer),
values_from = value) %>%
dplyr::right_join(n_sp, ., by = c("id" = "plot_id"))
fits_np <- multifit(data = data_np,
mod = "glm",
args = c("family = poisson"),
multief = colnames(data_np)[3:6],
formula = s ~ multief,
criterion = "AIC",
plot_est = TRUE)
# preparar os dados area
data_area <- metricas_multiplas_escalas %>%
dplyr::filter(metric == "area_mn",
class == 4) %>%
dplyr::mutate(value = ifelse(is.na(value), 0, value)) %>%
tidyr::pivot_wider(id_cols = plot_id,
names_from = c(metric, buffer),
values_from = value) %>%
dplyr::right_join(n_sp, ., by = c("id" = "plot_id"))
fits_area <- multifit(data = data_area,
mod = "glm",
multief = colnames(data_area)[3:6],
formula = s ~ multief,
criterion = "AIC",
plot_est = TRUE)
# preparar os dados pland
data_pland <- metricas_multiplas_escalas %>%
dplyr::filter(metric == "pland",
class == 4) %>%
dplyr::mutate(value = ifelse(is.na(value), 0, value)) %>%
tidyr::pivot_wider(id_cols = plot_id,
names_from = c(metric, buffer),
values_from = value) %>%
dplyr::right_join(n_sp, ., by = c("id" = "plot_id"))
fits_pland <- multifit(data = data_pland,
mod = "glm",
multief = colnames(data_pland)[3:6],
formula = s ~ multief,
criterion = "AIC",
plot_est = TRUE)
Vimos então que as melhores escolhas são: para o número de manchas - 100, para a área - 200, para a área core - 500.
# modelos finais
fits_np$models$np_100##
## Call: glm(formula = s ~ np_100, family = poisson, data = data_np)
##
## Coefficients:
## (Intercept) np_100
## 1.3455 0.2817
##
## Degrees of Freedom: 9 Total (i.e. Null); 8 Residual
## Null Deviance: 13.63
## Residual Deviance: 9.4 AIC: 45.85fits_area$models$area_mn_200##
## Call: glm(formula = s ~ area_mn_200, data = data_area)
##
## Coefficients:
## (Intercept) area_mn_200
## 3.4247 0.5665
##
## Degrees of Freedom: 9 Total (i.e. Null); 8 Residual
## Null Deviance: 62.4
## Residual Deviance: 21.29 AIC: 41.94fits_pland$models$pland_500##
## Call: glm(formula = s ~ pland_500, data = data_pland)
##
## Coefficients:
## (Intercept) pland_500
## 2.71793 0.08292
##
## Degrees of Freedom: 9 Total (i.e. Null); 8 Residual
## Null Deviance: 62.4
## Residual Deviance: 13.48 AIC: 37.36broom::tidy(fits_np$models$np_100)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 1.35 0.180 7.48 7.49e-14
## 2 np_100 0.282 0.128 2.20 2.80e- 2broom::tidy(fits_area$models$area_mn_200)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 3.42 0.596 5.74 0.000432
## 2 area_mn_200 0.567 0.144 3.93 0.00436broom::tidy(fits_pland$models$pland_500)## # A tibble: 2 x 5
## term estimate std.error statistic p.value
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 (Intercept) 2.72 0.539 5.04 0.000998
## 2 pland_500 0.0829 0.0154 5.39 0.000655# aicc
aic <- bbmle::ICtab(fits_np$models$np_100,
fits_area$models$area_mn_200,
fits_pland$models$pland_500,
type = "AICc",
weights = TRUE,
delta = TRUE,
logLik = TRUE,
sort = TRUE,
nobs = 10)Gráficos
Plotar gráficos que buscam uma relação entre a porcentagem de habitat (área core) e o número de espécies e entre a área média das manchas e o número de espécies. Observamos que há uma correlação positiva entre ambos.
# graficos
ggplot(data = data_pland) +
aes(x = pland_500, y = s) +
stat_smooth(method = "glm", method.args = list(family = "poisson"), col = "black", level = .95) +
geom_point(shape = 21, size = 5, col = "black", fill = "blue", alpha = .8) +
theme_classic() +
labs(x = "Porcentagem de habitat (%) - 500 m", y = "Número de espécies") +
theme(axis.title = element_text(size = 24),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))
ggplot(data = data_area) +
aes(x = area_mn_200, y = s) +
stat_smooth(method = "glm", method.args = list(family = "poisson"), col = "black", level = .95) +
geom_point(shape = 21, size = 5, col = "black", fill = "forestgreen", alpha = .8) +
theme_classic() +
labs(x = "Área média das manchas (ha) - 200 m", y = "Número de espécies") +
theme(axis.title = element_text(size = 24),
axis.text.x = element_text(size = 20),
axis.text.y = element_text(size = 20),
legend.title = element_text(size = 14),
legend.text = element_text(size = 12))