Sumário

Baixo Araguari

1 Apresentação

O objetivo é calcular métricas de paisagem, descrever a composição e a configuração da paisagem entre os rios Araguari e Amazonas nos municípios de Macapá e Cutias. Testando a hipótese de que a expansão do canal entre os rios foi precedida pela erosão da cobertura natural do solo.

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/baixo_araguari) aprenderemos sobre como analisar a cobertura da terra com métricas de paisagem em R.

2 Área de estudo

Antigamente o Rio Araguari deságuava no oceano Atlântico, onde o encontro das águas de rio com o oceano davam origem a Pororoca, formando as ondas mais longas do planeta utilizada para campeonatos de surfe. No entanto, atualmente, devido ao progressivo processo de expansão do Canal do Urucurituba, quase toda a água do rio é desviada para o rio Amazonas. Assim sendo, o Pororoca foi extinto, por causa do assoreamento na antiga foz do Rio Araguari.

Para visualizar um exemplo com as mudancas: https://earthengine.google.com/timelapse#v=1.07215,-50.49244,8.691,latLng&t=1.11&ps=50&bt=19840101&et=20201231&startDwell=0&endDwell=0

3 Pacotes necessarios

library(landscapemetrics)
library(tidyverse)
library(sf)
library(raster)
library(terra)
library(tmap)
library(gridExtra)
library(kableExtra)
library(leafpop)
library(mapview)
library(mgcv)
library(readxl)

4 Dados

4.1 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 1985 até 2020 “mapbiomas_1985_2020.tif” .

Link: https://github.com/darrennorris/baixo-araguari/blob/main/data/raster/mapbiomas_1985_2020.tif

Lembrando-se de salvar o arquivo (“utm_cover_AP_muni_cuitas_1985.tif”) em um local conhecido no seu computador. Agora, nós podemos carregar os dados de cobertura da terra “mapbiomas_1985_2020.tif” com a função rast.

# Selecionar e carregar arquivo "mapbiomas_1985_2020.tif"
rin <- file.choose()
mapbiomas_1985_2020 <- rast(rin)

Mapas com legenda de MapBiomas

# Aggregate para fazer a mapa
mapbiomas_1985_2020_ag <- aggregate(mapbiomas_1985_2020, fact=5, fun="modal")

# Reclassify
class_nomes <- read_excel("data/raster/AP_utm_muni_cutias/mapbiomas_6_legend.xlsx")
class_antropic <- class_nomes %>% 
  filter(type_class == "antropic") %>% pull(aid)
mapbiomas_1985_2020_ag <- classify(mapbiomas_1985_2020_ag, cbind(0, NA))
reclass_m <- as.matrix(data.frame(human = class_antropic, new_value = 18))
mapbiomas_1985_2020_ag <- classify(mapbiomas_1985_2020_ag, reclass_m)

# Make Mapbiomas legend
vals85a20 <- c(unique(values(subset(mapbiomas_1985_2020_ag, 1))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 2))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 3))),
               unique(values(subset(mapbiomas_1985_2020_ag, 4))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 5))),
               unique(values(subset(mapbiomas_1985_2020_ag, 6))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 7))),
               unique(values(subset(mapbiomas_1985_2020_ag, 8))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 9))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 10))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 11))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 12))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 13))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 14))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 15))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 16))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 17))),
               unique(values(subset(mapbiomas_1985_2020_ag, 18))),
               unique(values(subset(mapbiomas_1985_2020_ag, 19))),
               unique(values(subset(mapbiomas_1985_2020_ag, 20))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 21))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 22))),
               unique(values(subset(mapbiomas_1985_2020_ag, 23))),
               unique(values(subset(mapbiomas_1985_2020_ag, 24))),
               unique(values(subset(mapbiomas_1985_2020_ag, 25))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 26))),
               unique(values(subset(mapbiomas_1985_2020_ag, 27))),
               unique(values(subset(mapbiomas_1985_2020_ag, 28))),
               unique(values(subset(mapbiomas_1985_2020_ag, 29))),
               unique(values(subset(mapbiomas_1985_2020_ag, 30))),
               unique(values(subset(mapbiomas_1985_2020_ag, 31))),
               unique(values(subset(mapbiomas_1985_2020_ag, 32))),
               unique(values(subset(mapbiomas_1985_2020_ag, 33))),
               unique(values(subset(mapbiomas_1985_2020_ag, 34))),
               unique(values(subset(mapbiomas_1985_2020_ag, 35))), 
               unique(values(subset(mapbiomas_1985_2020_ag, 36)))
               )

class_vals <- unique(vals85a20[is.finite(vals85a20)])
class_color <- class_nomes %>% 
  dplyr::filter(aid %in% class_vals) %>% dplyr::pull(hexadecimal_code)
class_color <- paste("#", class_color, sep="")
names(class_color) <-  class_nomes %>% filter(aid %in% class_vals) %>% pull(aid)
#labelas
my_label <- class_nomes %>% 
  dplyr::filter(aid %in% class_vals) %>% dplyr::pull(classe_descricao)
my_label <- ifelse(my_label=="Agricultura", "Antropico", my_label)
names(my_label) <- class_nomes %>% filter(aid %in% class_vals) %>% pull(aid)

# Plot
tm_shape(mapbiomas_1985_2020_ag) +
  tm_raster(style = "cat", palette = class_color, labels = my_label, title="") + 
tm_scale_bar(breaks = c(0, 10), text.size = 1, 
             position=c("right", "bottom")) +
tm_facets(ncol=6) + 
tm_layout(legend.bg.color="white", legend.outside = TRUE,
            legend.outside.position = "right", 
          panel.labels = c('1985', '1986', '1987', '1988','1989', '1990',
                          '1991',  '1992', '1993', '1994', '1995','1996', 
                          '1997', '1998',
                           '1999', '2000', '2001', '2002', '2003', 
                          '2004', '2005', '2006', '2007',
                          '2008', '2009', '2010', '2011', '2012', 
                          '2013', '2014', '2015',
                           '2016', '2017', '2018', '2019', '2020'))
#> Some legend labels were too wide. These labels have been resized to 0.62. Increase legend.width (argument of tm_layout) to make the legend wider and therefore the labels larger.

4.2 Dados: Vector (rios, pontos de amostragem, buffers)

Para entender as mudanças precisamos tambem os pontos de amostra. Por isso, precisamos carregar os dados de rios e pontos de amostragem cada 0.75 km. Alem disso, poligonos representando of buffers. Vamos carregar as camadas necessarios. Baixar o arquivo Link: https://github.com/darrennorris/baixo-araguari/blob/main/data/vector/baixo_araguari.GPKG . Lembrando-se de salvar o arquivo (“baixo_araguari.GPKG”) em um local conhecido no seu computador.

Agora, com o proximo bloco de codigo, podemos selecionar o arquivo “baixo_araguari.GPKG”, e carregar camadas.

#  Selecionar o arquivo "baixo_araguari.GPKG"
meuSIG <- file.choose()

# pontos cada 0.75 km
rio_pontos <- sf::st_read(meuSIG, layer = "rio_pontos")
# linha central
rio_gurijuba <- sf::st_read(meuSIG, layer = "rio_gurijuba") 
rio_urucurituba <- sf::st_read(meuSIG, layer = "rio_urucurituba") 
# municipios
municipios_Macapa_Cutias <- sf::st_read(meuSIG, layer = "muni_poly") 
# buffers
rios_points_buffers <- sf::st_read(meuSIG, layer = "rios_points_buffers") 

# Metricás 
minhas_metricas <- c("lsm_c_cpland", 
                     "lsm_c_area_mn", "lsm_c_area_sd", "lsm_c_area_cv", 
                     "lsm_c_enn_mn", "lsm_c_enn_sd", "lsm_c_enn_cv", 
                     "lsm_c_pd", "lsm_c_cohesion")
# Buffers
mybuffers <- rios_points_buffers %>% 
  filter(buff_dist %in% c(500, 1000, 2000), 
         nome_rio == "Rio Urucurituba")
# Calcular as metricás (50 minutos mais ou menos)
# 20 anos X 9 metricas x 50 pontos X 3 buffers
# Numero de classes mude conforme ano, ponto e buffer
res_metricas <- sample_lsm(mapbiomas_1985_2020, 
                           y = mybuffers, 
                           what = minhas_metricas, 
                           plot_id = data.frame(mybuffers)[, 'aid_buff'], 
                           edge_depth=1)

# Incluir nomes para classes e ano para as camadas de raster
unique(res_metricas$layer) # 25 camadas de raster no resultados
names(mapbiomas_1985_2020) # 25 anos

res_metricas <- res_metricas %>% 
  left_join(class_nomes, by = c("class" = "aid") 
            ) %>% 
   dplyr::mutate(ano = case_when(layer == 1 ~ 1985, 
                          layer == 2 ~ 1986, 
                          layer == 3 ~ 1987, 
                          layer == 4 ~ 1988, 
                          layer == 5 ~ 1989, 
                          layer == 6 ~ 1992, 
                          layer == 7 ~ 1993, 
                          layer == 8 ~ 1994, 
                          layer == 9 ~ 1995, 
                          layer == 10 ~ 1996, 
                          layer == 11 ~ 1999, 
                          layer == 12 ~ 2000, 
                          layer == 13 ~ 2001, 
                          layer == 14 ~ 2002, 
                          layer == 15 ~ 2003, 
                          layer == 16 ~ 2008, 
                          layer == 17 ~ 2009, 
                          layer == 18 ~ 2010, 
                          layer == 19 ~ 2011, 
                          layer == 20 ~ 2012, 
                          layer == 21 ~ 2016,
                          layer == 22 ~ 2017,
                          layer == 23 ~ 2018,
                          layer == 24 ~ 2019,
                          layer == 25 ~ 2020 
                          ) 
          ) %>% 
  dplyr::select(ano, plot_id, class, classe_descricao, metric, value)

Arrumar dados.

expand.grid(ano = c(1985:2020), 
            aid_buff = unique(mybuffers$aid_buff), 
            stringsAsFactors = FALSE) %>% 
  mutate(class = 33, 
         newvalue = NA) -> df_all_points

df_all_points %>% left_join(res_metricas %>% 
                              filter(metric =="cohesion", class==33), 
                             by=c("aid_buff"="plot_id", "ano"="ano", "class"="class") 
                            ) %>% 
  left_join(data.frame(mybuffers) %>% 
              dplyr::filter(nome_rio == "Rio Urucurituba") %>% 
              dplyr::select(aid_buff, nome_rio, buff_dist, dist_amazonas_km), 
  by=c("aid_buff"="aid_buff") 
  ) %>% 
    mutate(value = replace_na(value, -1)) -> res_grafico

Grafico

res_grafico %>%
  ggplot(aes(x=ano, y = dist_amazonas_km)) + 
  geom_point(aes(colour=value)) + 
  #scale_color_gradient(low="#56B1F7", high="#132B43") + 
  scale_color_gradient2("cohesion", 
                        low = "green", mid = "yellow", high = "blue", 
                        midpoint = 50) +
  facet_wrap(~buff_dist) + 
  labs(y="distância até Rio Amazonas (km)")

Figure 4.1: Dinâmica da cobertura fluvial. Mudanças na configuração da cobertura hídrica na paisagem entre 1985 a 2020. Valores anuais de Coesão (Cohesion) de agua calculados em buffers (raios de 125, 250, 500, 1000, 2000 e 4000 metros) em torno de 50 pontos regularmente espaçados entre os rios Amazonas e Araguari.