mazagao.Rmd

---
title: "Mazagão"
author: "Darren Norris"
date: "`r Sys.Date()`"
output:
  bookdown::html_document2:
    toc: yes
    toc_depth: 3
    toc_float: yes
    fig_caption: yes
  bookdown::pdf_document2:
    toc: yes
    toc_depth: 3
    number_sections: yes
    extra_dependencies: flafter
    highlight: tango
    includes:
      in_header: preamble.txe
urlcolor: blue
toc-title: Sumário
header-includes:
  - \counterwithin{figure}{section}
---

```{r setup, include=FALSE}
library(knitr)
knitr::opts_chunk$set(
  echo = TRUE, collapse = TRUE,
  comment = "#>" 
  )
def_hook <- knit_hooks$get("output")
knit_hooks$set(output = function(x, options) {
  out <- def_hook(x, options)
  return(paste("\\begin{framed}\\begin{verbatim}", x, "\\end{verbatim}\\end{framed}", collapse = "\n"))
})
```

\newpage{}

# 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.

\newpage

## Organização do codigo no tutorial
O tutorial está organizado em etapas de processamento, com blocos de código em caixas cinzas:
```{r, eval=FALSE}
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 `r kableExtra::text_spec("#", bold = TRUE)` é usado para incluir comentarios sobre os passos no codgio. Ou seja, linhas começando com  `r kableExtra::text_spec("#", bold = TRUE)` são ignorados por R, e não é codigo de executar.

```{r, eval=FALSE}
# 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 `r kableExtra::text_spec("#>", bold = TRUE)` e/ou `r kableExtra::text_spec("[1]", bold = TRUE)` 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).


```{r, echo=TRUE, results='asis', evaluate = TRUE, collapse = TRUE}
# Passo 1
1+1

# Passo 2
x <- 1+1
# Passo 3
x

# Passo 4
x + 1

```


\newpage


# Pacotes necessarios

```{r, message=FALSE, warning=FALSE}
library(landscapemetrics)
library(tidyverse)
library(sf)
library(raster)
library(terra)
library(tmap)
library(gridExtra)
library(kableExtra)
library(leafpop)
library(mapview)
library(mgcv)
```

# Dados

## 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](https://github.com/darrennorris/mazagao/blob/main/vector/mazagao.GPKG){target="_blank"} .
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.


```{r, echo=FALSE, eval=TRUE, message=FALSE, warning=FALSE, results='hide'}
meuSIG <- "vector/mazagao.GPKG"

# distancias ate cidade (250m, 500m, 1km, 2km, 4km, 8km, 16km e 32km)
maza_buffers <- sf::st_read(meuSIG, layer = "maza_buffers")
# lines for plotting
maza_buffer_lines <- sf::st_read(meuSIG, layer = "maza_buffers") %>% 
  st_cast('LINESTRING')
# municipio
maza_poly <- sf::st_read(meuSIG, layer = "maza_poly") 
# ponto
maza_ponto <- sf::st_read(meuSIG, layer = "maza_ponto") 
```


```{r, eval=FALSE, message=FALSE, results = FALSE}
#  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")
```

\newpage

Visualizar as camadas

```{r}
  mapview::mapview(maza_buffers, z="dist_km") + 
  mapview::mapview(maza_ponto, col.regions= "black") 

```


## 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" .

Link: [https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2020.tif](https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2020.tif){target="_blank"} 

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 <code>rast</code>.
```{r eval=FALSE, echo=TRUE, message=FALSE, warning=FALSE}
# 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 
```

```{r, echo=FALSE}
rin <- "raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2020.tif"
mapbiomas_2020 <- rast(rin)
# 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).

```{r, eval = FALSE}
# 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:

```{r, echo = FALSE, message=FALSE, warning=FALSE, fig.width=5, fig.height=5, fig.cap="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."}
# 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"), 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")
  
```


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_1985.tif){target="_blank"} 

[https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2005.tif](https://github.com/darrennorris/mazagao/blob/main/raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2005.tif){target="_blank"} 
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.

```{r eval=FALSE, echo=TRUE, message=FALSE, warning=FALSE}
# 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

```

```{r, echo=FALSE}
rin <- "raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_1985.tif"
mapbiomas_1985 <- rast(rin)
# 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 

rin <- "raster/AP_utm_munis_maza_santana/utm_cover_AP_munis_maza_santana_2005.tif"
mapbiomas_2005 <- rast(rin)
# 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 
```



## Dados: IBGE

Disonivel para cada municipio.

```{r}
# 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")
```



# Métricas da paisagem e pacote "landscapemetrics"

## 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.

###  Região único, métrica única
Para a função `r kableExtra::text_spec("sample_lsm()", background = "#dedede")` funcionar, precisamos informar 
(i) a paisagem (arquivo de raster), (ii) região de interesse (polígono), 
e por final (v) a métrica  desejada.

```{r, echo=TRUE, message=FALSE, warning=FALSE}
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. 
```{r, eval=FALSE}
minha_amostra_250
```

Os dados deve ter os valores (coluna value) da métrica (coluna metric) de cada classe (coluna class):
```{r, echo=FALSE, message=FALSE, warning=FALSE}
minha_amostra_250 %>% 
  kbl() %>%
  kable_styling(full_width = F,  latex_options = "hold_position")
```


### 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.


```{r echo=TRUE, message=FALSE, warning=FALSE}
# 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) 
```

# 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.


```{r}
# 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......
```{r}
# 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 .......


```{r}
buffer_line <- st_cast(resultados, 'LINESTRING')
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.

```{r}
# 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) 

# 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.

```{r , warning=FALSE, message=FALSE}
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.

```{r , warning=FALSE, message=FALSE}
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")

```