INTRODUCCIÓN

La expansión de la frontera agrícola es uno de los mayores impulsores de degradación, fragmentación y deforestación de los bosques tropicales (^{Foley et al., 2011}; ^{Armenteras et al., 2017}; ^{Ceddia, 2019}). Para desarrollar planes de manejo del suelo y políticas que aborden estos procesos es necesario monitorear sus dinámicas espaciotemporales (^{Chowdhury, 2006}; ^{Pettorelli et al., 2014}; ^{Seymour y Harris, 2019}). La disponibilidad actual de instrumentos globales que producen datos espaciales, tales como MODIS Land Cover (^{Friedl et al., 2002}), GlobeLand30 (^{Jun et al., 2014}), o Global Forest Change (^{Hansen et al., 2013}), ha hecho que el análisis espacial y temporal de las dinámicas y procesos del suelo resulte ampliamente accesible. Particularmente, la base de datos de Global Forest Change (GFC) ha jugado un papel clave en la evaluación de tendencias de cobertura forestal, ya que puede ser usado para estimar la deforestación a escala de bioma (^{Hansen et al., 2013}), nacional (^{Sannier et al., 2016}; ^{Cunningham et al., 2019}; ^{Galiatsatos et al., 2020}) y subnacional (^{McRoberts et al., 2016}; ^{Linke et al., 2017}; ^{Milodowski et al., 2017}).

La base de datos de GFC contiene varias capas con una resolución espacial de 30 metros: la cobertura del dosel para el año 2000, un agregado de recuperación forestal para el periodo 2000-2012, y la pérdida anual de cobertura forestal. Estas características, junto con su alcance global y la actualización anual, han proporcionado a investigadores y tomadores de decisiones una herramienta confiable para monitorear las tendencias de cobertura forestal, especialmente en sitios donde las bases de datos regionales no existen o son insuficientes (^{Sannier et al., 2016}). Sin embargo, debido al algoritmo usado para producir la base de datos, en algunos casos hay diferencias locales en la estimación de áreas de pérdida o recuperación forestal, las cuales pueden llevar a una sub o sobrestimación del cambio forestal, así como a imprecisiones temporales en el año asignado a la pérdida forestal (^{Hansen et al., 2013}; ^{Galiatsatos et al., 2020}). Esta sobre/subestimación ha sido atribuida a características orográficas y climáticas del paisaje (^{Cunningham et al., 2019}); la presencia de cultivos arbóreos, y los clareos a pequeña escala, usualmente asociados a la agricultura minifundista (^{Tropek et al., 2014}; ^{Turubanova et al., 2018}; ^{Cunningham et al., 2019}). Las imprecisiones temporales han sido atribuidas principalmente a la presencia continua de nubes en los trópicos y a perturbaciones de pequeña escala (^{Milodowski et al., 2017}; ^{Turubanova et al., 2018}). Además, la definición de bosque y pérdida forestal usada en la base de datos de GFC no considera procesos de degradación forestal, los cuales frecuentemente superan a la deforestación en paisajes tropicales (^{Ellis et al., 2020}; ^{Matricardi et al., 2020}; ^{Qin et al., 2021}). Para reducir estas imprecisiones y evitar hacer inferencias incorrectas del área de cambio forestal a partir de estimaciones basadas en mapas, las bases de datos deben de ser validadas para cada caso particular utilizando estimaciones por muestreo (^{Hansen et al., 2014}; ^{Olofsson et al., 2014}; ^{Tropek et al., 2014}; ^{Sandker et al., 2021}).

Recientemente, ha habido un aumento en los estudios que usan la base de datos de GFC para evaluar la pérdida de cobertura forestal en México. Este ha sido usado para calcular la pérdida de hábitat y su efecto en la distribución de mamíferos (^{Klass et al., 2020}); para evaluar la deforestación asociada con la expansión agrícola y la extracción de madera (^{Blackman y Villalobos, 2021}; ^{Cho et al., 2021}); para evaluar la influencia de la tenencia de la tierra y el manejo del suelo sobre la extensión de cobertura forestal (^{Ellis et al., 2017}; ^{Lawrence et al., 2019}; ^{Haines, 2021}); y para evaluar la eficacia de distintas estrategias (por ejemplo, pago por servicios ambientales y áreas protegidas) en la conservación de la cobertura forestal (^{Sims y Alix-Garcia, 2017}; ^{Ramirez-Reyes et al., 2018}; ^{Izquierdo-Tort et al., 2019}; ^{Ellis, Sierra-Huelsz, et al., 2020}; ^{Powlen et al., 2021}; ^{González-Fernández et al., 2022}). Asimismo, la base de datos de GFC ha sido usada para evaluar tendencias e impulsores de deforestación a escalas nacionales (^{Portillo-Quintero y Smith, 2018}; ^{Figueroa et al., 2021}) y subnacionales (^{Krylov et al., 2018}). Sin embargo, de los estudios citados, solo tres utilizaron estimaciones por muestreo para calcular el área de pérdida forestal (^{Krylov et al., 2018}; ^{Ramirez-Reyes et al., 2018}; ^{González-Fernández et al., 2022}), mientras que en una hicieron una validación visual usando SIG (^{Portillo-Quintero y Smith, 2018}) y en otra llevaron a cabo una validación en campo (^{Izquierdo-Tort et al., 2019}). A pesar de que las estimaciones reportadas en la mayoría de estos estudios son confiables, la incertidumbre asociada con el uso de estimaciones basadas en mapas sin la validación a escala local podría dar resultados erróneos sobre la magnitud y tendencia del cambio forestal (^{Galiatsatos et al., 2020}).

En este estudio resaltamos la importancia de llevar a cabo estimaciones por muestreo a partir de la base de datos de GFC para reducir la incidencia de sub/sobreestimaciones del área de pérdida forestal y las imprecisiones temporales de cambio forestal. Realizamos la evaluación en cuatro regiones cafetaleras del este de México (Figura 1) donde el café es cultivado por minifundisttas bajo el dosel de agrobosques (~1 ha) en paisajes montañosos complejos (^{López-Arcadia, 2022}). Dadas las características orográficas y climáticas de la región, así como las prácticas de manejo de los cafeticultores, es muy probable que, en este paisaje, la precisión de la base de datos de GFC se vea comprometida (^{Milodowski et al., 2017}; ^{Turubanova et al., 2018}; ^{Cunningham et al., 2019}). Además, en 2013, la región fue afectada por la epidemia de la roya, lo que desencadenó un proceso de degradación y pérdida forestal ocasionada por la reducción de árboles de sombra en las parcelas de café (^{López-Arcadia, 2022}). Ya que esta epidemia coincide con el cambio en satélites usados para actualizar la base de datos de GFC (de Landsat 7 a Landsat 8), la evaluación de los datos puede también ayudar a separar los efectos de la epidemia de los ocasionados por la mayor sensibilidad del nuevo satélite.

Fuente: elaboración propia usando el censo cafetalero de SAGARPA para el año 2000.

Figura 1 Área de estudio. Las líneas continuas de color representan los límites administrativos de las regiones cafetaleras estudiadas, localizadas en el estado de Veracruz, México. El área sombreada es el área cafetalera de cada región donde la evaluación fue realizada. 

METODOLOGÍA

Siguiendo las recomendaciones de buenas prácticas para estimar áreas usando mapas (^{Olofsson et al., 2014}; ^{Tyukavina et al., 2015}), llevamos a cabo un muestreo probabilístico estratificado utilizando la versión 1.8 de la base de datos de GFC (Figura 2), la cual integra los datos originales con datos actualizados que abarcan desde el año 2000 hasta el 2020 (^{Hansen et al., 2013}). Para establecer el área inicial de cobertura forestal, consideramos como bosque aquellos pixeles que en el año 2000 tuvieron una cobertura del dosel mayor a 50% (^{Sannier et al., 2016}). Posteriormente, usando la capa lossyear, agrupamos los pixeles en tres categorías de acuerdo al cambio reportado en la cobertura forestal de 2001 a 2020: 1) No pérdida - pixeles donde no se detectó pérdida de cobertura forestal; 2) Posible pérdida - pixeles que no mostraron pérdida, pero son adyacentes a un pixel donde se detectó pérdida, y 3) Pérdida - pixeles donde se detectó pérdida de cobertura forestal. Posteriormente, generamos puntos de muestreo aleatorios para cada categoría usando la función ‘sampleRandom’ del paquete raster en R (^{Hijmans y van Etten, 2012}). Ya que nuestro mayor interés fue identificar áreas donde hubiera pérdida de bosque, para cada una de las cuatro regiones generamos 200 puntos para Pérdida, 100 puntos para Posible pérdida, y 100 puntos para No pérdida. Para verificar la pérdida de cobertura forestal, generamos mosaicos anuales (2001-2010) de los valores mínimos y máximos del índice de diferencia normalizada (NDVI) y los valores mínimos y máximos del índice de vegetación mejorada (EVI) para cada región utilizando imágenes Landsat 7 para el periodo 2001-2012 y Landsat 8 para el periodo 2013-2020. Todas las imágenes fueron generadas en la plataforma Google Earth Engine (^{Gorelick et al., 2017}). A partir de las imágenes resultantes, extrajimos el valor de cada punto de muestreo y elaboramos series de tiempo para ayudar en la verificación de la pérdida de cobertura forestal. Siguiendo a ^{Turubanova et al. (2018}), hicimos una verificación visual con imágenes de alta resolución de Google Earth para identificar el periodo en el que ocurrió la pérdida (Figura 3A). Después, usando las series de tiempo generadas para EVI y NDVI (Figura 3B), verificamos el año en el que ambos índices disminuyeron dentro del periodo identificado previamente y lo designamos como el año de pérdida observada. Consideramos la pérdida de cobertura forestal como la reducción de la cobertura del dosel de los agroecosistemas de café (por ejemplo, degradación), así como la remoción total de la cobertura forestal (i.e. deforestación), la cual fue usualmente remplazada por cultivos comerciales y pastizales. La verificación visual fue llevada a cabo por el primer autor dos veces para asegurar la precisión de la interpretación.

Figura 2 Metodología seguida para llevar a cabo las estimaciones de pérdida de cobertura forestal. En rojo se señalan las imágenes utilizadas y en azul las estimaciones calculadas. 

Figura 3 Ejemplo de imágenes usadas para verificar la pérdida de cobertura forestal. A. Establecimos visualmente el periodo en el cuál la pérdida de cobertura ocurrió (junio de 2014-julio de 2015 en este ejemplo). B. Usamos la serie de tiempo para verificar el año en el cual ambos índices disminuyeron dentro del rango de tiempo establecido en A (2015 en este ejemplo). 

Para calcular el área inicial de cobertura forestal para el año 2000 y el área de pérdida anual para el periodo 2001-2020 en cada región evaluada usamos la siguiente ecuación (^{Stehman, 2013}; ^{Turubanova et al., 2018}).

El error estándar fue estimado usando la siguiente ecuación:

Para evaluar si las estimaciones por muestreo de pérdida de cobertura forestal difirieron de las estimaciones basadas en el mapa, calculamos y comparamos el área de pérdida de cobertura forestal para el periodo 2001-2020 usando: 1) una estimación basada en el mapa a partir de los datos de GFC; 2) una estimación por muestreo utilizando solo las categorías de Pérdida y Posible pérdida; y 3) una estimación por muestreo considerando todas las categorías. Realizamos dos estimaciones por muestreo distintas para probar si la inclusión de la pérdida observada en las categorías de Pérdida posible y No pérdida tuvo un efecto en las estimaciones. Para comparar entre estimaciones, calculamos el cambio relativo del área estimada por los procedimientos por muestreo en relación con las estimaciones basadas en mapas. Para comparar los patrones temporales de pérdida, generamos series de tiempo de la pérdida anual de cobertura forestal para cada estimación.

Para las estimaciones basadas en muestreo, calculamos una media móvil de tres años para suavizar las fluctuaciones anuales de pérdida de cobertura forestal y ayudar a identificar las tendencias de cambio (^{Linke et al., 2017}). Para los años 2001 y 2020, la media móvil fue calculada usando solo el año posterior y el año anterior, respectivamente.

Dado que los datos usados fueron generados usando satélites con distintas sensibilidades para detectar cambios en la cobertura forestal (Landsat 7 de 2001 a 2012 y Landsat 8 de 2013 en adelante), verificamos: 1) si el año asignado a los pixeles marcados como Pérdida coincidieron con el año observado al realizar el muestreo; 2) si coincidieron dentro de ±1 año, o 3) si ocurrieron dos o más años antes o después de lo reportado. Con esta información evaluamos la precisión temporal de los datos de GFC para detectar la pérdida de cobertura forestal en el área de estudio.

RESULTADOS

El porcentaje de puntos correctamente asignado como Pérdida fue mayor del 90% para todas las regiones, con un valor mínimo de 92% (Zongolica) y un máximo de 98% (Huatusco). El porcentaje de puntos en la categoría de Posible pérdida (pixeles reportados como bosques adyacentes a pixeles de Pérdida) donde detectamos una pérdida de cobertura fue alto, entre 53% (Zongolica) y 65% (Córdoba). En cuanto a la categoría de No pérdida, el porcentaje de puntos donde observamos disminución fue considerable, con valores entre 23% y 26% (Huatusco y Coatepec, respectivamente; Cuadro 1).

El área de pérdida de cobertura forestal calculada a partir de la estimación basada en el mapa fue de 1433 ha (Coatepec), 699 ha (Córdoba), 1214 ha (Huatusco), y 388 ha (Zongolica). Al realizar las estimaciones por muestreo considerando solo la pérdida observada en las categorías de Pérdida y Pérdida posible, el área estimada se incrementó entre 348% y 396%. Al incluir la pérdida estimada en la categoría de No pérdida, la estimación de área de pérdida de cobertura forestal incrementó drásticamente, con valores que oscilaron entre 3057% y 7150% (Cuadro 2).

Al hacer la evaluación de las series de tiempo de pérdida de cobertura forestal, la estimación basada en el mapa mostró patrones similares para todas las regiones, con un periodo de pérdida relativamente estable de 2001 a 2012 y un incremento alto después de 2013 (Figura 4A). La estimación por muestreo usando las categorías de Pérdida y Pérdida posible (Figura 4B) muestra un patrón similar al de la estimación basada en el mapa, con un leve incremento principalmente entre 2001 y 2012. Al incluir la categoría de No pérdida, la estimación por muestreo mostró un incremento en la pérdida de cobertura forestal en el periodo de 2001-2012 (Figura 4C; Apéndice 1).

Figura 4 Estimaciones de área de pérdida de cobertura forestal (ha) para cada año para el periodo 2001-2020. A. Estimación basada en el mapa. B. Estimación por muestreo usando las categorías Pérdida y Posible Pérdida. C. Estimación por muestreo usando No pérdida, Pérdida y Posible Pérdida. La línea punteada vertical indica el año 2013, cuando hubo un cambio en el satélite usado para generar los datos de GFC. Para B y C, los puntos representan los valores estimados para cada año, la región sombreada el error estándar y la línea continua la media móvil calculada usando una ventana temporal de tres años. Nótese la diferencia en escalas para las tres estimaciones. 

El año de pérdida observada en el muestreo coincidió entre 49% y 64.5% con el año asignado en los datos del GFC (Cuadro 3). El porcentaje de coincidencias se incrementó entre 59.5% y 77.5% al considerar los años inmediatamente antes y después del año asignado en los datos originales. Hubo una cantidad considerable de puntos de muestreo (entre 15% y 28%) en los que la pérdida fue detectada dos o más años después de que ocurrió. Por otro lado, en una pequeña cantidad de puntos de muestreo, la pérdida fue detectada más de dos años antes de que tuviera lugar (de 1% a 5.5%).

DISCUSIÓN

La disponibilidad de datos espacialmente explícitos para el monitoreo de las dinámicas de cobertura del suelo ha permitido a los investigadores, profesionales y tomadores de decisiones a evaluar las tendencias locales, regionales y nacionales de cambio del suelo. Sin embargo, a causa de las metodologías usadas en la elaboración de estos datos (por ejemplo, la escala global oculta la heterogeneidad regional, las limitaciones de los datos de entrenamiento y el algoritmo utilizado), las estimaciones basadas en mapas pueden reportar tendencias de cambio imprecisas al ser usadas en paisajes tropicales complejos, donde el minifundismo es común. Al ser usadas en estos contextos, encontramos que las estimaciones basadas en mapas usando los datos de GFC pueden subestimar ampliamente el área de pérdida de cobertura forestal y los patrones temporales inferidos pueden estar sesgados.

La base de datos de GFC define la cobertura forestal como cualquier pixel con vegetación mayor a 5 m de altura, y la pérdida de bosque como la remoción de todos los árboles dentro de un pixel (^{Hansen et al., 2013}). Dada la importancia socioambiental de los agroecosistemas de café de sombra en nuestra área de estudio, consideramos como deforestación a la transición de agrobosque de café de sombra a otros cultivos (como en el caso de las plantaciones de plátano, cítricos, caña) y a la degradación forestal (por ejemplo, agroecosistemas de café degradados). Reconocemos que el incremento drástico en la pérdida de cobertura forestal que detectamos después de realizar las estimaciones por muestreo podría deberse a nuestra definición de deforestación. La definición, operacionalización y detección de la degradación forestal a través de técnicas de percepción remota representa un tema que ha sido pasado por alto en bases de datos globales, incluyendo el usado en nuestra evaluación (^{Hansen et al., 2014}; ^{Tropek et al., 2014}; ^{Ghazoul et al., 2015}; ^{Gao et al., 2020}). Sin embargo, dado el predominio de este proceso en los paisajes tropicales (^{Ellis et al., 2020}; ^{Matricardi et al., 2020}; ^{Qin et al., 2021}), su omisión podría llevar a conclusiones erróneas con respecto a las dinámicas de cambio en la cobertura forestal.

La escala subnacional de nuestra evaluación podría explicar las grandes diferencias entre las estimaciones basadas en el mapa y las estimaciones por muestreo. La base de datos fue elaborada con una perspectiva global, y hay evidencia de que se desempeña mejor a escalas supranacionales o nacionales (^{Hansen et al., 2013}; ^{Tyukavina et al., 2015}; ^{Turubanova et al., 2018}). Sin embargo, su implementación a escalas locales-regionales representa un reto importante (^{McRoberts et al., 2016}; ^{Milodowski et al., 2017}; ^{Galiatsatos et al., 2020}; ^{Kinnebrew et al., 2022}). Por consiguiente, el uso de estimaciones por muestreo puede ayudar a superar esta problemática, especialmente en áreas donde es escasa la disponibilidad de mapas locales de cobertura del suelo o están desactualizados (^{McRoberts et al., 2016}).

El porcentaje de falsos positivos (puntos de muestreo que indicaron deforestación en lugares donde no ocurrió) fue similar (2% a 8%) a lo reportado para Latinoamérica (alrededor del 4%; ^{Tyukavina et al., 2015}), Brasil (7.14%) e Indonesia (4.29%; ^{Turubanova et al., 2018}), pero menor a lo informado para el bioma tropical (13%; ^{Hansen et al., 2013}). Sin embargo, el número de puntos donde observamos pérdida en la categoría de Pérdida posible (pixeles adyacentes a aquellos reportados como Pérdida), fue mucho mayor que en otros estudios, con valores que oscilaron de 53% a 65%, en comparación con 25% a 39.9% en otros estudios (^{Turubanova et al., 2018}). Esto fue similar para el número de falsos negativos (puntos de muestreo que no reportaron pérdida forestal en lugares donde sí ocurrió) y, de hecho, fue mayor que lo observado en otros estudios: 23% a 26%, más alto que lo reportado en el bioma tropical (16.9%, de acuerdo con ^{Hansen et al., 2013}) y en Latinoamérica (17%; ^{Tyukavina et al., 2015}). Además, estos valores son drásticamente mayores que los reportados para Brasil (0.44%), la República Democrática del Congo (0.4%) e Indonesia (0.36%; ^{Turubanova et al., 2018}). La gran proporción de falsos negativos encontrados en este estudio es probablemente el resultado de las características del paisaje. La complejidad topográfica y la presencia continua de nubes en el área de estudio, características típicas de muchos paisajes montañosos en México y Centroamérica (^{Karmalkar et al., 2011}; ^{Cunningham et al., 2019}), podrían disminuir la capacidad de la base de datos de detectar la pérdida de cobertura forestal. Adicionalmente, el paisaje está compuesto de un mosaico de distintas coberturas creadas por pequeñas unidades de producción cuyos cambios producen perturbaciones muy pequeñas para ser detectadas con precisión (^{Milodowski et al., 2017}; ^{Turubanova et al., 2018}; ^{Cunningham et al., 2019}).

La congruencia de la escala temporal (años) de pérdida de cobertura forestal predicha y observada fue de 49%-64.5% para coincidencias exactas y aumentó a 59.5%-77.5% para coincidencias dentro ±1 año. Estos valores son notoriamente más bajos que lo reportado por ^{Hansen et al. (2013)}, quien encontró que el 75.2% de los puntos de muestreo fue asignado correctamente, y el 96.7% conferido dentro de ±1 de la pérdida observada. Asimismo, ^{Linke et al. (2017)} encontraron que en un bosque templado, 84.6% puntos fueron atribuidos correctamente para el año en que ocurrió la pérdida de cobertura, mientras que 92.1% fueron correctamente asignados dentro del rango de un año. Sin embargo, ^{Milodowski et al. (2017)} encontraron que los datos no lograron detectar perturbaciones pequeñas (<2 hectáreas), reportándolas como pérdida solo cuando las perturbaciones se habían expandido en años posteriores, sesgando por lo tanto la detección de pérdida de cobertura forestal a años posteriores. La mayoría de las imprecisiones temporales encontradas en nuestro estudio consistieron en la detección de pérdida más de dos años después de la pérdida observada. Por consiguiente, es posible que la presencia de pequeñas perturbaciones en la región, así como el uso de distintos satélites en los últimos años del periodo estudiado hayan inclinado los patrones temporales de pérdida forestal hacia años más recientes. Esto es visible en las series de tiempo, donde hay un cambio notable en el patrón de pérdida anual en el periodo 2001-2012, antes del cambio de satélites.

El muestreo probabilístico usado en este estudio permite reportar tendencias de pérdida de cobertura forestal más precisas; sin embargo, dicha metodología presenta algunos retos relacionados con el diseño muestral y la información espacial (^{Stehman y Foody, 2019}). Si bien es recomendable calcular el tamaño muestral por medio de fórmulas pertinentes, en la práctica, el tamaño muestral está influido por elecciones subjetivas relacionadas con los objetivos buscados y las capacidades técnicas y económicas del equipo de trabajo (^{Olofsson et al., 2014}; ^{Boschetti et al., 2016}; ^{Stehman y Foody, 2019}). Por otra parte, la asignación de los puntos de muestreo a las distintas categorías de suelo presentes en el mapa depende de un compromiso entre los distintos objetivos buscados (Stehman, 2012; ^{Wagner y Stehman, 2015}; ^{Stehman y Foody, 2019}). Así, el uso del muestreo probabilístico estratificado será más apropiado cuando el objetivo principal del trabajo sea el cálculo de estimaciones no sesgadas de cambio de uso del suelo, pero menos apropiado para calcular la precisión general del mapa de referencia (^{Stehman y Foody, 2019}). Por último, los análisis basados en muestreos probabilísticos implican la pérdida de la información espacial contenida en el mapa de referencia, lo que impide generar análisis espaciales a partir de las estimaciones calculadas (^{Turubanova et al., 2018}). Recomendamos tener en cuenta lo anterior al usar muestreos probabilísticos para el análisis de dinámicas de cobertura del suelo.