Como calcular a composição da paisagem no R?

R
ecologia de paisagem
Author

Ivana Cardoso

Published

February 26, 2024

Em ecologia de paisagens, uma paisagem pode ser definida como um mosaico de usos do solo que interagem entre si em uma determinada escala, e pode ser caracterizada por sua composição e configuração. A composição da paisagem diz respeito a quantidade dos diferentes elementos na paisagem, enquanto a configuração diz respeito ao arranjo espacial desses elementos.

Por exemplo, a cobertura florestal é uma medida de composição da paisagem que representa a quantidade de floresta, geralmente expressa em porcentagem, dentro da paisagem. Maior cobertura florestal pode indicar maior diversidade de espécies e disponibilidade de recursos (revisado por Arroyo-Rodríguez et al., 2020). Já o número de fragmentos é uma medida de configuração da paisagem que representa o nível de fragmentação dessa paisagem. Para um determinado valor de habitat, mais fragmentos podem facilitar a movimentação das espécies (revisado por Arroyo-Rodríguez et al., 2020).

A escala em que a composição e a configuração da paisagem são calculadas influencia fortemente a relação das espécies com a paisagem. Por isso, é indicado fazer uma análise multi-escalar para selecionar uma escala, chamada escala de efeito, que produza a relação mais forte entre sua variável resposta relacionada às espécies e sua variável preditora relacionada à paisagem. Por exemplo, em um estudo com aves foi utilizada uma análise multi-escalar para avaliar a relação entre número de espécies e quantidade de habitat disponível na paisagem.

No post de hoje, vamos aprender a calcular a composição da paisagem no software R utilizando o pacote landscapemetrics. Para isso, utilizaremos, para fins didáticos, poucas paisagens e apenas uma escala de 2000 m. No entanto, lembre-se de avaliar se uma análise multi-escalar seria adequada para o seu caso.

A área de estudo que escolhi foi a Reserva de Desenvolvimento Sustentável (RDS) do Rio Negro, no Amazonas. Então, baixei o raster de cobertura e uso do solo da área no MAPBIOMAS para o ano de 2022 (coleção 8). O arquivo está disponível no meu GitHub, mas caso prefira baixar diretamente do site original, sugiro utilizar o Toolkit no Google Engine e prestar atenção ao Sistema de Referência de Coordenadas (CRS) escolhido.

Calculando a composição da paisagem

# Baixe os pacotes necessarios
if (!requireNamespace(c("terra", "landscapemetrics", "sf"), quietly = TRUE)) {
  install.packages(c("terra", "landscapemetrics", "sf"))
}

# Carregue os pacotes necessarios
library(terra)
library(landscapemetrics)
library(sf)
# Defina seu diretorio
setwd("C:/Users/ivana/OneDrive/Ivana-Cardoso.github.io/posts/Composicao_paisagem") # Nao esqueca de mudar para a pasta onde voce ira salvar o raster baixado
# Importe o raster
# Você pode baixa-lo em:
# https://github.com/Ivana-Cardoso/Ivana-Cardoso.github.io/raw/main/posts/Composicao_paisagem/RDS_Rio_Negro.tif
RDS <- rast("RDS_Rio_Negro.tif")

# Veja qual o Sistema de Referencia de Coordenadas do raster
crs(RDS)
[1] "PROJCRS[\"WGS 84 / UTM zone 20S\",\n    BASEGEOGCRS[\"WGS 84\",\n        DATUM[\"World Geodetic System 1984\",\n            ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n                LENGTHUNIT[\"metre\",1]]],\n        PRIMEM[\"Greenwich\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433]],\n        ID[\"EPSG\",4326]],\n    CONVERSION[\"UTM zone 20S\",\n        METHOD[\"Transverse Mercator\",\n            ID[\"EPSG\",9807]],\n        PARAMETER[\"Latitude of natural origin\",0,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8801]],\n        PARAMETER[\"Longitude of natural origin\",-63,\n            ANGLEUNIT[\"degree\",0.0174532925199433],\n            ID[\"EPSG\",8802]],\n        PARAMETER[\"Scale factor at natural origin\",0.9996,\n            SCALEUNIT[\"unity\",1],\n            ID[\"EPSG\",8805]],\n        PARAMETER[\"False easting\",500000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8806]],\n        PARAMETER[\"False northing\",10000000,\n            LENGTHUNIT[\"metre\",1],\n            ID[\"EPSG\",8807]]],\n    CS[Cartesian,2],\n        AXIS[\"(E)\",east,\n            ORDER[1],\n            LENGTHUNIT[\"metre\",1]],\n        AXIS[\"(N)\",north,\n            ORDER[2],\n            LENGTHUNIT[\"metre\",1]],\n    USAGE[\n        SCOPE[\"Navigation and medium accuracy spatial referencing.\"],\n        AREA[\"Between 66°W and 60°W, southern hemisphere between 80°S and equator, onshore and offshore. Argentina. Bolivia. Brazil. Falkland Islands (Malvinas). Paraguay.\"],\n        BBOX[-80,-66,0,-60]],\n    ID[\"EPSG\",32720]]"

O nosso raster está no datum WGS-84, no formato de saída UTM, zona 20, hemisfério sul. É importante atentar para o formato de saída, pois queremos que as unidades do mapa sejam metros, e não graus. Isso facilita a compreensão do tamanho das nossas paisagens. É mais intuitivo entender que o tamanho das paisagens é de 2000 m do que 0.018 graus, não é?
Lembre-se de verificar a zona UTM da sua ára de estudo e garantir que seu raster esteja na projeção correta. 

# Cheque o raster
check_landscape(RDS) 
  layer       crs units   class n_classes OK
1     1 projected     m integer         9  ✔

Observe que nossas unidades estão em metros (coluna ‘units’), nosso raster possui nove classes (coluna ‘n_classes’) de uso e cobertura do solo, e que o valor dessas classes é especificado por números inteiros (coluna ‘class’). 

# Crie as coordenadas centrais das paisagens
coordenadas <- data.frame(X = c(774310, 757581, 745481),
                          Y = c(9658335, 9670747, 9688402))

# Identifique as paisagens
coordenadas$id <- c(1, 2, 3)

# Transforme a tabela de coordenadas em pontos espaciais
pontos <- st_as_sf(coordenadas, 
                        coords = c("X", "Y"), crs = crs(RDS))

Se você quiser visualizar seus dados diretamente no R antes de iniciar as análises, você pode utilizar os comandos abaixo:

# Crie as paisagens de 2000 m
paisagens <- st_buffer(pontos, dist = 2000)

# Visualize os dados
plot(RDS, col=gray.colors(34))
points(coordenadas, pch = 16)
plot(paisagens, add = TRUE, color="black")

Eu criei esse mapa no QGIS v.3.34.3, mas ele é essencialmente o mesmo do mapa gerado a partir dos comandos acima, com exceção das cores. Eu escolhi fazer no QGIS apenas para mudar as cores sem incluir e dificultar o código com ações que não são o objetivo desse tutorial. Os pontos das nossas coordenadas estão representadas no mapa pela bolinha preta com contorno branco e as paisagens de 2000 m estão representadas pelo círculo tracejado.

# Calcule a composição da paisagem
composicao <- landscapemetrics::sample_lsm(RDS, y = pontos, size = 2000,
                                           shape = "circle", plot_id = pontos$id,
                                           what = "lsm_c_pland")

Para amostrar paisagens em uma determinada escala ao redor de pontos de interesse utilizamos a função sample_lsm. Dentro dessa função especificamos nosso raster (RDS), os pontos centrais de interesse (y = pontos), a escala das nossas paisagens (size = 2000), o formato das nossas paisagens (shape = “circle”), o identificador de cada paisagem (plot_id = pontos$id) e o que gostaríamos de calcular dentro dessas paisagens (“lsm_c_pland”). A função lsm_c_pland calcula a porcentagem de cada classe do raster dentro das paisagens. No nosso caso, essa função calculou a porcentagem de cada uso do solo nas paisagens.

# Visualize os resultados
composicao
# A tibble: 16 × 8
   layer level class    id metric   value plot_id percentage_inside
   <int> <chr> <int> <int> <chr>    <dbl>   <dbl>             <dbl>
 1     1 class     3    NA pland  74.4          1              102.
 2     1 class     6    NA pland   6.75         1              102.
 3     1 class    11    NA pland   1.23         1              102.
 4     1 class    12    NA pland   0.0352       1              102.
 5     1 class    15    NA pland   0.458        1              102.
 6     1 class    33    NA pland  17.2          1              102.
 7     1 class     3    NA pland  97.8          2              102.
 8     1 class     6    NA pland   1.06         2              102.
 9     1 class    15    NA pland   0.810        2              102.
10     1 class    33    NA pland   0.374        2              102.
11     1 class     3    NA pland  68.1          3              102.
12     1 class     6    NA pland   2.26         3              102.
13     1 class    11    NA pland   2.18         3              102.
14     1 class    12    NA pland   0.507        3              102.
15     1 class    15    NA pland   9.73         3              102.
16     1 class    33    NA pland  17.2          3              102.

Essa tabela mostra nossos resultados. Preste atenção na coluna três, chamada “class”, na coluna seis, chamada “value”, e na coluna sete chamada “plot_id”.

A coluna três indica os valores dos pixels que representam as classes. Já que utilizamos o raster do MAPBIOMAS, temos que consultar o código de legenda para entender o que é cada classe do raster.

A coluna seis indica o valor, em porcentagem, de cada classe na paisagem.

A coluna sete indica a qual paisagem pertencem os valores. Com nossos resultados podemos concluir que:

A composição da paisagem 1 é:

  • 74.4% de Formação Florestal (classe 3),

  • 6.75% de Floresta Alagável (classe 6),

  • 1.23% Campo Alagado e Área Pantanosa (classe 11),

  • 0.0352% de Formação Campestre (classe 12),

  • 0.458% de Pastagem (classe 15),

  • 17.2% de Rio (classe 33).

A composição da paisagem 2 é:

  • 97.8% de Formação Florestal (classe 3),

  • 1.06% de Floresta Alagável (classe 6),

  • 0.810% de Pastagem (classe 15),

  • 0.374% de Rio (classe 33).

A composição da paisagem 3 é:

  • 68.1% de Formação Florestal (classe 3),

  • 2.26% de Floresta Alagável (classe 6),

  • 2.18% Campo Alagado e Área Pantanosa (classe 11),

  • 0.507% de Formação Campestre (classe 12),

  • 9.73% de Pastagem (classe 15),

  • 17.2% de Rio (classe 33).