Dinâmica de mudança de cobertura da terra e Diagrama de Sankey
Contabilização das mudanças de quantidade e área de classes de cobertura da terra ao longo do tempo e sua representação através de um Diagrama de Sankey
By Maurício Vancine in Blog
August 27, 2022
Contextualização
A coleção 7 do MapBiomas acaba de ser lançada com grandes melhorias em diversos aspectos e aumentando a janela temporal de 1985 a 2021. O MapBiomas é projeto composto por uma rede colaborativa, formada por ONGs, universidades e startups de tecnologia que produz mapeamentos anuais da cobertura e uso do solo, além de validação e elaboração de relatórios de desmatamento no Brasil.
Os mapas anuais de cobertura da terra são produzidos a partir da classificação pixel a pixel de imagens das satélites Landsat. Todo processo é feito com algoritmos de aprendizagem de máquina através da plataforma Google Earth Engine, descrito nessa figura que surrupiei do site.
Com toda certeza a parte mais interessante do MapBiomas é a analisar a variabilidade temporal. Pensando nisso, e como tive que escrever um script para ajudar o Francielson Barbosa, uma pessoa que cooriento, que está usando esses dados para entender as dinâmicas da paisagem em duas bacias do Piauí, aproveitei para escrever esse post.
Esses scrits são modificações de um trabalho que estamos submetendo do Vítor Abner Borges Dutra, onde analisamos a dinâmica da paisagem em algumas bacias hidrográficas na Amazônia Oriental no Estado do Pará. Eu achei os resultados muito interessantes que pensei que seria muito legal se outras pessoas pudessem utilizar para analisar sua própria região. Assim que ele ser publicado, eu coloco um link aqui.
Claro, vou fazer tudo no R e para Rio Claro/SP. Vou mostrar duas análises: uma tabela que mostra as mudanças das classes entre dois anos e como representar isso no formato de um Diagrama de Sankey.
Dados
Vamos iniciar instalando e carregando os pacotes.
# packages
# remotes::install_github(c("ropensci/tabulizerjars", "ropensci/tabulizer"), INSTALL_opts = "--no-multiarch")
library(tidyverse)
library(terra)
library(tabulizer)
library(tmap)
library(networkD3)
library(webshot2)
# options
options(timeout = 1e5)
Agora, façamos o download das camadas da versão 7 do MapBiomas.
# download data
download.file(url = "https://storage.googleapis.com/mapbiomas-public/brasil/collection-7/lclu/coverage/brasil_coverage_1985.tif", destfile = "/media/mude/Seagate Expansion Drive/data/mapbiomas_v07/mapbiomas_07_brasil_coverage_1985.tif")
download.file(url = "https://storage.googleapis.com/mapbiomas-public/brasil/collection-7/lclu/coverage/brasil_coverage_2021.tif", destfile = "/media/mude/Seagate Expansion Drive/data/mapbiomas_v07/mapbiomas_07_brasil_coverage_2021.tif")
Vamos também fazer o download da legenda do MapBiomas e preparar esses dados para serem usados em mapas e figuras.
# legend
legend <- tabulizer::extract_tables("https://mapbiomas-br-site.s3.amazonaws.com/C%C3%B3digos_Classes_Legenda_Cole%C3%A7%C3%A3o_7_-__PT__.docx__1_.pdf")[[1]][-(1:3), ] %>%
tibble::as_tibble() %>%
dplyr::slice(-1) |>
dplyr::mutate(legend_pt = stringr::str_trim(stringr::str_replace_all(V1, "[0-9.]", "")),
legend_en = stringr::str_trim(stringr::str_replace_all(V2, "[0-9.]", "")),
number = as.numeric(V3),
color = paste0("#", V4)) %>%
dplyr::select(legend_pt:color)
legend
## # A tibble: 35 × 4
## legend_pt legend_en number color
## <chr> <chr> <dbl> <chr>
## 1 Formação Florestal Forest Formation 3 #006400
## 2 Formação Savânica Savanna Formation 4 #32CD32
## 3 Mangue Mangrove 5 #687537
## 4 Restinga Arborizada Wooded Sandbank Vegetation 49 #6b9932
## 5 Formação Natural não Florestal Non Forest Natural Formation 10 #BBFCAC
## 6 Campo Alagado e Área Pantanosa Wetland 11 #45C2A5
## 7 Formação Campestre Grassland 12 #B8AF4F
## 8 Apicum Hypersaline Tidal Flat 32 #968c46
## 9 Afloramento Rochoso Rocky Outcrop 29 ##FF8C00
## 10 Restinga Herbácea Herbaceous Sandbank Vegetation 50 #66ffcc
## # … with 25 more rows
Vamos também fazer o download do município de Rio Claro/SP.
# rio claro municipality
rc <- geobr::read_municipality(code_muni = 3543907, year = 2020, showProgress = FALSE)
Agora sim, podemos importar os dados para o R (aqui o path do meu HD externo com os dados).
# import data
mapbiomas_1985 <- terra::rast("/media/mude/Seagate Expansion Drive/data/mapbiomas_v07/brasil_coverage_1985.tif")
mapbiomas_2021 <- terra::rast("/media/mude/Seagate Expansion Drive/data/mapbiomas_v07/brasil_coverage_2021.tif")
Fazer essas análises para o Brasil todo é extremamente demorado. Então, vamos analisar para uma área menor, apenas para Rio Claro/SP. Eu notei que na coleção 7 tem alguns pixels com zero. Então, para rodar as análises, vou substituir por NA.
# adjust to rio claro limit
mapbiomas_1985_rc <- mapbiomas_1985 %>%
terra::crop(rc) %>%
terra::mask(rc)
mapbiomas_1985_rc[mapbiomas_1985_rc == 0] <- NA
mapbiomas_2021_rc <- mapbiomas_2021 %>%
terra::crop(rc) %>%
terra::mask(rc)
mapbiomas_2021_rc[mapbiomas_2021_rc == 0] <- NA
Agora vou prepara as legendas dos mapas.
# legends
leg_1985 <- legend |>
dplyr::filter(number %in% freq(mapbiomas_1985_rc)[, 2]) |>
dplyr::arrange(number) |>
dplyr::pull(color)
leg_2021 <- legend |>
dplyr::filter(number %in% freq(mapbiomas_2021_rc)[, 2]) |>
dplyr::arrange(number) |>
dplyr::pull(color)
E por fim, plotar os dois mapas.
map_1985 <- tm_shape(mapbiomas_1985_rc) +
tm_raster(style = "cat", pal = leg_1985, title = "Legenda") +
tm_layout(main.title = "1985",
legend.position = c("left", "bottom"),
legend.title.fontface = "bold")
map_2021 <- tm_shape(mapbiomas_2021_rc) +
tm_raster(style = "cat", pal = leg_2021, title = "Legenda") +
tm_layout(main.title = "2021",
legend.position = c("left", "bottom"),
legend.title.fontface = "bold")
tmap::tmap_arrange(map_1985, map_2021)
Tabela com as mudanças das classes
Com os dados, agora podemos trabalhar nas tabelas de mudanças das classes. Aqui vou fazer três deles: quantidade de pixels, proporção e área. A classe do objeto resultante é xtable, novo para mim e meio chatinho de manipular. Devo fazer atualizações quando tiver tempo sobre essa classe.
# pixel table
table_freq <- addmargins(terra::crosstab(c(mapbiomas_1985_rc, mapbiomas_2021_rc)))
# row and column names
table_rownames_freq_class <- as.numeric(rownames(table_freq)[-length(rownames(table_freq))])
table_colnames_freq_class <- as.numeric(colnames(table_freq)[-length(colnames(table_freq))])
legend_class_rownames <- legend %>%
dplyr::filter(number %in% table_rownames_freq_class) %>%
dplyr::pull(legend_pt)
legend_class_colnames <- legend %>%
dplyr::filter(number %in% table_colnames_freq_class) %>%
dplyr::pull(legend_pt)
rownames(table_freq) <- c(legend_class_rownames, "Soma")
colnames(table_freq) <- c(legend_class_colnames, "Soma")
as.matrix(table_freq)
## brasil_coverage_2021
## brasil_coverage_1985 Formação Florestal Formação Savânica
## Formação Florestal 57854 59
## Formação Savânica 256 501
## Campo Alagado e Área Pantanosa 2023 0
## Formação Campestre 425 15
## Pastagem 0 0
## Cana 4874 709
## Outras Lavouras Temporárias 821 0
## Citrus 20190 16
## Silvicultura (monocultura) 0 0
## Mosaico de Usos 365 0
## Área Urbanizada 0 0
## Mineração 235 0
## Outras Áreas não Vegetadas 7 0
## Rio, Lago e Oceano 0 0
## Soma 87050 1300
## brasil_coverage_2021
## brasil_coverage_1985 Campo Alagado e Área Pantanosa
## Formação Florestal 1926
## Formação Savânica 0
## Campo Alagado e Área Pantanosa 1526
## Formação Campestre 0
## Pastagem 0
## Cana 959
## Outras Lavouras Temporárias 244
## Citrus 677
## Silvicultura (monocultura) 0
## Mosaico de Usos 0
## Área Urbanizada 0
## Mineração 0
## Outras Áreas não Vegetadas 0
## Rio, Lago e Oceano 0
## Soma 5332
## brasil_coverage_2021
## brasil_coverage_1985 Formação Campestre Pastagem Soja Cana
## Formação Florestal 98 0 3367 5259
## Formação Savânica 56 0 403 0
## Campo Alagado e Área Pantanosa 0 0 100 124
## Formação Campestre 839 0 98 7
## Pastagem 0 59 12 1
## Cana 234 58 55204 31798
## Outras Lavouras Temporárias 0 0 3750 56077
## Citrus 214 8 13327 50273
## Silvicultura (monocultura) 0 0 0 0
## Mosaico de Usos 0 0 456 868
## Área Urbanizada 0 0 0 0
## Mineração 31 0 16 0
## Outras Áreas não Vegetadas 0 0 43 567
## Rio, Lago e Oceano 0 0 0 109
## Soma 1472 125 76776 145083
## brasil_coverage_2021
## brasil_coverage_1985 Outras Lavouras Temporárias Café Citrus
## Formação Florestal 15052 691 267
## Formação Savânica 321 0 0
## Campo Alagado e Área Pantanosa 306 0 0
## Formação Campestre 666 0 0
## Pastagem 2 0 0
## Cana 76611 9856 1762
## Outras Lavouras Temporárias 14150 1705 1001
## Citrus 91017 20220 2590
## Silvicultura (monocultura) 2 20427 0
## Mosaico de Usos 2501 3858 459
## Área Urbanizada 0 41 0
## Mineração 888 47 23
## Outras Áreas não Vegetadas 126 0 1
## Rio, Lago e Oceano 119 0 0
## Soma 201761 56845 6103
## brasil_coverage_2021
## brasil_coverage_1985 Outras Lavouras Perenes
## Formação Florestal 42
## Formação Savânica 0
## Campo Alagado e Área Pantanosa 0
## Formação Campestre 0
## Pastagem 0
## Cana 372
## Outras Lavouras Temporárias 1634
## Citrus 1969
## Silvicultura (monocultura) 0
## Mosaico de Usos 207
## Área Urbanizada 2695
## Mineração 69
## Outras Áreas não Vegetadas 1
## Rio, Lago e Oceano 0
## Soma 6989
## brasil_coverage_2021
## brasil_coverage_1985 Silvicultura (monocultura) Mosaico de Usos
## Formação Florestal 11 128
## Formação Savânica 7 0
## Campo Alagado e Área Pantanosa 0 0
## Formação Campestre 9 0
## Pastagem 0 0
## Cana 74 642
## Outras Lavouras Temporárias 45 355
## Citrus 337 2038
## Silvicultura (monocultura) 0 0
## Mosaico de Usos 50 2
## Área Urbanizada 0 0
## Mineração 628 0
## Outras Áreas não Vegetadas 0 0
## Rio, Lago e Oceano 0 0
## Soma 1161 3165
## brasil_coverage_2021
## brasil_coverage_1985 Área Urbanizada Mineração
## Formação Florestal 18 0
## Formação Savânica 0 0
## Campo Alagado e Área Pantanosa 0 0
## Formação Campestre 0 0
## Pastagem 0 0
## Cana 261 16
## Outras Lavouras Temporárias 255 0
## Citrus 374 14
## Silvicultura (monocultura) 0 0
## Mosaico de Usos 0 0
## Área Urbanizada 0 0
## Mineração 0 0
## Outras Áreas não Vegetadas 0 0
## Rio, Lago e Oceano 0 0
## Soma 908 30
## brasil_coverage_2021
## brasil_coverage_1985 Outras Áreas não Vegetadas Rio, Lago e Oceano
## Formação Florestal 1437 80
## Formação Savânica 0 0
## Campo Alagado e Área Pantanosa 126 2
## Formação Campestre 0 0
## Pastagem 0 0
## Cana 1737 194
## Outras Lavouras Temporárias 1810 1
## Citrus 5083 46
## Silvicultura (monocultura) 0 0
## Mosaico de Usos 4 0
## Área Urbanizada 0 0
## Mineração 0 0
## Outras Áreas não Vegetadas 53 0
## Rio, Lago e Oceano 188 0
## Soma 10438 323
## brasil_coverage_2021
## brasil_coverage_1985 Soma
## Formação Florestal 86289
## Formação Savânica 1544
## Campo Alagado e Área Pantanosa 4207
## Formação Campestre 2059
## Pastagem 74
## Cana 185361
## Outras Lavouras Temporárias 81848
## Citrus 208393
## Silvicultura (monocultura) 20429
## Mosaico de Usos 8770
## Área Urbanizada 2736
## Mineração 1937
## Outras Áreas não Vegetadas 798
## Rio, Lago e Oceano 416
## Soma 604861
Agora com essa tabela pronta, podemos calcular a proporção.
# proportion table
table_freq_per <- round(table_freq/sum(freq(mapbiomas_1985_rc)[, 3]), 4) * 100
rownames(table_freq_per) <- c(legend_class_rownames, "Soma")
colnames(table_freq_per) <- c(legend_class_colnames, "Soma")
E fazendo uma pequena gambiarra, calcular a área em hectares.
# area table
area_pixel <- (res(mapbiomas_2021_rc)[1] * 3600 * 30^2)/10000
table_freq_area <- round(table_freq * area_pixel, 1)
rownames(table_freq_area) <- c(legend_class_rownames, "Soma")
colnames(table_freq_area) <- c(legend_class_colnames, "Soma")
Agora podemos exportar essas tabelas.
# export
write.table(table_freq, "content/blog/geo-cob-terra-sankey/table_freq.csv", quote = FALSE, row.names = TRUE, sep = ";", fileEncoding = "UTF-16LE")
write.table(table_freq_per, "content/blog/geo-cob-terra-sankey/table_freq_per.csv", quote = FALSE, row.names = TRUE, sep = ";", fileEncoding = "UTF-16LE")
write.table(table_freq_area, "content/blog/geo-cob-terra-sankey/table_freq_area.csv", quote = FALSE, row.names = TRUE, sep = ";", fileEncoding = "UTF-16LE")
Diagrama de Sankey
Rapaz, te falar, esse não é um código bonito, mas usa uma integração com JavaScript, então me dá um desconto. E mais importante, ele funciona (mais ou menos, tem um erro na parte das cores de cada classe que não consigui resolver, aceito ajuda…). Não vou explicar os detalhes, ainda estou processando…
# sankey plot
# freq
freq <- terra::crosstab(c(mapbiomas_1985_rc, mapbiomas_2021_rc))
# links
link_info <- freq %>%
tibble::as_tibble() %>%
dplyr::mutate(mapbiomas_1985 = as.numeric(brasil_coverage_1985),
mapbiomas_2021 = as.numeric(brasil_coverage_2021)) |>
dplyr::arrange(mapbiomas_1985, mapbiomas_2021) |>
dplyr::select(mapbiomas_1985, mapbiomas_2021, n) |>
dplyr::left_join(legend[, c(3, 1)], by = c("mapbiomas_1985" = "number")) |>
dplyr::rename(mapbiomas_1985_leg = legend_pt) |>
dplyr::left_join(legend[, c(3, 1)], by = c("mapbiomas_2021" = "number")) |>
dplyr::rename(mapbiomas_2021_leg = legend_pt) |>
dplyr::mutate(mapbiomas_1985_leg = paste("1985", mapbiomas_1985_leg),
mapbiomas_2021_leg = paste("2021", mapbiomas_2021_leg))
# nodes
nodes_info <- tibble::tibble(node_name = c(as.character(link_info$mapbiomas_1985_leg),
as.character(link_info$mapbiomas_2021_leg)) %>% unique(),
node_id = 0:(length(unique(c(as.character(link_info$mapbiomas_1985_leg),
as.character(link_info$mapbiomas_2021_leg))))-1),
node_column = c(rep(1, length(unique(as.character(link_info$mapbiomas_1985_leg)))),
rep(2, length(unique(as.character(link_info$mapbiomas_2021_leg))))),
node_group = as.character(c(1:length(unique(as.character(link_info$mapbiomas_1985_leg))),
1:length(unique(as.character(link_info$mapbiomas_2021_leg))))))
# link and nodes
link_info_node_1985 <- link_info |>
dplyr::left_join(nodes_info[, c("node_name", "node_id")], by = c("mapbiomas_1985_leg" = "node_name")) |>
dplyr::rename(source = node_id)
link_info_node_2021 <- link_info |>
dplyr::left_join(nodes_info[, c("node_name", "node_id")], by = c("mapbiomas_2021_leg" = "node_name")) |>
dplyr::rename(target = node_id)
link_info_node <- dplyr::bind_cols(link_info_node_1985, link_info_node_2021[, "target"]) |>
dplyr::select(source, target, n) |>
dplyr::mutate(source = as.numeric(source),
target = as.numeric(target),
n = as.integer(n)) |>
dplyr::arrange(source)
# color scale
cover <- nodes_info |>
dplyr::mutate(node_name = stringr::str_replace(node_name, "1985", ""),
node_name = stringr::str_replace(node_name, "2021", ""),
node_name = stringr::str_trim(node_name)) |>
dplyr::select(node_name)
col <- dplyr::left_join(cover, legend[, c("legend_pt", "color")], by = c("node_name" = "legend_pt")) |>
dplyr::pull(color)
color_scale <- paste0('d3.scaleOrdinal().range([', paste0('"', paste(col, collapse = '", "'),'"'),'])')
# gambiarra
nodes_info$node_group <- col
# sankey plot
sn <- networkD3::sankeyNetwork(Links = link_info_node,
Nodes = nodes_info,
Source = "source",
Target = "target",
Value = "n",
NodeID = "node_name",
NodeGroup = "node_group",
fontSize = 12,
nodeWidth = 25,
colourScale = color_scale)
sn
Esse gráfico vai aparecer na aba Viewer do seu RStudio. Para exportar como algo interativo que você pode abrir no seu brownse de internet ou colocar num site ou usar no RMarkdown ou quanto, exportamos em .html. Caso queira algo fixo, o código mostra como exportar em .png, mas acho que isso pode requerer algumas configurações a mais que não vou entrar em detalhes.
Outra possibilidade é simplesmente usar a opção de Export da aba de Viewer RStudio.
# export .html
networkD3::saveNetwork(sn, "content/blog/geo-sankey/sn.html")
# export .png
webshot2::webshot(url = "/home/mude/data/github/mauriciovancine/apero/content/blog/geo-sankey/sn.html",
file = "/home/mude/data/github/mauriciovancine/apero/content/blog/geo-sankey/sn.png",
vwidth = 1000, vheight = 900)
Fonte da imagem: Zbyněk Skrčený