Sumário

1 Apresentaçã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 tutorial (https://rpubs.com/darren75/Impactos-das-barragens) aprenderemos sobre como analisar a cobertura da terra com métricas de paisagem em R. As tecnicas será ilustrada através cálculos usando a cobertura florestal ao redor do Rio Araguari. Ao longo do caminho, aprenderemos sobre manipulação de dados em R e aprenderemos como criar gráficos com o pacote ggplot2.

Durante a aula você aprenderá a:

  • Importar e plotar dados raster em R e mapear com os pacotes terra, sf e tmap.

  • Calcular métricas de paisagem com o pacote landscapemetrics.

  • Calcular métricas de paisagem em regiões de interesse e dentro de um buffer ao redor deles.

  • Construir gráficos com o pacote ggplot2.

Mudanças na cobertura da terra ao redor da barragem Cahoeira Caldeirao. Mostrando a Área Diretamente Afetada (linha amarela).

Figura 1.1: Mudanças na cobertura da terra ao redor da barragem Cahoeira Caldeirao. Mostrando a Área Diretamente Afetada (linha amarela).

2 Codigo e R

Objetivo não é de apresentar detalhes sobre os métodos ou os funções no R. Existem diversos exemplos disponíveis para que vocês podem obter ajudar e mais detalhes sobre manipulação, visualização de dados e análises no R, por exemplo:

  • livros online:

    • Análises Ecológicas no R desenvolvidos pelos professores durante aulas de graduação e cursos de pós-graduação nas universidades de UNESP, UFSCAR, UFRPE, UFMS, e UFRN. Capitulo 15 Dados geoespaciais tem exemplos de geoprocessamento passo a passo com o pacote mais antiga “raster”, visualização de dados geoespaciais e exemplos de aplicações de análises geoespaciais para dados ecológicos.

    • Introdução ao R curso de capacitação para servidores do Tribunal de Contas do Estado da Paraíba.

    • Ciência de Dados com R desenvolvido pelo Instituto Brasileiro de Pesquisa e Análise de Dados (IBPAD).

    • Ciência de Dados em R desenvolvido durante Cursos de Verão do Instituto de Matemática e Estatística da Universidade de São Paulo (IME-USP).

  • Tutoriais e exemplos online através buscas com Google, sempre da preferência para exemplos mais recentes porque o código e os pacotes estão sempre avançando. Exemplo buscando com “r cran introdução tutorial”….. tem mais de 50 mil resultados com exemplos, imagens e videos ……

  • Grupos de ajuda, como por exemplo:

O objetivo é de apresentar um tutorial mostrado os capacidades e opções para desenvolver e integrar pesquisas de ecologia da paisagem no ambiente estatística de R

Porque use R? R tem a capacidade (baseada em codigo) para alternar entre tarefas de processamento, modelagem e visualização de dados geográficos e não geográficos. Além disso, como é possível importar, modificar, analisar e visualizar dados espaciais no mesmo ambiente com script/código, o R permite fluxos de trabalho transparentes e reproduzíveis (A Ciência Aberta).

Aliás, atualmente a grande maioria dos artigos científicos publicados na revista Landscape Ecology incluir análises usando R.

2.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 executar

Para 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 3

Alem 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

3 Métricas da paisagem e pacote “landscapemetrics”

As métricas de paisagem são a forma que os ecólogos de paisagem usam para descrever os padrões espaciais de paisagens para depois avaliar a influência destes padrões espaciais nos padrões e processos ecológicos.

landscapemetrics tem funções para calcular métricas de paisagem em paisagems categóricos (onde tem uma classificação de cobertura de terra/habitat - modelo mancha-corredor-matriz), em um fluxo de trabalho organizado. O pacote pode ser usado como um substituto do FRAGSTATS (McGarigal et al. 1995, https://doi.org/10.2737/PNW-GTR-351), pois oferece um fluxo de trabalho reproduzível para análise de paisagem em um único ambiente (Professor McGarigal se aposentou, então FRAGSTATS não é mais apoiado). landscapemetrics também permite cálculos de quatro métricas teóricas de complexidade da paisagem: entropia marginal, entropia condicional, entropia conjunta e informação mútua (Nowosad e Stepinski 2019 https://doi.org/10.1007/s10980-019-00830-x).

Nesse pacote o formato geral para uma função é o seguinte “lsm_nível_métrica”:

  • A primeira parte é sempre lsm_ (“landscapemetric”), seguinda do “nível_” e por fim a “métrica”. Ou seja, todas as funções que calculam métricas começam com lsm_ …….

  • Daí você deve incluir o nível da análise “p” para patch (ou seja, para a mancha/fragmento), “c” para classe e “l” para landscape ou seja, métricas para a paisagem como um todo.

  • E daí existem inúmeras métricas, como por exemplo a cpland que é o percentual de área central - “core area” na paisagem, como vimos na aula teórica. Assim sendo, a função lsm_c_cpland vai calacular a métrica porcentagem da área central em cada classe. Lembrando existem metricas que podem se calculados nos trés niveis, e metricas que so pode se calculados somente para um nivel espcifco.

3.1 Pacotes

Além do “landscapemetrics”, precisamos carregar alguns pacotes a mais para facilitar a organização e apresentação de dados espaciais (vector e raster) e os resultados.

Carregar pacotes (que deve esta instalado antes):

library(landscapemetrics)
library(tidyverse)
library(sf)
library(raster)
library(terra)
library(tmap)
library(mapview)
library(kableExtra)

4 Dados

Como esperado houve mudancas na paisagem devido a barragem Cachoeira Caldeirao. Para quantificar as mudanças vamos trabalhar com os dados de MapBiomas, anos 2012 - 2020. Anos representado 4 anos antes e 4 anos depois da barragem.

Primeirmante vamos olhar o ano de 2020. Assim para verificar que todo os codigos estão validos antes de rodar o processo para todos os anos.

Existem varios formas de importar e exportar dados geoespaciais. Aqui, precisamos o arquivo com os dados de MapBiomas “utm_cover_AP_rio_2020.tif”, que voces baixaram no tutorial anterior (tutorial Escala https://rpubs.com/darren75/escala).

Link: https://github.com/darrennorris/gisdata/blob/master/inst/raster/mapbiomas_AP_utm_rio/utm_cover_AP_rio_2020.tif

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

# Selecionar e carregar arquivo "utm_cover_AP_rio_2020.tif"
mapbiomas_2020 <- rast(file.choose())
# Reclassificação - 
# Criar uma nova camada de floresta e agua 
# (floresta_2020 novo objeto de raster copiando mapbiomas_2020, 
# assim para ter os mesmos coordenados, resolução e extensão)
floresta_2020 <- mapbiomas_2020
# Todos os pixels com valor de 0
values(floresta_2020) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2020[mapbiomas_2020==3 | mapbiomas_2020==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2020[mapbiomas_2020==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2020[mapbiomas_2020==11] <- 11

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")
# Plot
tm_shape(floresta_2020_modal) +
  tm_raster(style = "cat", 
            palette = c("0" = "#E974ED", "1" ="#129912", 
                        "33" = "#0000FF", "11"="#45C2A5"), legend.show = FALSE) + 
  tm_add_legend(type = "fill", labels = c("outra", "floresta", "agua", "alagado"),
    col = c("#E974ED", "#129912", "#0000FF", "#45C2A5"), title = "Classe") + 
tm_layout(legend.bg.color = "white")

Se esta todo certo, voces devem ter uma imagem assim:

Floresta ao redor do Rio Araguari. MapBiomas 2020 reclassificado em floresta e não-floresta.

Figura 4.1: Floresta ao redor do Rio Araguari. MapBiomas 2020 reclassificado em floresta e não-floresta.

4.1 Exibir dados raster e sobreposição com uma região de interesse

Agora temos a paisagem, precisamos tambem a região de interesse (area de estudo). Por isso, precisamos carregar os dados de rios que usamos no tutorial Escala - arquivo “rivers.gpkg”.

Baixar o arquivo Link: https://github.com/darrennorris/gisdata/blob/master/inst/vector/rivers.gpkg . Lembrando-se de salvar o arquivo (“rivers.GPKG”) em um local conhecido no seu computador.

Agora, com o proximo bloco de codigo, podemos selecionar o arquivo “rivers.gpkg”, e carregar as camadas necessarios.

No exemplo, usamos %>%, que estabelece a ligação entre os passos do processo. Ou seja, %>% passa o objeto resultante automaticamente para a próxima função como primeiro argumento. Primeiramente carregamos os dados e em seguida converter (reprojeção) as coordenadas para o mesmo sistema de referência que o arquivo raster (com a função st_transform).

#  Selecionar o arquivo "rivers.GPKG",
rivers <- file.choose()
# Carregar linha central de rios, camada centerline
rio_31976 <- sf::st_read(rivers, layer = "centerline") %>% 
  st_transform(31976) 
# Carregar ponto da barragem
barragem_31976 <- sf::st_read(rivers, layer = "cachoeira_caldeirao") %>% 
  st_transform(31976) 
# Carregar Área Diretamente Afetada (ADA) poligono
ada_31976 <- sf::st_read(rivers, layer = "direct_affect") %>% 
  st_transform(31976) 
# Carregar Área Diretamente Afetada (ADA) linha
ada_31976_linha <- sf::st_read(rivers, layer = "direct_affect_line") %>% 
  st_transform(31976)

Visualizer para verificar.


# Plot
tm_shape(floresta_2020_modal) +
tm_raster(style = "cat", 
            palette = c("0" = "#E974ED", "1" ="#129912", 
                        "33" = "#0000FF", "11"="#45C2A5"), legend.show = FALSE) + 
  tm_add_legend(type = "fill", labels = c("outra", "floresta", "agua", "alagado"),
    col = c("#E974ED", "#129912", "#0000FF", "#45C2A5"), title = "Classe") +
tm_shape(rio_31976) + 
  tm_lines(col="blue") + 
tm_shape(ada_31976_linha) + 
  tm_lines(col = "yellow") + 
  tm_shape(barragem_31976) + 
  tm_dots(size=0.5, col = "yellow") + 
tm_layout(legend.bg.color="white")

Depois de executar (“run”) o código acima, você deverá ver a figura a seguir.

Cobertura da terra ao redor do Rio Araguari em 2020. Mostrando a barragem cachoeira Caldeirão (ponto amarela), Área Diretamente Afetada (linha amarela) e o rio (linha azul).

Figura 4.2: Cobertura da terra ao redor do Rio Araguari em 2020. Mostrando a barragem cachoeira Caldeirão (ponto amarela), Área Diretamente Afetada (linha amarela) e o rio (linha azul).

5 Calculo de métricas

Para ilustrar como rodar as funções e cálculos com landscapemetrics, vamos calcular a área central na paisagem. 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.).

Para calcular as métricas de paisagem dentro de um certo buffer em torno de pontos de amostra, existe a função sample_lsm(). Através da função sample_lsm() podemos calcular mais de 50 métricas da paisagem, dentro de extensões (raios/distancias) diferentes.

5.1 Região único, distâncias variados, métrica única

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”) no entorno de um ponto de amostragem. Isso seria, a percentual de áreas centrais (excluídas as bordas de 30 m) de cada classe em relação à área total da paisagem.

Vamos represntar a região de interesse com um sequencia de buffers ao redor da Área Diretamente Afetada. Por isso, precisamos carregar os dados de buffers - arquivo “barragem_impactos.GPKG”.

Baixar o arquivo Link: https://github.com/darrennorris/impactos-das-barragens/blob/main/data/vector/barragem_impactos.GPKG . Lembrando-se de salvar o arquivo (“barragem_impactos.GPKG”) em um local conhecido no seu computador.

Agora, com o proximo bloco de codigo, podemos selecionar o arquivo “barragem_impactos.GPKG”, e carregar as camada necessario.

#  Selecionar o arquivo "barragem_impactos.GPKG"
buffers <- file.choose()
# Carregar buffers ao redor da Área Diretamente Afetada (ADA)
ada_buff_31976 <- sf::st_read(buffers, layer = "faixas")

Precisamos entender exatamente o que os “buffers” representem, assim para interpretar os resltados coretamente. Visualizar com o codigo:

mapview::mapview(ada_buff_31976, zcol = "buff_dist")

Se tudo da certo, depois de executar (“run”) o código acima, você deverá ver a figura a seguir. Mostrando os cinco buffers, em faixas de distancia diferentes. Onde “0” = Área Diretamente Afetada sem buffer, “1000” = 1 a 1000 metros além do ADA, “2000” = 1000 a 2000 metros além do ADA, “4000” = 3000 a 4000 metros além do ADA e “8000” = 7000 e 8000 metros além do ADA.

mapview::mapview(ada_buff_31976, zcol = "buff_dist")

Para os calculos aqui usamos a função sample_lsm() .

Para a função sample_lsm() funcionar, precisamos informar (i) a paisagem (arquivo de raster), (ii) regiao de interesse (poligono) e (iii) a métrica desejada. Cada opção tem especificações particulares assim para que a função pode receber dados em formatos diferentes e produzir resultados conforme as necessidades de diversos estudos.

Aqui usamos os dados de 2020 para calcular uma metrica para cada buffer no obejeto (ada_buff_31976). Isso pode demorar uns minutos.

metricas_1 <- sample_lsm(landscape = floresta_2020, 
                                 y = ada_buff_31976, 
                                 plot_id = ada_buff_31976$minhaid,
                            metric = "cpland", 
                            edge_depth = 1)

Depois que executar (“run”), podemos olhar os dados com o codigo a seguir. Os dados deve ter os valores (coluna “value”) da métrica (coluna “metric”) de cada classe (coluna “class”):

metricas_1
layer level class id metric value plot_id percentage_inside
1 class 0 NA cpland 13.4064852 diferença_1000 54.70852
1 class 1 NA cpland 67.0410578 diferença_1000 54.70852
1 class 33 NA cpland 1.1871515 diferença_1000 54.70852
1 class 0 NA cpland 16.4598876 diferença_2000 25.45186
1 class 1 NA cpland 62.8098636 diferença_2000 25.45186
1 class 11 NA cpland 0.0000000 diferença_2000 25.45186
1 class 33 NA cpland 1.6539616 diferença_2000 25.45186
1 class 0 NA cpland 18.3450787 diferença_4000 19.79118
1 class 1 NA cpland 61.3108993 diferença_4000 19.79118
1 class 11 NA cpland 0.0165518 diferença_4000 19.79118
1 class 33 NA cpland 1.7493936 diferença_4000 19.79118
1 class 0 NA cpland 13.9436222 diferença_8000 12.03599
1 class 1 NA cpland 65.3349904 diferença_8000 12.03599
1 class 11 NA cpland 0.0124280 diferença_8000 12.03599
1 class 33 NA cpland 3.9481415 diferença_8000 12.03599
1 class 0 NA cpland 13.6594430 original_0 99.99878
1 class 1 NA cpland 34.7983124 original_0 99.99878
1 class 33 NA cpland 25.6864982 original_0 99.99878

5.2 Região único, distâncias variados, métricas variadas, anos variados

No exemplo anterior olhamos uma métrica da paisagem em distancias diferentes. Aqui mostraremos como incluir cálculos de diferentes métricas de paisagem ao mesmo tempo.

Primeiramente precisamos carregar dados de MapBiomas referentes os anos desejados. Neste caso vamos incluir dados de 2012 a 2020, ou seja quantro anos antes e quatro anos depois o enchimento do reservatorio em 2016.

Deve baixar os anos desjados (arqivos .tif) aqui link: https://github.com/darrennorris/impactos-das-barragens/tree/main/data/mapbiomas_AP_utm_rio

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

# limites 
myexent <- ext(vect(ada_buff_31976)) 
# arquivos
# utm_cover_AP_rio_2012.tif
rin_12 <- file.choose()
# utm_cover_AP_rio_2013.tif
rin_13 <- file.choose()
# utm_cover_AP_rio_2014.tif
rin_14 <- file.choose()
# utm_cover_AP_rio_2015.tif
rin_15 <- file.choose()
# utm_cover_AP_rio_2016.tif
rin_16 <- file.choose()
# utm_cover_AP_rio_2017.tif
rin_17 <- file.choose()
# utm_cover_AP_rio_2018.tif
rin_18 <- file.choose()
# utm_cover_AP_rio_2019.tif
rin_19 <- file.choose()
# utm_cover_AP_rio_2020.tif
rin_20 <- file.choose()

#Agora carregar. Inclundo somente a região de interesse.
mapbiomas_2012 <- rast(rin_12) %>% 
  crop(myexent, snap="out")
mapbiomas_2013 <- rast(rin_13) %>% 
  crop(myexent, snap="out")
mapbiomas_2014 <- rast(rin_14) %>% 
  crop(myexent, snap="out")
mapbiomas_2015 <- rast(rin_15) %>% 
  crop(myexent, snap="out")
mapbiomas_2016 <- rast(rin_16) %>% 
  crop(myexent, snap="out")
mapbiomas_2017 <- rast(rin_17) %>% 
  crop(myexent, snap="out")
mapbiomas_2018 <- rast(rin_18) %>% 
  crop(myexent, snap="out")
mapbiomas_2019 <- rast(rin_19) %>% 
  crop(myexent, snap="out")
mapbiomas_2020 <- rast(rin_20) %>% 
  crop(myexent, snap="out")

E agora precisamos reptir o processo de reclassificação para cada ano.

# Reclassificação - 
# Criar uma nova camada "floresta" (novo objeto de raster copiando mapbiomas_ , 
# assim para ter os mesmos coordenados, resolução e extensão). Cada ano pode incluir
# até 4 classes: 1 = floresta, 33 = agua, 11 =campo alagado, 0 = outro
floresta_2012 <- mapbiomas_2012
# Todos os pixels com valor de 0
values(floresta_2012) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2012[mapbiomas_2012==3 | mapbiomas_2012==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2012[mapbiomas_2012==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2012[mapbiomas_2012==11] <- 11 
#2013
floresta_2013 <- mapbiomas_2013
# Todos os pixels com valor de 0
values(floresta_2013) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2013[mapbiomas_2013==3 | mapbiomas_2013==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2013[mapbiomas_2013==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2013[mapbiomas_2013==11] <- 11 
#2014
floresta_2014 <- mapbiomas_2014
# Todos os pixels com valor de 0
values(floresta_2014) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2014[mapbiomas_2014==3 | mapbiomas_2014==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2014[mapbiomas_2014==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2014[mapbiomas_2014==11] <- 11 
#2015
floresta_2015 <- mapbiomas_2015
# Todos os pixels com valor de 0
values(floresta_2015) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2015[mapbiomas_2015==3 | mapbiomas_2015==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2015[mapbiomas_2015==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2015[mapbiomas_2015==11] <- 11 
#2016
floresta_2016 <- mapbiomas_2016
# Todos os pixels com valor de 0
values(floresta_2016) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2016[mapbiomas_2016==3 | mapbiomas_2016==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2016[mapbiomas_2016==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2016[mapbiomas_2016==11] <- 11 
#2017
floresta_2017 <- mapbiomas_2017
# Todos os pixels com valor de 0
values(floresta_2017) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2017[mapbiomas_2017==3 | mapbiomas_2017==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2017[mapbiomas_2017==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2017[mapbiomas_2017==11] <- 11 
#2018
floresta_2018 <- mapbiomas_2018
# Todos os pixels com valor de 0
values(floresta_2018) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2018[mapbiomas_2018==3 | mapbiomas_2018==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2018[mapbiomas_2018==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2018[mapbiomas_2018==11] <- 11 
#2019
floresta_2019 <- mapbiomas_2019
# Todos os pixels com valor de 0
values(floresta_2019) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2019[mapbiomas_2019==3 | mapbiomas_2019==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2019[mapbiomas_2019==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2019[mapbiomas_2019==11] <- 11 
#2020
floresta_2020 <- mapbiomas_2020
# Todos os pixels com valor de 0
values(floresta_2020) <- 0
# Atualizar com valor de 1 quando pixels originais são de floresta (classe 3 e 4)
floresta_2020[mapbiomas_2020==3 | mapbiomas_2020==4] <- 1 
# Atualizar com valor de 33 quando pixels originais são de agua
floresta_2020[mapbiomas_2020==33] <- 33 
# Atualizar com valor de 11 quando pixels originais são de Campo Alagado
floresta_2020[mapbiomas_2020==11] <- 11

Agora podemos reunir em uma raster unica para os proximos passos não preciso se reptir para cada ano.

floresta_12a20 <- c(floresta_2012, floresta_2013, floresta_2014, 
                     floresta_2015, floresta_2016, floresta_2017, 
                     floresta_2018, floresta_2019, floresta_2020)

Para incluir cálculos de métricas da paisagem diferentes ao mesmo tempo, precisamos acrescentar somente uma nova linha de codigo. Uma nova linha, que cria um objeto com os nomes das funções para as métricas que queremos calcular….. Também precisamos usar a opção “what” na função para aceitar os nomes das funções. Os calculos pode levar um tempo (15 minutos até uma hora dependendo o computador) para rodar devido a area (poligono de cada distancia tem > 140 km2).

# Objeto com os nomes das funções para calcular as métricas desejadas.
minhas_metricas <- c("lsm_c_pland", "lsm_c_ed", "lsm_c_cpland", 
                     "lsm_c_enn_mn", "lsm_c_enn_sd", "lsm_c_enn_cv", 
                     "lsm_c_pd","lsm_c_cohesion")

# 8 Métricas calculadas para cada distancia e cada ano.
metricas_12a20 <- sample_lsm(floresta_12a20, 
           y = ada_buff_31976,
           plot_id = ada_buff_31976$minhaid,
           what = minhas_metricas)

Depois que executar (“run”), podemos olhar os dados com o codigo a seguir.

metricas_12a20
Os dados deve ter os valores (coluna value) das métricas (coluna metric) de cada classe (coluna class) para cada distância (coluna raio):
layer level class id metric value plot_id percentage_inside
1 class 0 NA cohesion 97.27 diferença_1000 54.7
1 class 1 NA cohesion 99.43 diferença_1000 54.7
1 class 33 NA cohesion 96.97 diferença_1000 54.7
1 class 0 NA cpland 13.04 diferença_1000 54.7
1 class 1 NA cpland 69.44 diferença_1000 54.7
1 class 33 NA cpland 1.14 diferença_1000 54.7
1 class 0 NA ed 32.47 diferença_1000 54.7
1 class 1 NA ed 32.76 diferença_1000 54.7
1 class 33 NA ed 1.37 diferença_1000 54.7
1 class 0 NA enn_cv 112.98 diferença_1000 54.7

Primeiramente vamos comparar as mudancas na cobertura entre 2012 e 2020. Precisamos fazer alguns ajustes e melhorias para que os dados ficam claros nos graficos.


res_metricas <- metricas_12a20 %>% 
  left_join(data.frame(ada_buff_31976) %>% 
              dplyr::select(minhaid, tipo, buff_dist, area_km2), 
            by = c("plot_id" = "minhaid") 
            ) %>% 
   dplyr::mutate(ano = case_when(layer == 1 ~ 2012, 
                          layer == 2 ~ 2013, 
                          layer == 3 ~ 2014, 
                          layer == 4 ~ 2015, 
                          layer == 5 ~ 2016, 
                          layer == 6 ~ 2017, 
                          layer == 7 ~ 2018,
                          layer == 8 ~ 2019, 
                          layer == 9 ~ 2020
                          ), 
                 classe = case_when( 
                   class == 0 ~"outra", 
                   class == 1 ~"floresta", 
                   class == 33 ~"agua", 
                   class == 11 ~"alagado")
          )

Plotar. Aqui apresentando resultados da metrica cpland que é o percentual de área central - “core area” para cada classe, como vimos na aula teórica. Percentual de áreas centrais (excluídas as bordas de 30 m) em relação à área total da paisagem.

res_metricas %>% 
  filter(tipo == "original", metric %in% c("pland" ,"cpland")) %>% 
  mutate(nome_metrica = case_when(metric == "pland" ~"total", 
                                  metric== "cpland" ~"área central")) %>%
ggplot(aes(x=ano, y=value)) + 
  geom_vline(aes(xintercept=2016), linewidth=2, colour="yellow") +
  stat_smooth(aes(colour=classe), linetype="dashed", se=FALSE) +
  geom_point(aes(fill=classe), shape=21, size =2) + 
  scale_fill_manual(values = c("blue", "green", "grey70")) + 
  scale_colour_manual(values = c("blue", "green", "grey70")) + 
  coord_cartesian(y=c(0,60)) + 
  facet_wrap(~nome_metrica) +
  labs(x="Ano", 
       y="Cobertura (% da paisagem)")
#> `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
Mudanças de cobertura florestal na Área Diretamente Afetada (ADA).

Figura 5.1: Mudanças de cobertura florestal na Área Diretamente Afetada (ADA).

Agora, vamos fazer uma comparação enter as distancias diferentes. Para isso, precisamos fazer pequenos mudanças somente em 4 linhas do codigo anterior.


res_metricas %>% 
  filter(metric %in% c("pland" ,"cpland")) %>% 
  mutate(nome_metrica = case_when(metric == "pland" ~"total", 
                                  metric== "cpland" ~"área central")) %>%
ggplot(aes(x=ano, y=value)) + 
  geom_vline(aes(xintercept=2016), linewidth=2, colour="yellow") +
  stat_smooth(aes(colour=classe), linetype="dashed", se=FALSE) +
  geom_point(aes(fill=classe), shape=21, size =2) + 
  scale_fill_manual(values = c("blue", "aquamarine2", "#129912", "grey70")) + 
  scale_colour_manual(values = c("blue", "aquamarine2","#129912", "grey70")) + 
  facet_grid(buff_dist ~ nome_metrica, scales = "free_y") +
  labs(x="Ano", 
       y="Cobertura (% da paisagem)")


res_metricas %>% 
  filter(metric == "cpland", classe %in% c("agua", "floresta")) %>% 
  mutate(nome_metrica = case_when(metric== "cpland" ~"área central")) %>% 
  pivot_wider(id_cols = c(nome_metrica, buff_dist, ano), 
              names_from = classe, values_from = value) %>% 
  mutate(quando = case_when(ano %in% c(2012, 2013, 2014, 2015) ~"antes", 
                           ano %in% c(2017, 2018, 2019, 2020) ~"depois")) %>%
  filter(!is.na(quando)) -> antes_depois

Calcular medianos.

medianos <- antes_depois %>% 
  group_by(buff_dist, quando) %>% 
  summarise(mediano_floresta = median(floresta), 
          mediano_agua = median(agua))
#> `summarise()` has grouped output by 'buff_dist'. You can override using the
#> `.groups` argument.

Plot.

antes_depois %>%
ggplot(aes(x=agua, y=floresta)) + 
  stat_smooth(method="lm", colour="black", linetype="dashed", se=FALSE) +
  geom_hline(data = medianos, aes(yintercept = mediano_floresta, 
                                  colour=quando), 
             linetype="dashed") +
  geom_point(aes(fill=quando), shape=21, size =2) + 
  facet_wrap(~buff_dist, scales = "free") +
  labs(x="Agua (% da paisagem)", 
       y="Floresta (% da paisagem)")
#> `geom_smooth()` using formula = 'y ~ x'
Resultados apresentados incorretamente.

Figura 5.2: Resultados apresentados incorretamente.

Agora a forma de apresentação mais adequado para comunicar os resultados obtidos. Mudancas simples no eixos para facilitar interpretação mais clara e robusto.

antes_depois %>%
ggplot(aes(x=agua, y=floresta)) + 
  stat_smooth(method="lm", colour="black", linetype="dashed", se=FALSE) +
  geom_hline(data = medianos, aes(yintercept = mediano_floresta, 
                                  colour=quando), 
             linetype="dashed") +
  geom_point(aes(fill=quando), shape=21, size =2) + 
  coord_cartesian(y=c(0,70)) +
  facet_wrap(~buff_dist, scales = "free_x") +
  labs(x="Agua (% da paisagem)", 
       y="Floresta (% da paisagem)")
#> `geom_smooth()` using formula = 'y ~ x'
Comparação entre floresta e agua nos anos antes (2012-2015) e depois (2017-2020) o enchemento do reservatorio barragem Cachoeira Caldeirão. Valores anuais apresentados (pontos) são a proporção de aréa central na paisagem (%) para as classes de floresta e agua. Comparando valores na Àrea Diretamente Afetada (ADA) e em 4 distancias (metros) alem da ADA. Linhas horizontais pontilados representam os valores medianos dos anos antes e depois para cada distancia. As linhas pretas pontilados são da regressao linear entre a aréa central de agua e floresta na paisagem.

Figura 5.3: Comparação entre floresta e agua nos anos antes (2012-2015) e depois (2017-2020) o enchemento do reservatorio barragem Cachoeira Caldeirão. Valores anuais apresentados (pontos) são a proporção de aréa central na paisagem (%) para as classes de floresta e agua. Comparando valores na Àrea Diretamente Afetada (ADA) e em 4 distancias (metros) alem da ADA. Linhas horizontais pontilados representam os valores medianos dos anos antes e depois para cada distancia. As linhas pretas pontilados são da regressao linear entre a aréa central de agua e floresta na paisagem.

6 Conclusões e próximos passos

Os resultados apresentados na figura anterior não segam os resultados esperados que a cobertura de agua ia aumentar alem dos limites do ADA. Para entender melhor os padores, precisamos:

  • Aplicar testes estatisticos
  • Verficar padrões nas outras metricas calculados
  • Verficar padrões com mais anos
  • Verficar padrões com abordagens alternativas como buffers ao redor dos rios:
    • Cochrane et al. 2017. Landsat-based analysis of mega dam flooding impacts in the Amazon compared to associated environmental impact assessments: Upper Madeira River example 2006–2015. https://doi.org/10.1016/j.rsase.2017.04.005
    • Li et al. 2020. Examining Water Area Changes Accompanying Dam Construction in the Madeira River in the Brazilian Amazon https://doi.org/10.3390/w12071921