Sumário
Mazagão
Darren Norris
2022-11-23
1 Apresentação
O objetivo é calcular métricas de paisagem, descrever a composição e a configuração da paisagem no município de Mazagão.
As métricas da paisagem nos ajudam a entender as mudanças na paisagem de diferentes perspectivas (visual, ecológica, cultural). Asssim sendo, análises com métricas de paisagem é um atividade fundamental na ecologia da paisagem. Nesta exemplo (https://rpubs.com/darren75/mazagao) aprenderemos sobre como analisar a cobertura da terra com métricas de paisagem em R.
1.1 Organização do codigo no tutorial
O tutorial está organizado em etapas de processamento, com blocos de código em caixas cinzas:
codigo de R para executarPara segue os passos, os blocos de código precisam ser executados em sequência. Se você pular uma etapa, ou rodar fora de sequência o próximo bloco de código provavelmente não funcionará.
As linhas de codigo de R dentro de cada caixa tambem preciso ser executados em sequência. O simbolo # é usado para incluir comentarios sobre os passos no codgio. Ou seja, linhas começando com # são ignorados por R, e não é codigo de executar.
# Passo 1
codigo de R passo 1 # texto e numeros tem cores diferentes
# Passo 2
codigo de R passo 2
# Passo 3
codigo de R passo 3Alem disso, os simbolos #> e/ou [1] no início de uma linha indica o resultado que você verá no console de R depois de rodar o codigo, como no proximo exemplo. Digite o código abaixo e veja o resultados (passos 1 a 4).
# Passo 1
1+1[1] 2
# Passo 2
x <- 1+1
# Passo 3
x[1] 2
# Passo 4
x + 1[1] 3
2 Pacotes necessarios
library(landscapemetrics)
library(tidyverse)
library(sf)
library(raster)
library(terra)
library(tmap)
library(gridExtra)
library(kableExtra)
library(leafpop)
library(mapview)
library(mgcv)3 Dados
3.1 Dados: Mazagão
Para entender as mudanças precisamos estabelecer os espaços de interesse. Vamos carregar as camadas necessarios. Baixar o arquivo Link: https://github.com/darrennorris/mazagao/blob/main/vector/mazagao.GPKG . Lembrando-se de salvar o arquivo (“mazagao.GPKG”) em um local conhecido no seu computador.
Agora, com o proximo bloco de codigo, podemos selecionar o arquivo “mazagao.GPKG”, e carregar as camadas.
# Selecionar o arquivo "mazagao.GPKG"
meuSIG <- file.choose()
# distancias ate cidade (250m, 500m, 1km, 2km, 4km, 8km, 16km e 32km)
maza_buffers <- sf::st_read(meuSIG, layer = "maza_buffers")
# municipio
maza_poly <- sf::st_read(meuSIG, layer = "maza_poly")
# ponto
maza_ponto <- sf::st_read(meuSIG, layer = "maza_ponto")Visualizar as camadas
mapview::mapview(maza_buffers, z="dist_km") +
mapview::mapview(maza_ponto, col.regions= "black") 3.2 Dados: MapBiomas
Existem varios formas de importar e exportar dados geoespaciais. Precisamos o arquivo com os dados de MapBiomas referente a região de estudo. Aqui vamos usar dados de 2020 “utm_cover_AP_munis_maza_santana_2020.tif” .
Lembrando-se de salvar o arquivo (“utm_cover_AP_munis_maza_santana_2020.tif”) em um local conhecido no seu computador. Agora, nós podemos carregar os dados de cobertura da terra “utm_cover_AP_munis_maza_santana_2020.tif” com a função rast.
# Selecionar e carregar arquivo "utm_cover_AP_munis_maza_santana_2020.tif"
mapbiomas_2020 <- rast(file.choose())
# Reclassificação - criar uma nova camada de floresta
floresta_2020 <- mapbiomas_2020
# Com valor de 0
values(floresta_2020) <- 0
# Atualizar categorias florestais agrupados com valor de 1
# Juntando Formação Florestal, Formação Savânica, Mangue na mesma classe = 1
floresta_2020[mapbiomas_2020==3 | mapbiomas_2020==4 | mapbiomas_2020==5] <- 1 Plotar para verificar, incluindo nomes e os cores para classes de floresta (valor = 1) e não-floresta (valor = 0).
# Passo necessario para agilizar o processamento
floresta_2020_modal<-aggregate(floresta_2020, fact=10, fun="modal")
# Mapa
tm_shape(floresta_2020_modal) +
tm_raster(style = "cat",
palette = c("0" = "#E974ED", "1" ="#129912"), legend.show = FALSE) +
tm_add_legend(type = "fill", labels = c("não-floresta", "floresta"),
col = c("#E974ED", "#129912"), title = "Classe") +
tm_shape(maza_ponto) +
tm_dots(size = 0.2, col = "yellow") +
tm_shape(maza_buffers) +
tm_borders(col = "white", lwd = 2.5) +
tm_shape(maza_buffers) +
tm_borders(col = "black", lwd = 2, lty = "dashed") +
tm_scale_bar(breaks = c(0, 20, 40), text.size = 1,
position=c("right", "bottom")) +
tm_layout(legend.bg.color="white")Se esta todo certo, voces devem ter uma imagem assim:
Figure 3.1: MapBiomas 2020 reclassificado em floresta e não-floresta. Linhas tracejadas representam buffers em torno do centro de Mazagão, com raios de 250, 500, 1000, 2000, 4000, 8000, 16000 e 32000 metros.
Vamos repetir o mesmo processo para carregar arquivos de mais anos (1985 e 2005). Baixar os arquivos link: https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_1985.tif
https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2005.tif Lembrando-se de salvar o arquivos (“utm_cover_AP_munis_maza_santana_1985.tif” e “utm_cover_AP_munis_maza_santana_2005.tif”) em um local conhecido no seu computador.
# Selecionar e carregar arquivo "utm_cover_AP_munis_maza_santana_1985.tif"
mapbiomas_1985 <- rast(file.choose())
# Reclassificação - criar uma nova camada de floresta
floresta_1985 <- mapbiomas_1985
# Com valor de 0
values(floresta_1985 <- 0
# Atualizar categorias florestais agrupados com valor de 1
# Juntando Formação Florestal, Formação Savânica, Mangue na mesma classe = 1
floresta_1985[mapbiomas_1985==3 | mapbiomas_1985==4 | mapbiomas_1985==5] <- 1
# Selecionar e carregar arquivo "utm_cover_AP_munis_maza_santana_2005.tif"
mapbiomas_2005 <- rast(file.choose())
# Reclassificação - criar uma nova camada de floresta
floresta_2005 <- mapbiomas_2005
# Com valor de 0
values(floresta_2005 <- 0
# Atualizar categorias florestais agrupados com valor de 1
# Juntando Formação Florestal, Formação Savânica, Mangue na mesma classe = 1
floresta_2005[mapbiomas_2005==3 | mapbiomas_2005==4 | mapbiomas_2005==5] <- 1
3.3 Dados: IBGE
Disonivel para cada municipio.
# IBGE
# Geral https://cidades.ibge.gov.br/brasil/ap/mazagao/panorama
# IBGE SIDRA https://sidra.ibge.gov.br/pesquisa/censo-demografico/series-temporais/series-temporais/
# Censo: Tabela 136; 1996 = Tabela 475; e 2007 = Tabela 793.
censo_pop <- c(8894, 11353, 11986, 13862, 17032)
censo_ano <- c(1991, 1996, 2000, 2007, 2010)
df_censo <- data.frame(tipo = "censo", ano = censo_ano, pop = censo_pop)
# Estimativas anuais IBGE SIDRA Tabela 6579
est_pop <- c(12410, 12633, 12933, 13139, 13913, 14259, 14418, 14655, 17420, 17794, 18739, 19157, 19571, 19981, 20387, 21206, 21632, 22053, 22468)
est_ano <- c(2001, 2002, 2003, 2004, 2005, 2006, 2008, 2009, 2011, 2012, 2013, 2014, 2015, 2016, 2017, 2018, 2019, 2020, 2021)
df_est <- data.frame(tipo = "estimativa", ano = est_ano, pop = est_pop)
# Calcular densidade demográfica
df_pop <- bind_rows(df_censo, df_est) %>%
mutate(pop_dens_km2 = pop/13294.78) %>% arrange(ano)
# Vizualizar
df_pop %>%
ggplot(aes(x=ano, y=pop_dens_km2)) +
geom_point(aes(shape=tipo), size = 4) +
stat_smooth(method = "gam") +
labs(title = "População residente em Mazagão",
y = "Densidade demográfica (hab/km²)",
x = "Ano")
#> `geom_smooth()` using formula = 'y ~ s(x, bs = "cs")'
4 Métricas da paisagem e pacote “landscapemetrics”
4.1 Calculo de métricas
Para ilustrar como rodar as funções e cálculos com landscapemetrics, vamos calcular a área central na paisagem que usamos no tutorial de Escala. Vamos estudar uma classe (floresta), portanto vamos incluir as métricas para nível de classe. Além disso, as métricas de paisagem em nível de classe são mais eficazes na definição de processos ecológicos (Tischendorf, L. Can landscape indices predict ecological processes consistently?. Landscape Ecology 16, 235–254 (2001). https://doi.org/10.1023/A:1011112719782.).
Métricas de área central (“core area”) são consideradas medidas da qualidade de hábitat, uma vez que indica quanto existe realmente de área efetiva de um fragmento, após descontar-se o efeito de borda. Vamos calcular a percentual de área central (“core area”). Isso seria, a percentual de áreas centrais (excluídas as bordas de 30 m) de cada classe em relação à área total da paisagem.
4.1.1 Região único, métrica única
Para a função sample_lsm() funcionar, precisamos informar (i) a paisagem (arquivo de raster), (ii) região de interesse (polígono), e por final (v) a métrica desejada.
minha_amostra_250 <- sample_lsm(landscape = floresta_2020,
y = maza_buffers[1, ],
plot_id = data.frame(maza_buffers)[1, 'dist_km'],
metric = "cpland",
edge_depth = 1) Depois que executar (“run”), podemos olhar os dados com o codigo a seguir.
minha_amostra_250| layer | level | class | id | metric | value | plot_id | percentage_inside |
|---|---|---|---|---|---|---|---|
| 1 | class | 0 | NA | cpland | 46.11872 | 0.25 | 99.66742 |
| 1 | class | 1 | NA | cpland | 9.13242 | 0.25 | 99.66742 |
4.1.2 Distâncias variados, métrica única
Aqui no exmplo vamos quantificar a mesma métrica para 8 distancias diferentes. Usando buffer com extensões diferentes de 250, 500, 1000, 2000, 4000, 8000, 16000 e 32000 metros distantes.
# Metricá para as distancias
minha_metrica <- sample_lsm(landscape = floresta_2020,
y = maza_buffers,
plot_id = data.frame(maza_buffers)[, 'dist_km'],
metric = "cpland",
edge_depth = 1) 5 Gráficos
Antes de apresentar os resultados, precisamos primeiramente organizar os dados. Os dados de entrada eram no objeto “maza_buffers”, formato vetor (polígonos). O formato de vetor tem uma tabela de atributos e precisamos acrescentar os valores da métrica (no objeto “minha_metrica”) associado com cada polígono (buffer) para apresentar e analisar os resultados.
Escolhendo a classe de floresta através de um filtro e selecionando as colunas desejadas e colocando-as na sequência desejada (“select”), usando a função “left_join” para adicionar as métricas aos buffers.
# Organizar dados
resultados <- maza_buffers %>%
left_join(minha_metrica %>% dplyr::filter(class==1) %>%
dplyr::select(plot_id, class, metric, value),
by=c("dist_km"="plot_id")) %>%
mutate(value = round(value,2))Agora, usando o codigo a seguir podemos visualizar os resultados em um gráfico……
# Gráfico
resultados %>%
ggplot(aes(x=dist_km, y=value)) +
geom_point() +
geom_line() +
labs(title = "Cobertura florestal ao redor do Mazagão",
subtitle = "MapBiomas ano 2020",
x = "Distância até cidade (quilômetros)",
y = "Área central de floresta\n(porcentagem da paisagem)")
Além disso, como os resultados ficaram na tabela de atributos podemos também visualizar em um mapa, neste caso um mapa interativa …….
buffer_line <- st_cast(resultados, 'LINESTRING')
#> Warning in st_cast.sf(resultados, "LINESTRING"): repeating attributes for all
#> sub-geometries for which they may not be constant
mapview::mapview(buffer_line, zcol="value", lwd = 8,
popup = popupTable(buffer_line,
zcol = c("dist_km",
"metric",
"value"))
) +
mapview::mapview(maza_ponto, col.regions= "black")Agora vamos repetir o mesmo process para varios anos e tres métricas.
# Objeto com os nomes das funções para calcular as métricas desejadas.
minhas_metricas <- c("lsm_c_cpland",
"lsm_c_enn_mn", "lsm_c_enn_sd", "lsm_c_enn_cv",
"lsm_c_pd")
# Objeto com o mesmo paisagem em anos diferentes.
floresta_anos <- c(floresta_1985, floresta_2005, floresta_2020)
# Metricás para as distancias
minha_metricas_anos <- sample_lsm(landscape = floresta_anos,
y = maza_buffers,
plot_id = data.frame(maza_buffers)[, 'dist_km'],
what = minhas_metricas,
edge_depth = 1)
#> Warning: Class 0: ENN = NA for class with only 1 patch.
#> Warning: Class 0: ENN = NA for class with only 1 patch.
#> Warning: Class 0: ENN = NA for class with only 1 patch.
#> Warning: Class 1: ENN = NA for class with only 1 patch.
# Organizar dados
resultados_anos <- maza_buffers %>%
left_join(minha_metricas_anos %>% dplyr::filter(class==1) %>%
dplyr::mutate(value = round(value,2),
ano = case_when(layer==1 ~1985,
layer==2~2005,
layer==3~2020)) %>%
dplyr::select(ano, plot_id, class, metric, value),
by=c("dist_km"="plot_id"))Grafico.
resultados_anos %>%
mutate(ext_km = (2*dist_km)) %>%
# fazer grafico
ggplot(aes(x=ano, y=value)) +
stat_smooth(linetype="dashed", colour = "magenta") +
geom_point(aes(group=ext_km, colour = factor(ext_km))) +
geom_line(aes(group=ext_km, colour = factor(ext_km))) +
scale_color_viridis_d("extensão\n(quilômetros)") +
facet_wrap(~metric, scales = "free_y") +
labs(title = "Comparação multiescala de várias métricas",
x = "ano",
y = "metric value") +
theme_bw() + theme(legend.position="top")
Gráfico alternativa.
resultados_anos %>%
mutate(ext_km = (2*dist_km)) %>%
# fazer grafico
ggplot(aes(x=ext_km, y=value)) +
stat_smooth(method="gam",
formula = y ~ s(x, bs = "cs", k=5),
linetype="dashed", colour = "magenta") +
geom_point(aes(group=ano, colour = factor(ano))) +
geom_line(aes(group=ano, colour = factor(ano))) +
scale_color_viridis_d("ano") +
facet_wrap(~metric, scales = "free_y") +
labs(title = "Comparação multiescala de várias métricas",
x = "extensão (quilômetros)",
y = "metric value") +
theme_bw() + theme(legend.position="top")