Métricas de Paisaje con R

¿Se ha fragmentado la cuenca del río Pilpilco entre 2000 y 2022?

1 Introducción

1.1 Contexto

La fragmentación del paisaje puede modificar el tamaño de los hábitats, la conectividad, los efectos de borde y la diversidad biológica. En esta actividad utilizaremos datos de MapBiomas para evaluar cambios en la cuenca del río Pilpilco entre 2000 y 2022.

NoteObjetivos

• Preparar datos espaciales para análisis de paisaje.
• Comprender la importancia de los sistemas de referencia (CRS).
• Calcular métricas a nivel de parche, clase y paisaje.
• Interpretar evidencia de fragmentación.
• Comunicar resultados ecológicos.

2 La pregunta ecológica

¿Se ha fragmentado la cuenca del río Pilpilco entre 2000 y 2022?

3 Preparación

gc()

          used (Mb) gc trigger (Mb) max used (Mb)
Ncells  593801 31.8    1353351 72.3   686445 36.7
Vcells 1083323  8.3    8388608 64.0  1875920 14.4

rm(list = ls())

library(landscapemetrics)

Warning: package 'landscapemetrics' was built under R version 4.4.3

library(tidyverse)

Warning: package 'tidyverse' was built under R version 4.4.3

Warning: package 'ggplot2' was built under R version 4.4.3

Warning: package 'dplyr' was built under R version 4.4.3

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.2.1     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   4.0.3     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(terra)

Warning: package 'terra' was built under R version 4.4.3

terra 1.9.11

Attaching package: 'terra'

The following object is masked from 'package:tidyr':

    extract

library(tidyterra)

Warning: package 'tidyterra' was built under R version 4.4.3

Attaching package: 'tidyterra'

The following object is masked from 'package:stats':

    filter

library(cowplot)

Attaching package: 'cowplot'

The following object is masked from 'package:lubridate':

    stamp

4 Carga de datos

cuencas_biobio <- vect("cuencas_biobio.shp")

lc_2000 <- rast("chile_coverage_2000.tif")
lc_2022 <- rast("chile_coverage_2022.tif")

4.1 Exploración inicial

biobio_00 <- crop(lc_2000, cuencas_biobio)
biobio_22 <- crop(lc_2022, cuencas_biobio)

plot(biobio_00)

5 Un problema frecuente: CRS

Intentemos aplicar una máscara.

biobio_m <- mask(biobio_00, cuencas_biobio)

Warning: [mask] CRS do not match

¿Por qué falla?

crs(cuencas_biobio)

[1] "GEOGCRS[\"GCS_SIRGAS-Chile\",\n    DATUM[\"SIRGAS-Chile realization 1 epoch 2002\",\n        ELLIPSOID[\"GRS 1980\",6378137,298.257222101,\n            LENGTHUNIT[\"metre\",1]],\n        ID[\"EPSG\",1064]],\n    PRIMEM[\"Greenwich\",0,\n        ANGLEUNIT[\"Degree\",0.0174532925199433]],\n    CS[ellipsoidal,2],\n        AXIS[\"longitude\",east,\n            ORDER[1],\n            ANGLEUNIT[\"Degree\",0.0174532925199433]],\n        AXIS[\"latitude\",north,\n            ORDER[2],\n            ANGLEUNIT[\"Degree\",0.0174532925199433]]]"

crs(cuencas_biobio) <- crs(biobio_00)

5.1 Aplicar máscara

biobio_00_m <- mask(biobio_00, cuencas_biobio)
biobio_22_m <- mask(biobio_22, cuencas_biobio)

6 Selección de la cuenca

6.1 ¿Qué subcuencas existen?

cuencas_biobio$NOM_SUBC

 [1] "Costeras entre R. Itata y R.  Pingueral (Incl.)"                
 [2] "Costeras entre Rio Pingueral Y Rio Andalien"                    
 [3] "Isla Quiriquina"                                                
 [4] "Costeras entre Rio Andalien Y Rio Bio-Bio"                      
 [5] "Rio Andalien"                                                   
 [6] "Rio Bio-Bio Bajo"                                               
 [7] "Rio Itata Alto (Hasta Rio Diguillin)"                           
 [8] "Laja Bajo"                                                      
 [9] "Costeras entre Rio Bio-Bio y Rio Manco"                         
[10] "Rio Laja Alto (hasta bajo junta Rio Rucue)"                     
[11] "Isla Santa Maria"                                               
[12] "Costeras entre R. Manco (incl.) y R. Laraquete"                 
[13] "Rio Lia"                                                        
[14] "Costeras entre R. Laraquete (incl.) y R. Carampangue"           
[15] "Costeras entre Rio Carampangue y Punta Lavapie"                 
[16] "Costeras entre Punta Lavapie y Rio Quiapo"                      
[17] "Rio Bio-Bio entre Rio Vergara y Rio Laja"                       
[18] "R. Carampangue entre arriba R. Colorado y desembocadura"        
[19] "Rio Duqueco"                                                    
[20] "Rio Quiapo"                                                     
[21] "Rio Carampangue entre Estero Animas y Rio Colorado"             
[22] "Rio Curanilahue"                                                
[23] "Rio Carampangue hasta bajo junta Estero Animas"                 
[24] "Rio Bio-Bio entre Rio Duqueco y Rio Vergara"                    
[25] "Costeras entre Rio Quiapo y Rio Lebu"                           
[26] "Rio Bio-Bio entre Rio Ranquil y Rio Duqueco"                    
[27] "R. Lebu entre junta rios Curanilahue y Pilpilco y desembocadura"
[28] "Rio Pilpilco"                                                   
[29] "Costeras entre R. Pangue (incl.) y R. Paicavi"                  
[30] "Costeras entre Rio Lebu y Estero Pangue"                        
[31] "Rio Paicavi"                                                    
[32] "Costeras entre R. Paicavi y R. Lleullen"                        
[33] "R. Lleullen"                                                    
[34] "Costeras entre R. Lleullen y R. Tirua"                          
[35] "Rio Tirua"                                                      
[36] "Isla Mocha"                                                     

Trabajaremos con la cuenca del río Pilpilco.

pilpilco <- subset(
  cuencas_biobio,
  cuencas_biobio$NOM_SUBC == "Rio Pilpilco"
)

ventana_zoom <- ext(pilpilco)

pilpilco_00 <- crop(biobio_00, ventana_zoom)
pilpilco_22 <- crop(biobio_22, ventana_zoom)

pilpilco_00 <- mask(pilpilco_00, pilpilco)
pilpilco_22 <- mask(pilpilco_22, pilpilco)


#con pipe (más avanzado)
# pilpilco_00<-biobio_00 %>% 
#   crop(ventana_zoom) %>% 
#   mask(pilpilco)

#con pipe (más avanzado)
# pilpilco_22<-biobio_22 %>% 
#   crop(ventana_zoom) %>% 
#   mask(pilpilco)

TipPregunta

¿Cuál es la diferencia entre crop() y mask()?

7 ¿Está listo el paisaje para landscapemetrics?

check_landscape(pilpilco_00)

Warning: Caution: Coordinate reference system not metric - Units of results
based on cellsizes and/or distances may be incorrect.

  layer        crs   units   class n_classes OK
1     1 geographic degrees integer         8  ✖

check_landscape(pilpilco_22)

Warning: Caution: Coordinate reference system not metric - Units of results
based on cellsizes and/or distances may be incorrect.

  layer        crs   units   class n_classes OK
1     1 geographic degrees integer         7  ✖

Warning

Dice “No OK”. En este caso, las métricas de paisaje requieren coordenadas proyectadas (metros), no coordenadas geográficas (grados).

7.1 Reproyección a UTM 18S

target_crs <- "EPSG:32718" #UTM 18S

pil_00 <- project(pilpilco_00, target_crs, method = "near")
pil_22 <- project(pilpilco_22, target_crs, method = "near")

check_landscape(pil_00)

  layer       crs units   class n_classes OK
1     1 projected     m integer         8  ✔

check_landscape(pil_22)

  layer       crs units   class n_classes OK
1     1 projected     m integer         7  ✔

7.1.1 Discusión

¿Por qué utilizamos method = "near" y no interpolación bilineal?

8 Explorando el paisaje

8.1 Mapas con la leyenda oficial de MapBiomas

Primero, cargar la leyenda, averiguar formato y filtrar por clases presentes.

##leyenda completa, si lo necesitas:
# 
# leyenda_completa <- data.frame(
#   class = c(1, 3, 4, 9, 11, 12, 13, 21, 22, 23, 24, 25, 29, 31, 32, 33, 34, 66),
#   descripcion = c(
#     "Bosque", 
#     "Formación Forestal", 
#     "Formación Sabánica", 
#     "Silvicultura (Plantación)", 
#     "Humedal", 
#     "Pastizal", 
#     "Otras Formaciones No Forestales", 
#     "Mosaico de Agricultura y Pastura", 
#     "Área No Vegetada", 
#     "Playa, Duna y Arena", 
#     "Infraestructura Urbana", 
#     "Otras Áreas No Vegetadas", 
#     "Afloramiento Rocoso", 
#     "Cuerpo de Agua", 
#     "Glaciar y Nieve", 
#     "Cuerpo de Agua Antrópico", 
#     "Acuicultura", 
#     "Matorral"
#   ),
#   color = c(
#     "#1f8d49", "#7dc242", "#77a605", "#7a5900", "#519799", "#d6bc74", 
#     "#d89f5c", "#ffefc3", "#ffaa00", "#ffa07a", "#d4271e", "#dbd83d", 
#     "#ff8c00", "#2532e3", "#a3fbfb", "#00ffff", "#022c43", "#a89358"
#   ),
#   stringsAsFactors = FALSE
# )
# 


leyenda_completa <- read.csv("leyenda_completa.csv")

pil_00 <- as.int(pil_00)
pil_22 <- as.int(pil_22)

valores_presentes <- unique(values(pil_00))
leyenda_filtrada1 <- leyenda_completa %>%
  filter(class %in% valores_presentes)

levels(pil_00) <- leyenda_filtrada1

valores_presentes <- unique(values(pil_22))
leyenda_filtrada2 <- leyenda_completa %>%
  filter(class %in% valores_presentes)

levels(pil_22) <- leyenda_filtrada2

Después, graficar con ggplot:

# 2000
p1 <- ggplot() +
  geom_spatraster(data = pil_00, aes(fill = descripcion)) +
  scale_fill_manual(
    values = leyenda_filtrada1$color,
    breaks = leyenda_filtrada1$descripcion,
    name = "Cobertura MapBiomas",
    na.value = "transparent"
  ) +
  theme_minimal() +
  labs(
    title = "Cobertura de Suelo - Cuenca Pilpilco 2000",
    subtitle = "Datos de MapBiomas"
  )

<SpatRaster> resampled to 500888 cells.

#2022

p2 <- ggplot() +
  geom_spatraster(data = pil_22, aes(fill = descripcion)) +
  scale_fill_manual(
    values = leyenda_filtrada2$color,
    breaks = leyenda_filtrada2$descripcion,
    name = "Cobertura MapBiomas",
    na.value = "transparent"
  ) +
  theme_minimal() +
  labs(
    title = "Cobertura de Suelo - Cuenca Pilpilco 2022",
    subtitle = "Datos de MapBiomas"
  )

<SpatRaster> resampled to 500888 cells.

Visualizar individualmente, y en un solo gráfico

p1

p2

#combinar en uno:
plot_grid(p1, p2, ncol = 1, nrow = 2)

Tip

Antes de calcular métricas:

• ¿Qué coberturas parecen aumentar?
• ¿Qué coberturas parecen disminuir?
• ¿Dónde observas indicios de fragmentación?

9 Análisis a nivel de parche

Un parche corresponde a un conjunto continuo de píxeles de una misma cobertura.

parches_00 <- calculate_lsm(pil_00, level = "patch", metric = c("area", "cai")) %>% 
  mutate(year = "2000")
parches_22 <- calculate_lsm(pil_22, level = "patch", metric = c("area", "cai")) %>% 
  mutate(year = "2022")

analisis_parches <- bind_rows(parches_00, parches_22) %>%
  left_join(leyenda_completa, by = "class")

9.1 Tamaño de parches por cobertura

# Gráfico 1: distribución del tamaño de parches para todas las clases
g_parches_todo <- analisis_parches %>% 
  filter(metric == "area") %>% # Quitamos el filtro de clases para incluir todas
  ggplot(aes(x = year, y = value, fill = year)) +
  geom_boxplot(outlier.shape = 16, outlier.alpha = 0.5, width = 0.6) +
  scale_y_log10() + # escala logarítmica porque los tamaños varían de <1 ha a miles de ha
  scale_fill_manual(values = c("2000" = "#457764", "2022" = "#cccccc")) +
  facet_wrap(~descripcion, scales = "free_y") + #permite que cada panel adapte su propio eje Y
  theme_bw() +
  labs(
    title = "Análisis Área (ha) por tipo de cobertura Parches, Cuenca Pilpilco, 2000-2022",
    x = "Año",
    y = "Área del Parche (hectáreas)"
  ) +
  theme_bw()

g_parches_todo

9.1.1 Interpretación

¿Qué observaciones tienes sobre cambios en las métricas?

9.2 Formación forestal (Área)

Vamos a observar una cobertura (clase) en particular, al nivel de parche:

# Gráfico 1a: distribución del tamaño de parches para sola una clase
g_parches_forestal <- analisis_parches %>% 
  filter(metric == "area", class == 3) %>% #aqui el codigo de la leyenda, formación forestal
  ggplot(aes(x = year, y = value, fill = year)) +
  geom_boxplot(outlier.shape = 16, outlier.alpha = 0.5, width = 0.6) +
  scale_fill_manual(values = c("2000" = "#457764", "2022" = "#cccccc")) +
  scale_y_log10() + 
  theme_bw() +
  labs(
    title = "Análisis Tamaño Parches Formación Forestal, Cuenca Pilpilco, 2000-2022",
    x = "Año",
    y = "Área del Parche (hectáreas)"
  ) +
  theme_bw()

g_parches_forestal

9.3 Formación forestal, Área núcleo (CAI)

# Gráfico 1b: otra métrica
g_parches_forestal_cai <- analisis_parches %>% 
  filter(metric == "cai", class == 3)  %>%
  ggplot(aes(x = year, y = value, fill = year)) +
  geom_boxplot(outlier.shape = 16, outlier.alpha = 0.5, width = 0.6) +
  scale_fill_manual(values = c("2000" = "#457764", "2022" = "#cccccc")) +
  # scale_y_log10() + 
  theme_bw() +
  labs(
    title = "Análisis CAI Parches Formación Forestal, Cuenca Pilpilco, 2000-2022",
    x = "Año",
    y = "Área del Parche (hectáreas)"
  ) +
  theme_bw()
g_parches_forestal_cai

Note

CAI mide la proporción del parche que permanece alejada de los bordes.

10 Análisis a nivel de clase

clase_00 <- calculate_lsm(pil_00, level = "class", metric = c("ca", "np", "pland")) %>% mutate(year = "2000")
clase_22 <- calculate_lsm(pil_22, level = "class", metric = c("ca", "np", "pland")) %>% mutate(year = "2022")

analisis_clases <- bind_rows(clase_00, clase_22) %>%
  left_join(leyenda_completa, by = "class")

10.1 Dominancia del paisaje (PLAND)

g_clase_pland <- analisis_clases %>%
  filter(metric == "pland") %>%
  ggplot(aes(x = reorder(descripcion, value), y = value, fill = year)) +
  geom_bar(stat = "identity", position = "dodge") +
  coord_flip() +
  scale_fill_manual(values = c("2000" = "#457764", "2022" = "#cccccc")) +
  theme_bw() +
  labs(
    title = "Cambios en la Dominancia del Paisaje (2000 vs 2022)",
    subtitle = "Nivel de Clase: Porcentaje de la cuenca ocupado por cobertura",
    x = "Tipo de Cobertura (MapBiomas)",
    y = "Porcentaje del Paisaje (%)",
    fill = "Año"
  )
g_clase_pland

10.1.1 Pregunta

¿Qué coberturas ganaron o perdieron importancia relativa?

10.2 Evidencia de fragmentación

Un patrón clásico es:

• CA ↓
• NP ↑

g_clase_frag <- analisis_clases %>%
  filter(metric %in% c("ca", "np")) %>%
  mutate(
    year = as.factor(year),
    value = as.numeric(value),
    metric = recode(metric, ca = "Área Total (ha)", np = "N° Parches")
  ) %>%
  ggplot(aes(x = year, y = value, group = metric, color = metric, linetype = metric)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  facet_wrap(~descripcion, scales = "free_y") +
  scale_color_manual(values = c("Área Total (ha)" = "#555555", "N° Parches" = "#4a7aab")) +
  scale_linetype_manual(values = c("Área Total (ha)" = "solid", "N° Parches" = "dashed")) +
  theme_bw() +
  labs(
    title = "Cambio en Fragmentación 2000–2022",
    x = NULL, y = NULL, color = NULL, linetype = NULL
  ) +
  theme(legend.position = "bottom")
g_clase_frag

`geom_line()`: Each group consists of only one observation.
ℹ Do you need to adjust the group aesthetic?

11 Análisis a nivel de paisaje

paisaje_00 <- calculate_lsm(pil_00, level = "landscape", metric = c("shdi", "np", "pr", "contag")) %>% mutate(year = "2000")
paisaje_22 <- calculate_lsm(pil_22, level = "landscape", metric = c("shdi", "np", "pr", "contag")) %>% mutate(year = "2022")

analisis_paisaje <- bind_rows(paisaje_00, paisaje_22)

11.1 Resumen numérico

# Crear una tabla comparativa resumida 
tabla_paisaje <- analisis_paisaje %>%
  select(metric, value, year) %>%
  pivot_wider(names_from = year, values_from = value) %>%
  mutate(cambio_neto = `2022` - `2000`)
tabla_paisaje

# A tibble: 4 × 4
  metric  `2000`  `2022` cambio_neto
  <chr>    <dbl>   <dbl>       <dbl>
1 contag  75.4    77.1         1.67 
2 np     377     403          26    
3 pr       8       7          -1    
4 shdi     0.797   0.683      -0.115

11.2 Cambios en la estructura del paisaje

g_paisaje <- analisis_paisaje %>%
  filter(metric %in% c("shdi", "contag", "np", "pr")) %>%
  mutate(
    year = as.factor(year),
    value = as.numeric(value),
    metric = recode(metric, "shdi" = "Diversidad (Shannon)", "contag" = "Contagio (Agregación %)",
                    "np" = "Número de Parches", "pr" = "Riqueza de Parches")
  ) %>%                                                                          
  ggplot(aes(x = year, y = value, group = metric, color = metric)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 4) +
  facet_wrap(~metric, scales = "free_y") +
  scale_color_brewer(palette = "Set1") +
  theme_bw() +
  labs(
    title = "Análisis escala paisaje, Cuenca Pilpilco, 2000 y 2022",
    x = "Año",
    y = "Valor del Índice"
  ) +
  theme(legend.position = "none")
g_paisaje

11.2.1 Guía de interpretación

Métrica	Significado
SHDI	Diversidad del paisaje
NP	Número total de parches
PR	Riqueza de clases
CONTAG	Agregación espacial

12 Síntesis final

Completa:

Entre 2000 y 2022 la cuenca del río Pilpilco se experimentó …

Utiliza:

• una métrica de parche,
• una métrica de clase,
• una métrica de paisaje.