SIMOCUTE - Ejercicio metodológico de interpretación de puntos de muestreo

Resumen

FIESTA es una biblioteca de software para análisis de datos de inventarios forestales. En este documento se explica el uso de FIESTA mediante un caso de ejemplo basado en un conjunto de datos interpretados por funcionarios e investigadores de instituciones que forman parte de la Mesa de Monitoreo por Puntos de SIMOCUTE.

Introducción

FIESTA (Forest Inventory Estimation and Analysis) es una biblioteca de software, desarrollada en el lenguaje de programación R, para el análisis de datos de inventarios forestales basados en muestras. Fue desarrollada por el Programa de Inventario y Análisis Forestal (FIA) del Servicio Forestal del Departamento de Agricultura (USDA) de los Estados Unidos.

En este documento se detalla el análisis, mediante FIESTA, de un conjunto de datos de 26125 puntos de muestreo de 1045 parcelas de monitoreo, ubicadas en la zona norte de Costa Rica. Estos datos provienen de un ejercicio de interpretación de puntos de muestreo, realizado en el marco del Sistema Nacional de Monitoreo del Uso, Cobertura y Ecosistemas de Costa Rica (SIMOCUTE), en el cual participaron 26 investigadores y funcionarios de diferentes instituciones y organizaciones que forman parte de la Mesa de Monitoreo por Puntos. La interpretación se realizó con la plataforma Collect Earth Online, un sistema colaborativo para la interpretación de imágenes satelitales.

Para cada punto se interpretaron las variables:

• Cobertura de la tierra.
• Uso del suelo.

La interpretación se realizó en dos tiempos:

• t1: entre 2005 y 2007.
• t2: en 2019.

El documento fue desarrollado en el sistema de publicación técnica y científica Quarto, el cual combina código en R y sus salidas (tablas, gráficos, mapas) con texto en Markdown.

Trabajo previo al análisis

Instalación de software

Sistema base de R y herramientas de desarrollo

Para trabajar con FIESTA, se debe instalar:

• El sistema base de R. R es un lenguaje de programación enfocado en análisis estadístico y visualización de datos.
• La interfaz de desarrollo integrada RStudio Desktop, la cual proporciona un editor de texto y otras herramientas para escribir programas en R y visualizar sus resultados, entre otras facilidades. También pueden emplearse otros ambientes de desarrollo que ofrecen características similares (ej. Visual Studio Code).

Versión mínima de R

De acuerdo con su documentación, la versión mínima de R que requiere FIESTA, a la fecha de escritura de este documento (2023-12-02), es la 4.2.0. Puede consultar la versión de su instalación al ejecutar el siguiente comando en la consola de R:

# Versión de R
R.version.string

[1] "R version 4.3.2 (2023-10-31 ucrt)"

La salida del comando anterior debe indicar que la versión de su instalación de R es mayor o igual a 4.2.0.

Paquetes de R

FIESTA

El paquete FIESTA está disponible en CRAN (Comprehensive R Archive Network), un repositorio en línea que alberga una amplia colección de paquetes y extensiones para el lenguaje de programación R, lo que facilita su instalación y actualización a nuevas versiones.

Para usar FIESTA, debe instalarse primero. Puede utilizar la función install.packages().

# Instalación del paquete FIESTA
install.packages("FIESTA")

El resultado de la instalación puede verificarse al cargar el paquete con la función library().

# Carga del paquete FIESTA
library(FIESTA)

Si el comando anterior no genera ningún mensaje de error, FIESTA debe haberse instalado adecuadamente.

Otros

Además de FIESTA, se recomienda instalar los siguientes paquetes para procesamiento y visualización de datos.

# Paquete para el desarrollo de documentos computacionales
install.packages("rmarkdown")

# Colección de paquetes para análisis de datos
install.packages("tidyverse")

# Estilos para gráficos de tidyverse
install.packages("ggthemes")

# Paquete para limpieza de datos
install.packages("janitor")

# Paquete para tablas interactivas
install.packages("DT")

# Paquete para graficación interactiva
install.packages("plotly")

# Paquete para mapas interactivos
install.packages("leaflet")

# Funciones adicionales para leaflet
install.packages("leaflet.extras")

# Funciones adicionales para leaflet
install.packages("leafem")

Luego de instalarlos, debe cargar los paquetes con la función library().

# Carga de paquetes adicionales
library(rmarkdown)
library(tidyverse)
library(ggthemes)
library(janitor)
library(DT)
library(plotly)
library(leaflet)
library(leaflet.extras)
library(leafem)
library(sf) # se instala con FIESTA

Obtención de este repositorio

Este documento, junto con otros similares, forma parte de un repositorio en GitHub, una plataforma en línea para compartir código fuente de aplicaciones, basada en el sistema de control de versiones Git. El repositorio contiene el código fuente del documento y los datos que se utilizan en los ejemplos. Su dirección es https://github.com/mesa-monitoreo-puntos/fiesta.

Puede descargar el repositorio, como un archivo ZIP, de https://github.com/mesa-monitoreo-puntos/fiesta/archive/refs/heads/main.zip

También puede “clonar” el repositorio mediante el comando clone de Git:

# Clonación de este repositorio
git clone https://github.com/mesa-monitoreo-puntos/fiesta.git

Una vez que el repositorio haya sido descargado o clonado, puede abrirse con RStudio o con otra herramienta de desarrollo.

Variables generales

En esta sección se definen algunas variables generales del proceso, incluyendo:

• La ruta al archivo con los datos de puntos de muestreo.
• Los códigos de los colores de las clases de cobertura de la tierra y uso del suelo.
• El área de la zona de estudio.

# Ruta a los datos de puntos de muestreo
ARCHIVO_PUNTOS <- "../../datos/Resul_Fin_2Grupos_SinReplicas_csv_Fix.csv"

# Colores
COLOR_VEGETACION         = "#00a600" # Corine EU - Coniferous forest
COLOR_SIN_VEGETACION     = "#dcdcc8" # Corine CR - Tierras desnudas o degradadas
COLOR_AGUA               = "#00ccf2" # Corine EU - Water courses
COLOR_NUBES_Y_SOMBRAS    = "#800080" # Púrpura
COLOR_SIN_INFORMACION    = "#000000" # Negro

COLOR_MC_BOSQUES         = "#266900" # Corine CR - Bosque denso
COLOR_AGRICULTURA        = "#Becd05" # Corine CR - Mosaico de cultivos
COLOR_GANADERIA_Y_PASTOS = "#ffffa6" # Corine CR - Pastos limpios
COLOR_ZONAS_HUMEDALES    = "#00ccf2" # Corine EU - Water courses
COLOR_NO_CLASIFICACION   = "#000000" # Negro
COLOR_INFRAESTRUCTURA    = "#f6d9df" # Corine CR - Zona urbana continua
COLOR_OTRAS_TIERRAS      = "#005acf" # Corine CR - Canales

COLOR_NULO               = "#808080" # Gris

# Paleta de colores de tipos de cobertura de la tierra
COLORES_COBERTURA <- 
    c(
        "1000-Vegetacion"         = COLOR_VEGETACION,
        "T1-1000-Vegetacion"      = COLOR_VEGETACION,
        "T2-1000-Vegetacion"      = COLOR_VEGETACION,
        "2000-Sin vegetacion"     = COLOR_SIN_VEGETACION, 
        "T1-2000-Sin vegetacion"  = COLOR_SIN_VEGETACION,
        "T2-2000-Sin vegetacion"  = COLOR_SIN_VEGETACION,
        "3000-Agua"               = COLOR_AGUA, 
        "T1-3000-Agua"            = COLOR_AGUA,
        "T2-3000-Agua"            = COLOR_AGUA,
        "4000-Nubes y sombras"    = COLOR_NUBES_Y_SOMBRAS,
        "T1-4000-Nubes y sombras" = COLOR_NUBES_Y_SOMBRAS,
        "T2-4000-Nubes y sombras" = COLOR_NUBES_Y_SOMBRAS,
        "5000-Sin informacion"    = COLOR_SIN_INFORMACION,
        "T1-5000-Sin informacion" = COLOR_SIN_INFORMACION,
        "T2-5000-Sin informacion" = COLOR_SIN_INFORMACION,
        "Nulo"                    = COLOR_NULO,
        "T1-Nulo"                 = COLOR_NULO,
        "T2-Nulo"                 = COLOR_NULO
    )

# Paleta de colores de tipos de uso de la tierra
COLORES_USO <- 
    c(
        "1000-MC bosques"       = COLOR_MC_BOSQUES,
        "T1-1000-MC bosques"    = COLOR_MC_BOSQUES,
        "T2-1000-MC bosques"    = COLOR_MC_BOSQUES,
        "2000-Agricultura"      = COLOR_AGRICULTURA,
        "T1-2000-Agricultura"   = COLOR_AGRICULTURA,
        "T2-2000-Agricultura"   = COLOR_AGRICULTURA,
        "3000-Ganad y past"     = COLOR_GANADERIA_Y_PASTOS,
        "T1-3000-Ganad y past"  = COLOR_GANADERIA_Y_PASTOS,
        "T2-3000-Ganad y past"  = COLOR_GANADERIA_Y_PASTOS,
        "4000-Zonas humed"      = COLOR_ZONAS_HUMEDALES,
        "T1-4000-Zonas humed"   = COLOR_ZONAS_HUMEDALES,
        "T2-4000-Zonas humed"   = COLOR_ZONAS_HUMEDALES,
        "5000-Infraest"         = COLOR_INFRAESTRUCTURA,
        "T1-5000-Infraest"      = COLOR_INFRAESTRUCTURA,
        "T2-5000-Infraest"      = COLOR_INFRAESTRUCTURA,
        "6000-Otras tierras"    = COLOR_OTRAS_TIERRAS,
        "T1-6000-Otras tierras" = COLOR_OTRAS_TIERRAS,
        "T2-6000-Otras tierras" = COLOR_OTRAS_TIERRAS,
        "7000 No clasif"        = COLOR_NO_CLASIFICACION,
        "T1-7000 No clasif"     = COLOR_NO_CLASIFICACION,
        "T2-7000 No clasif"     = COLOR_NO_CLASIFICACION,
        "Nulo"                  = COLOR_NULO,
        "T1-Nulo"               = COLOR_NULO,
        "T2-Nulo"               = COLOR_NULO
    )

# Área de estudio en hectáreas
AREA_ESTUDIO = 318167.67

Carga y limpieza de datos

Los datos de los puntos de muestreo se proporcionaron en un archivo CSV. En el siguiente bloque de código se cargan en el dataframe puntos. Los nombres de las columnas se “limpian” para evitar la presencia de números al inicio y otros problemas que dificultan su manejo. Los valores nulos (NA) se convierten a hileras de texto.

# Carga de datos de puntos de muestreo
puntos <- read_delim(ARCHIVO_PUNTOS)

# Limpieza de los nombres de columnas
puntos <- clean_names(puntos)

# Reemplazo de valores NA por la hilera "Nulo"
puntos <- 
    puntos |> 
    mutate(
        t1_cobertura = replace_na(t1_cobertura, "Nulo"),
        t2_cobertura = replace_na(t2_cobertura, "Nulo"),
        t1_uso = replace_na(t1_uso, "Nulo"),
        t2_uso = replace_na(t2_uso, "Nulo")
    )

# Coversión a de t1_cobertura, t2_cobertura, t1_uso y t2_uso a factores
puntos <-
    puntos |>
    mutate(
        t1_cobertura = factor(t1_cobertura, levels = unique(t1_cobertura)),
        t2_cobertura = factor(t2_cobertura, levels = unique(t2_cobertura)),
        t1_uso = factor(t1_uso, levels = unique(t1_uso)),
        t2_uso = factor(t2_uso, levels = unique(t2_uso))        
    )

El dataframe puntos contiene más de 40 columnas, provenientes del archivo CSV. Algunas de las más importantes para efectos de este análisis son:

• plot_id: identificador de la parcela.
• sample_id: identificador del punto de muestreo.
• t1_cobertura: interpretación de la variable de cobertura de la tierra en t1.
• t2_cobertura: interpretación de la variable de cobertura de la tierra en t2.
• t1_uso: interpretación de la variable de uso del suelo en t1.
• t2_uso: interpretación de la variable de uso del suelo en t2.

La siguiente tabla muestra las columnas mencionadas.

# Despliegue de los datos de puntos de muestreo en una tabla
puntos |>
    select(plot_id, sample_id, t1_cobertura, t2_cobertura, t1_uso, t2_uso) |>
    datatable(
        caption = "Puntos de muestreo de cobertura y uso de la tierra",
        rownames = FALSE,
        colnames = c(
            "plot_id", "sample_id", 
            "t1_cobertura", "t2_cobertura",
            "t1_uso", "t2_uso"
        ),
        options = list(
            pageLength = 5,
            language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
        )
    )

Análisis con paquetes básicos de R

En esta sección, se analizan los datos con paquetes de R como los de la colección de Tidyverse, Plotly R para gráficos estadísticos interactivos y Leaflet R para mapas interactivos.

Primero se analiza la variable de cobertura de la tierra y posteriormente la variable de uso del suelo.

Cobertura de la tierra

Mapa

El siguiente mapa interactivo, implementado con el paquete Leaflet, muestra la distribución espacial de la interpretación de la variable de cobertura de la tierra en t1 y t2.

# Conversión de datos de puntos de muestreo a objeto sf (vectorial de puntos)
geo_puntos <-
    puntos |>
    select(plot_id, sample_id, lon, lat, t1_cobertura, t2_cobertura) |>
    st_as_sf(
        coords = c("lon", "lat"),
        crs = 4326
  )

colores <- 
    colorFactor(
        palette = c(
            COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
            COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO
        ), 
    levels = c(
        "1000-Vegetacion", "2000-Sin vegetacion", "3000-Agua",
        "4000-Nubes y sombras", "5000-Sin informacion", "Nulo"
    )
)

leaflet() |>
    addTiles(group = "OSM") |>
    addProviderTiles(
        provider = providers$Esri.WorldImagery, 
        group = "ESRI World Imagery"
    ) |>
    addProviderTiles(
        provider = providers$CartoDB.DarkMatter,
        group = "Dark Matter"
    ) |>        
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t1_cobertura),
        color = ~colores(geo_puntos$t1_cobertura),      
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Cobertura en t1: </strong>", geo_puntos$t1_cobertura),
      paste0("<strong>Cobertura en t2: </strong>", geo_puntos$t2_cobertura),
      sep = '<br/>'
    ),      
        group = "Cobertura en t1"
    ) |>
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t1_cobertura,
    title = "Cobertura en t1",
    group = "Cobertura en t1"    
  ) |>      
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t2_cobertura),
        color = ~colores(geo_puntos$t2_cobertura),      
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Cobertura en t1: </strong>", geo_puntos$t1_cobertura),
      paste0("<strong>Cobertura en t2: </strong>", geo_puntos$t2_cobertura),
      sep = '<br/>'
    ),      
        group = "Cobertura en t2"
    ) |>    
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t2_cobertura,
    title = "Cobertura en t2",
    group = "Cobertura en t2"    
  ) |>  
    addLayersControl(
        baseGroups = c("OSM", "ESRI World Imagery", "Dark Matter"),
        overlayGroups = c("Cobertura en t1", "Cobertura en t2"),
    )

Distribución de puntos de muestreo en clases de cobertura

La distribución de los puntos de muestreo en clases de cobertura de la tierra, en t1 y t2, se muestra en un gráfico de barras y en un gráfico de pastel, ambos elaborados con el paquete Plotly.

t1

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t1_cobertura), fill = fct_infreq(t1_cobertura))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_COBERTURA) +
    xlab("Cobertura") +
    ylab("Cantidad de puntos") +
    labs(fill = "Cobertura") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

# Total de puntos
conteo <- puntos |>
    count(t1_cobertura)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_COBERTURA[as.character(conteo$t1_cobertura)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t1_cobertura, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

t2

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t2_cobertura), fill = fct_infreq(t2_cobertura))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_COBERTURA) +
    xlab("Cobertura") +
    ylab("Cantidad de puntos") +
    labs(fill = "Cobertura") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

# Total de puntos
conteo <- puntos |>
    count(t2_cobertura)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_COBERTURA[as.character(conteo$t2_cobertura)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t2_cobertura, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Uso del suelo

Mapa

El siguiente mapa interactivo, implementado con el paquete Leaflet, muestra la distribución espacial de la interpretación de la variable de uso del suelo en t1 y t2.

# Conversión de datos de puntos de muestreo a objeto sf (vectorial de puntos)
geo_puntos <-
    puntos |>
    select(plot_id, sample_id, lon, lat, t1_uso, t2_uso) |>
    st_as_sf(
        coords = c("lon", "lat"),
        crs = 4326
  )

colores <- 
    colorFactor(
        palette = c(
            COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
            COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS, 
            COLOR_NO_CLASIFICACION, COLOR_NULO
        ), 
    levels = c(
        "1000-MC bosques", "2000-Agricultura", "3000-Ganad y past",
        "4000-Zonas humed", "5000-Infraest", "6000-Otras tierras",
        "7000 No clasif", "Nulo"
    )
)

leaflet() |>
    addTiles(group = "OSM") |>
    addProviderTiles(
        provider = providers$Esri.WorldImagery, 
        group = "ESRI World Imagery"
    ) |>
    addProviderTiles(
        provider = providers$CartoDB.DarkMatter,
        group = "Dark Matter"
    ) |>        
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t1_uso),
        color = ~colores(geo_puntos$t1_uso),        
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Uso en t1: </strong>", geo_puntos$t1_uso),
      paste0("<strong>Uso en t2: </strong>", geo_puntos$t2_uso),
      sep = '<br/>'
    ),      
        group = "Uso en t1"
    ) |>
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t1_uso,
    title = "Uso en t1",
    group = "Uso en t1"    
  ) |>      
    addCircleMarkers(
        data = geo_puntos,
        radius = 2,
        fillColor = ~colores(geo_puntos$t2_uso),
        color = ~colores(geo_puntos$t2_uso),        
        # clusterOptions = markerClusterOptions(),
    popup = paste(
      paste0("<strong>Parcela: </strong>", geo_puntos$plot_id),
      paste0("<strong>Muestra: </strong>", geo_puntos$sample_id),
      paste0("<strong>Uso en t1: </strong>", geo_puntos$t1_uso),
      paste0("<strong>Uso en t2: </strong>", geo_puntos$t2_uso),
      sep = '<br/>'
    ),      
        group = "Uso en t2"
    ) |>    
  addLegend(
    position = "bottomleft",    
    pal = colores,
    values = geo_puntos$t2_uso,
    title = "Uso en t2",
    group = "Uso en t2"    
  ) |>  
    addLayersControl(
        baseGroups = c("OSM", "ESRI World Imagery", "Dark Matter"),
        overlayGroups = c("Uso en t1", "Uso en t2"),
    )

Distribución de puntos de muestreo en clases de uso

La distribución de los puntos de muestreo en clases de uso del suelo, en t1 y t2, se muestra en un gráfico de barras y en un gráfico de pastel, ambos elaborados con el paquete Plotly.

t1

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t1_uso), fill = fct_infreq(t1_uso))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_USO) +
    xlab("Uso") +
    ylab("Cantidad de puntos") +
    labs(fill = "Uso") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

# Total de puntos
conteo <- puntos |>
    count(t1_uso)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_USO[as.character(conteo$t1_uso)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t1_uso, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

t2

# Gráfico de barras ggplot2
grafico_barras <-
    puntos |>
    ggplot(aes(x = fct_infreq(t2_uso), fill = fct_infreq(t2_uso))) +
    geom_bar(
    aes(
      text = paste0(
        "Cantidad de puntos: ", after_stat(count)
      )     
    )
    ) +
    scale_fill_manual(values = COLORES_USO) +
    xlab("Uso") +
    ylab("Cantidad de puntos") +
    labs(fill = "Uso") +
    theme_clean() +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    legend.title = element_blank(),
    legend.box.background = element_blank(),
    legend.key = element_blank()
  )

# Gráfico de barras plotly
ggplotly(grafico_barras, tooltip = "text") |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)' # Establece el color del borde a transparente
    )
  ) 

# Total de puntos
conteo <- puntos |>
    count(t2_uso)

# Asegurarse de que los colores se asignen en el orden correcto
colores_para_plotly <- COLORES_USO[as.character(conteo$t2_uso)]

# Gráfico de pastel plotly
plot_ly(
    conteo, 
    labels = ~ t2_uso, 
    values = ~ n,
    marker = list(colors = colores_para_plotly)
) |>
add_pie()

Análisis con el paquete FIESTA

En esta sección se analizan los datos con el paquete FIESTA y su módulo de inventarios basados en fotografías (Photo-Based). Este módulo calcula estimaciones de población y errores de muestreo asociados. A diferencia de los estimadores tradicionales del libro verde de FIA (utilizados en otros módulos de FIESTA), que se construyeron basándose en el paradigma de muestreo finito utilizando parcelas de muestra con área distinta, los estimadores basados en fotos se construyeron en el contexto del paradigma de muestreo infinito, junto con el concepto de una región de soporte. FIESTA incluye estimadores no proporcionales para estimaciones de área y cobertura porcentual por dominio, y estimadores de media de cocientes para estimaciones de área y cobertura porcentual dentro del dominio, y soporta la post-estratificación para reducir la varianza.

Primero se analiza la variable de cobertura de la tierra y posteriormente la variable de uso del suelo. En ambos casos, se calcula el porcentaje de cada clase en cada parcela y el porcentaje total de cada clase en el área de estudio, tanto en t1 como en t2. Por último, se analiza el flujo de cambios de t1 a t2.

Cobertura de la tierra

Porcentaje de coberturas por parcela

Se utiliza la función datPBpnt2pct para transponer los datos de puntos a porcentajes en cada parcela de monitoreo.

t1

# Transposición de puntos a porcentajes de coberturas por parcela en t1
porcentaje_cobertura_parcela_t1 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t1_cobertura")

# Despliegue de la transposición de puntos a porcentajes de coberturas por parcela en t1
porcentaje_cobertura_parcela_t1 |>
  datatable(
    caption = "Porcentaje de coberturas por parcela en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

t2

# Transposición de puntos a porcentajes de coberturas por parcela en t2
porcentaje_cobertura_parcela_t2 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t2_cobertura")

# Despliegue de la transposición de puntos a porcentajes de coberturas por parcela en t2
porcentaje_cobertura_parcela_t2 |>
  datatable(
    caption = "Porcentaje de coberturas por parcela en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de coberturas

La función modPBpop() genera datos de poblaciones, los cuales se utilizan posteriormente para calcular el porcentaje total (en todas las parcelas) en t1 y t2.

# Generación de datos de poblaciones
PBpopdat <- modPBpop(pnt = puntos,
    pltassgnid = "plot_id",
    pntid = "sample_id")

t1

La función modPB() genera las estimaciones de los porcentajes y los errores de muestreo.

# Estimación de la distribución de coberturas en t1
LCt1 <- modPB(PBpopdat = PBpopdat, rowvar = "t1_cobertura")

results.LCt1 <- LCt1$est

# Despliegue de la distribución de coberturas en t1
results.LCt1 |>
  datatable(
    caption = "Distribución de coberturas en t1",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

t2

# Estimación de la distribución de coberturas en t2
LCt2 <- modPB(PBpopdat = PBpopdat, rowvar = "t2_cobertura")

results.LCt2 <- LCt2$est

# Despliegue de la distribución de coberturas en t2
results.LCt2 |>
  datatable(
    caption = "Distribución de coberturas en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Cambios porcentuales en coberturas de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id"
  )

# Cobertura en t1 vs cobertura en t2
coberT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_cobertura",
      colvar = "t2_cobertura"
    )

# Despliegue de datos en una tabla
coberT1vT2$est |>
  datatable(
    caption = "Cobertura en t1 vs cobertura en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Parte de la información de la tabla anterior se muestra en el siguiente gráfico de líneas.

LCt1.extract <- data.frame(LCt1$est$t1_cobertura, LCt1$est$Estimate)

# Eliminación de la fila "Total"
LCt1.extract <- 
    LCt1.extract |>
  filter(LCt1.est.t1_cobertura != "Total")

Year2006 <- as.data.frame(c(2006, 2006, 2006, 2006, 2006, 2006))
names(Year2006)[1] <- "Year"
LCt1.extract <- data.frame(LCt1.extract, Year2006)
names(LCt1.extract)[1] <- "Variable"
names(LCt1.extract)[2] <- "Value"


LCt2.extract <- data.frame(LCt2$est$t2_cobertura, LCt2$est$Estimate)

# Eliminación de la fila "Total"
LCt2.extract <- 
    LCt2.extract |>
  filter(LCt2.est.t2_cobertura != "Total")

Year2019 <- as.data.frame(c(2019, 2019, 2019, 2019, 2019, 2019))
names(Year2019)[1] <- "Year"
LCt2.extract <- data.frame(LCt2.extract, Year2019)
names(LCt2.extract)[1] <- "Variable"
names(LCt2.extract)[2] <- "Value"

LCt.1y2 <- rbind(LCt2.extract, LCt1.extract)

LCt.1y2$Value <- as.numeric(LCt.1y2$Value)

# Gráfico de líneas en ggplot2
grafico_lineas <-
    ggplot(data = LCt.1y2, aes(
        x = Year,
        y = Value,
        color = factor(Variable)
    )) +
    geom_line(size = 1) +
    scale_x_continuous(limits = c(2006, 2019), breaks = c(2006, 2019)) +
    scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
    xlab("Año de análisis") + 
    ylab ("Porcentaje de cobertura") +
  scale_colour_manual(values = 
    c("1000-Vegetacion"      = COLOR_VEGETACION,
      "2000-Sin vegetacion"  = COLOR_SIN_VEGETACION,
      "3000-Agua"            = COLOR_AGUA,
        "4000-Nubes y sombras" = COLOR_NUBES_Y_SOMBRAS,
        "5000-Sin informacion" = COLOR_SIN_INFORMACION,
        "Nulo"                 = COLOR_NULO
    )
  ) +
    labs(colour = "Cobertura de la tierra") +
    theme_clean() +
    theme(legend.position = "top")

ggplotly(grafico_lineas) |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)'
    )
  ) 

Cambios de áreas de coberturas de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Cobertura en t1 vs cobertura en t2
coberT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_cobertura",
      colvar = "t2_cobertura",
        tabtype = "AREA"
    )

# Despliegue de datos en una tabla
coberT1vT2$est |>
  datatable(
    caption = "Cobertura en t1 vs cobertura en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Las cantidades de la tabla anterior están expresadas en acres. La misma información se muestra en el siguiente gráfico de barras, pero en hectáreas (ha).

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Area de cobertura en t1
T1.COBERTURA.area <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_cobertura",
        tabtype = "AREA"
    )

# Area de cobertura en t2
T2.COBERTURA.area <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t2_cobertura",
        tabtype = "AREA"
    )

T1.COBERTURA.area$est <-
    T1.COBERTURA.area$est |>
    mutate(Estimate = as.numeric(Estimate))

T2.COBERTURA.area$est <-
    T2.COBERTURA.area$est |>
    mutate(Estimate = as.numeric(Estimate))

CambioNetoT1vT2 <- 
  data.frame(
    cobertura=T2.COBERTURA.area$est$t2_cobertura,
    estimacion_t1=T1.COBERTURA.area$est$Estimate,
    estimacion_t2=T2.COBERTURA.area$est$Estimate,
    cambio_neto.t1at2=T1.COBERTURA.area$est$Estimate-T2.COBERTURA.area$est$Estimate
  )

# Eliminación de la fila "Total"
CambioNetoT1vT2 <- 
    CambioNetoT1vT2 |>
  filter(cobertura != "Total")

# Conversión de acres a hectáreas
CambioNetoT1vT2 <- 
    CambioNetoT1vT2 |>
  mutate(cambio_neto.t1at2 = cambio_neto.t1at2 * 0.40468564)

# Gráfico de barras ggplot2
grafico_barras <-
  ggplot(
    data=CambioNetoT1vT2,
    aes(
      x=reorder(cobertura, cambio_neto.t1at2), 
      y=cambio_neto.t1at2,
      fill=cobertura
    )
  ) +
  geom_bar(stat="identity") +
  scale_y_continuous(limits = c(-5000, 7000), breaks= seq(-6000, 6000, 1000)) +
  coord_flip() +
  theme_igray() +
  scale_fill_manual(values = COLORES_COBERTURA) +
  theme(axis.text.x = element_text(angle = 40, hjust = 0.9))+
  labs(title = "", x = "Coberturas", y = "Área (ha)")
  
# Gráfico de barras plotly
ggplotly(grafico_barras) |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)'
    )
  )   

Transiciones de t1 a t2

El siguiente gráfico de Sankey muestra las cantidades de puntos que han cambiado en las clases entre t1 y t2.

puntos_sankey <-
    puntos |>
    select(sample_id, t1_cobertura, t2_cobertura)

# write_csv(puntos_sankey, "puntos_sankey.csv")
# puntos_sankey <- read_csv("puntos_sankey.csv")

# El prefijo es para diferenciar t1 y t2 en el gráfico de Sankey
puntos_sankey <- 
    puntos_sankey |>
  mutate(
    t1_cobertura = paste0("T1-", t1_cobertura),
    t2_cobertura = paste0("T2-", t2_cobertura)
  )

label <- c("T1-1000-Vegetacion", "T1-2000-Sin vegetacion", "T1-3000-Agua", 
                     "T1-4000 Nubes y sombras", "T1-5000-Sin informacion", "T1-Nulo",
                     "T2-1000-Vegetacion", "T2-2000-Sin vegetacion", "T2-3000-Agua", 
                     "T2-4000 Nubes y sombras", "T2-5000-Sin informacion", "T2-Nulo"
         )

color <- c(COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
           COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO,
           COLOR_VEGETACION, COLOR_SIN_VEGETACION, COLOR_AGUA, 
           COLOR_NUBES_Y_SOMBRAS, COLOR_SIN_INFORMACION, COLOR_NULO
         )

# Agregar columnas de códigos para usar en el gráfico de Sankey
puntos_sankey <- puntos_sankey |>
  mutate(
    t1_codigo = case_when(
      t1_cobertura == 'T1-1000-Vegetacion' ~ 0,
      t1_cobertura == 'T1-2000-Sin vegetacion' ~ 1,
      t1_cobertura == 'T1-3000-Agua' ~ 2,
      t1_cobertura == 'T1-4000 Nubes y sombras' ~ 3,
      t1_cobertura == 'T1-5000-Sin informacion' ~ 4,
      t1_cobertura == 'T1-Nulo' ~ 5
    ),
    t2_codigo = case_when(
      t2_cobertura == 'T2-1000-Vegetacion' ~ 6,
      t2_cobertura == 'T2-2000-Sin vegetacion' ~ 7,
      t2_cobertura == 'T2-3000-Agua' ~ 8,
      t2_cobertura == 'T2-4000 Nubes y sombras' ~ 9,
      t2_cobertura == 'T2-5000-Sin informacion' ~ 10,
      t2_cobertura == 'T2-Nulo' ~ 11
    )
  )

# Contar las combinaciones de t1_codigo y t2_codigo
conteo_combinaciones <- 
    puntos_sankey |>
  count(t1_codigo, t2_codigo)

# Crear las listas separadas
source <- conteo_combinaciones$t1_codigo
target <- conteo_combinaciones$t2_codigo
value <- conteo_combinaciones$n

# Gráfico de Sankey
grafico_sankey <- 
    plot_ly(
    type = "sankey",
    orientation = "h",
    node = list(
      label = label,
      color = color,
      pad = 15,
      thickness = 20,
      line = list(
        color = "black",
        width = 0.5
      ),
      font = list(
        size = 12,
        color = "#000000",
        weight = "bold"
      )
    ),
    link = list(
      source = source,
      target = target,
      value =  value
    )
  )

grafico_sankey <- 
    grafico_sankey |> 
    layout(
        title = list(
            text = "Transiciones entre t1 y t2"
        ),
    font = list(
      size = 12,
      color = "#000000",
      weight = "bold"
    )  
  )

grafico_sankey

Uso del suelo

Porcentaje de usos por parcela

t1

# Transposición de puntos a porcentajes de usos por parcela en t1
porcentaje_uso_parcela_t1 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t1_uso")

# Despliegue de la transposición de puntos a porcentajes de uso por parcela en t1
porcentaje_uso_parcela_t1 |>
  datatable(
    caption = "Porcentaje de usos por parcela en t1",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

t2

# Transposición de puntos a porcentajes de uso por parcela en t2
porcentaje_uso_parcela_t2 <- 
    datPBpnt2pct(puntos, uniqueid="plot_id", tvar="t2_uso")

# Despliegue de la transposición de puntos a porcentajes de usos por parcela en t2
porcentaje_uso_parcela_t2 |>
  datatable(
    caption = "Porcentaje de usos por parcela en t2",
    rownames = FALSE,
    options = list(
      pageLength = 5,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Porcentaje total de usos

# Generación de datos de poblaciones
PBpopdat <- modPBpop(pnt = puntos,
    pltassgnid = "plot_id",
    pntid = "sample_id")

t1

# Estimación de la distribución de usos en t1
LCt1 <- modPB(PBpopdat = PBpopdat, rowvar = "t1_uso")

results.LCt1 <- LCt1$est

# Despliegue de la distribución de usos en t1
results.LCt1 |>
  datatable(
    caption = "Distribución de usos en t1",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

t2

# Estimación de la distribución de usos en t2
LCt2 <- modPB(PBpopdat = PBpopdat, rowvar = "t2_uso")

results.LCt2 <- LCt2$est

# Despliegue de la distribución de usos en t2
results.LCt2 |>
  datatable(
    caption = "Distribución de usos en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Cambios porcentuales en usos de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id"
  )

# Uso en t1 vs uso en t2
usoT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_uso",
      colvar = "t2_uso"
    )

# Despliegue de datos en una tabla
usoT1vT2$est |>
  datatable(
    caption = "Uso en t1 vs uso en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Parte de la información de la tabla anterior se muestra en el siguiente gráfico de líneas.

LCt1.extract <- data.frame(LCt1$est$t1_uso, LCt1$est$Estimate)

# Eliminación de la fila "Total"
LCt1.extract <- 
    LCt1.extract |>
  filter(LCt1.est.t1_uso != "Total")

Year2006 <- as.data.frame(c(2006, 2006, 2006, 2006, 2006, 2006, 2006, 2006))
names(Year2006)[1] <- "Year"
LCt1.extract <- data.frame(LCt1.extract, Year2006)
names(LCt1.extract)[1] <- "Variable"
names(LCt1.extract)[2] <- "Value"


LCt2.extract <- data.frame(LCt2$est$t2_uso, LCt2$est$Estimate)

# Eliminación de la fila "Total"
LCt2.extract <- 
    LCt2.extract |>
  filter(LCt2.est.t2_uso != "Total")

Year2019 <- as.data.frame(c(2019, 2019, 2019, 2019, 2019, 2019, 2019, 2019))
names(Year2019)[1] <- "Year"
LCt2.extract <- data.frame(LCt2.extract, Year2019)
names(LCt2.extract)[1] <- "Variable"
names(LCt2.extract)[2] <- "Value"

LCt.1y2 <- rbind(LCt2.extract, LCt1.extract)

LCt.1y2$Value <- as.numeric(LCt.1y2$Value)

# Gráfico de líneas en ggplot2
grafico_lineas <-
    ggplot(data = LCt.1y2, aes(
        x = Year,
        y = Value,
        color = factor(Variable)
    )) +
    geom_line(size = 1) +
    scale_x_continuous(limits = c(2006, 2019), breaks = c(2006, 2019)) +
    scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10)) +
    xlab("Año de análisis") + 
    ylab ("Porcentaje de uso") +
  scale_colour_manual(values = 
    c("1000-MC bosques"    = COLOR_MC_BOSQUES,
      "2000-Agricultura"   = COLOR_AGRICULTURA,
      "3000-Ganad y past"  = COLOR_GANADERIA_Y_PASTOS,
        "4000-Zonas humed"   = COLOR_ZONAS_HUMEDALES,
        "5000-Infraest"      = COLOR_INFRAESTRUCTURA,
        "6000-Otras tierras" = COLOR_OTRAS_TIERRAS,
        "7000 No clasif"     = COLOR_NO_CLASIFICACION,
        "Nulo"               = COLOR_NULO
    )
  ) +
    labs(colour = "Uso del suelo") +
    theme_clean() +
    theme(legend.position = "top")

ggplotly(grafico_lineas) |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)'
    )
  ) 

Cambios de áreas de usos de t1 a t2

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Uso en t1 vs uso en t2
usoT1vT2 <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_uso",
      colvar = "t2_uso",
        tabtype = "AREA"
    )

# Despliegue de datos en una tabla
usoT1vT2$est |>
  datatable(
    caption = "Uso en t1 vs uso en t2",
    rownames = FALSE,
    options = list(
      pageLength = 10,
      language = list(url = '//cdn.datatables.net/plug-ins/1.10.11/i18n/Spanish.json')
    )
  )

Las cantidades de la tabla anterior están expresadas en acres. La misma información se muestra en el siguiente gráfico de barras, pero en hectáreas (ha).

PBpoparea <- 
    modPBpop(
      pnt = puntos,
      pltassgnid = "plot_id",
      pntid = "sample_id",
      unitarea = AREA_ESTUDIO
  )

# Area de uso en t1
T1.USO.area <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t1_uso",
        tabtype = "AREA"
    )

# Area de uso en t2
T2.USO.area <- 
    modPB(
        PBpopdat = PBpoparea,
      rowvar = "t2_uso",
        tabtype = "AREA"
    )

T1.USO.area$est <-
    T1.USO.area$est |>
    mutate(Estimate = as.numeric(Estimate))

T2.USO.area$est <-
    T2.USO.area$est |>
    mutate(Estimate = as.numeric(Estimate))

CambioNetoT1vT2 <- 
  data.frame(
    uso=T2.USO.area$est$t2_uso,
    estimacion_t1=T1.USO.area$est$Estimate,
    estimacion_t2=T2.USO.area$est$Estimate,
    cambio_neto.t1at2=T1.USO.area$est$Estimate-T2.USO.area$est$Estimate
  )

# Eliminación de la fila "Total"
CambioNetoT1vT2 <- 
    CambioNetoT1vT2 |>
  filter(uso != "Total")

# Conversión de acres a hectáreas
CambioNetoT1vT2 <- 
    CambioNetoT1vT2 |>
  mutate(cambio_neto.t1at2 = cambio_neto.t1at2 * 0.40468564)

# Gráfico de barras ggplot2
grafico_barras <-
  ggplot(
    data=CambioNetoT1vT2,
    aes(
      x=reorder(uso, cambio_neto.t1at2), 
      y=cambio_neto.t1at2,
      fill=uso
    )
  ) +
  geom_bar(stat="identity") +
  scale_y_continuous(limits = c(-5000, 7000), breaks= seq(-6000, 6000, 1000)) +
  coord_flip() +
  theme_igray() +
  scale_fill_manual(values = COLORES_USO) +
  theme(axis.text.x = element_text(angle = 40, hjust = 0.9))+
  labs(title = "", x = "Usos", y = "Área (ha)")
  
# Gráfico de barras plotly
ggplotly(grafico_barras) |> 
  config(locale = 'es') |>
  layout(
    showlegend = TRUE,
    legend = list(
      title = list(text = ''),
      bordercolor = 'rgba(0,0,0,0)'
    )
  )   

Transiciones de t1 a t2

puntos_sankey <-
    puntos |>
    select(sample_id, t1_uso, t2_uso)

# write_csv(puntos_sankey, "puntos_sankey.csv")
# puntos_sankey <- read_csv("puntos_sankey.csv")

# El prefijo es para diferenciar t1 y t2 en el gráfico de Sankey
puntos_sankey <- 
    puntos_sankey |>
  mutate(
    t1_uso = paste0("T1-", t1_uso),
    t2_uso = paste0("T2-", t2_uso)
  )

label <- c("T1-1000-MC bosques", "T1-2000-Agricultura", "T1-3000-Ganad y past", 
                     "T1-4000-Zonas humed", "T1-5000-Infraest", "T1-6000-Otras tierras",
                     "T1-7000 No clasif", "T1-Nulo",
                     "T2-1000-MC bosques", "T2-2000-Agricultura", "T2-3000-Ganad y past", 
                     "T2-4000-Zonas humed", "T2-5000-Infraest", "T2-6000-Otras tierras",
                     "T2-7000 No clasif", "T2-Nulo"
                    )

color <- c(COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
           COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS,
           COLOR_NO_CLASIFICACION, COLOR_NULO,
           COLOR_MC_BOSQUES, COLOR_AGRICULTURA, COLOR_GANADERIA_Y_PASTOS, 
           COLOR_ZONAS_HUMEDALES, COLOR_INFRAESTRUCTURA, COLOR_OTRAS_TIERRAS,
           COLOR_NO_CLASIFICACION, COLOR_NULO
                    )

# Agregar columnas de códigos para usar en el gráfico de Sankey
puntos_sankey <- puntos_sankey |>
  mutate(
    t1_codigo = case_when(
      t1_uso == 'T1-1000-MC bosques' ~ 0,
      t1_uso == 'T1-2000-Agricultura' ~ 1,
      t1_uso == 'T1-3000-Ganad y past' ~ 2,
      t1_uso == 'T1-4000-Zonas humed' ~ 3,
      t1_uso == 'T1-5000-Infraest' ~ 4,
      t1_uso == 'T1-6000-Otras tierras' ~ 5,
      t1_uso == 'T1-7000 No clasif' ~ 6,
      t1_uso == 'T1-Nulo' ~ 7
    ),
    t2_codigo = case_when(
      t2_uso == 'T2-1000-MC bosques' ~ 8,
      t2_uso == 'T2-2000-Agricultura' ~ 9,
      t2_uso == 'T2-3000-Ganad y past' ~ 10,
      t2_uso == 'T2-4000-Zonas humed' ~ 11,
      t2_uso == 'T2-5000-Infraest' ~ 12,
      t2_uso == 'T2-6000-Otras tierras' ~ 13,
      t2_uso == 'T2-7000 No clasif' ~ 14,
      t2_uso == 'T2-Nulo' ~ 15
    )
  )

# Contar las combinaciones de t1_codigo y t2_codigo
conteo_combinaciones <- 
    puntos_sankey |>
  count(t1_codigo, t2_codigo)

# Crear las listas separadas
source <- conteo_combinaciones$t1_codigo
target <- conteo_combinaciones$t2_codigo
value <- conteo_combinaciones$n

# Gráfico de Sankey
grafico_sankey <- 
    plot_ly(
    type = "sankey",
    orientation = "h",
    node = list(
      label = label,
      color = color,
      pad = 15,
      thickness = 20,
      line = list(
        color = "black",
        width = 0.5
      ),
      font = list(
        size = 12,
        color = "#000000",
        weight = "bold"
      )
    ),
    link = list(
      source = source,
      target = target,
      value =  value
    )
  )

grafico_sankey <- 
    grafico_sankey |> 
    layout(
    font = list(
      size = 12,
      color = "#000000",
      weight = "bold"
    )
  )

grafico_sankey

Recursos de referencia y consulta

• Forest Inventory Estimation and Analysis (FIESTA) (sitio oficial de FIESTA)
• FIESTA - Photo-Based Module (módulo de inventarios basados en fotografías de FIESTA)