INTRODUCCIÓN
Una gran variedad de ecosistemas se encuentran amenazados por el cambio climático, la contaminación, la deforestación, la degradación de suelos, la escasez de agua, la generación de residuos, así como pérdida de biodiversidad relacionada con las actividades humanas (Beraldi-Campesi, 2014). El análisis de acervos de imágenes satelitales representa una alternativa económica, aunque indirecta, para monitorear a mediano y largo plazo el estado que guarda la diversidad de plantas y animales en diferentes ecosistemas terrestres. Por ejemplo, a nivel global se han documentado las tendencias de deforestación y degradación a partir de analizar índices -proxy del estado real- de vegetación, y mediciones de temperatura, precipitación y otras variables climáticas derivadas de imágenes satelitales (Yang et al., 2013; Vicente-Serrano et al., 2013; Liu et al., 2015).
Debido a su resolución temporal (16 días) y espacial (hasta 30 m), la familia de satélites Landsat ha brindado múltiples casos de éxito para analizar e interpretar fenómenos naturales o impactos ambientales a nivel regional durante más de cuatro décadas (Loveland y Dwyer, 2012). Por ejemplo, Duveiller et al. (2008), Broich et al. (2011) y Ernst et al. (2013) estudiaron zonas tropicales con el interés de monitorear -y evaluar cuando fue posible- deforestación, reforestación, degradación y regeneración de la vegetación. DeVries et al. (2015) desarrollaron un sistema de monitoreo de cambios en la vegetación en algunas selvas de Etiopía, mientras que Zhu et al. (2016) utilizaron imágenes Landsat 5, 7 y 8, para determinar tendencias de verdor en China.
Distintas áreas de México también han sido objeto de estudios basados en imágenes Landsat. En particular, se ha estudiado la dinámica espacio-temporal del NDVI y sus aplicaciones. Se conocen, por ejemplo, estudios de tendencias de NDVI Landsat 5 sobre manglares del noroeste de México (Flores-Cárdenas et al., 2017); análisis de cambio en series anuales de NDVI con Landsat 4, 5 y 7 sobre la Sierra Madre Occidental en Durango (Novo-Fernández et al., 2018); análisis de cambios en el uso de suelo y cobertura en el área metropolitana de Tepic-Xalisco establecidos a través de imágenes Landsat multi-anuales (Ávalos-Jiménez et al., 2018); y estudios de la variabilidad de NDVI Landsat-7 en algunas zonas de Baja California Sur (Maldonado-Enríquez et al., 2020).
En años recientes, la investigación científica sobre regiones selváticas y boscosas tropicales ha experimentado un creciente interés. Se ha cuantificado, por ejemplo, el elevado peligro de incendio al que están sometidos estos ecosistemas, derivado por una parte de huracanes que aportan carga de combustibles (pedazos de madera, hojarasca, etc.) y por otra parte los factores antropogénicos (cambio de uso de suelo, mal uso del fuego) que resulta en deforestación y fragmentación (Rodríguez-Trejo et al., 2011). Aunque el fuego no es la única causa de la degradación del suelo y la vegetación, para países tropicales la degradación de sus selvas se traduce en importantes impactos ecológicos y económicos (Armenteras et al., 2020).
Es plausible que el primer paso para el establecimiento de sistemas robustos de monitoreo de la degradación y recuperación de la vegetación de las selvas y bosques tropicales sea la clasificación cuantitativa de tendencias de NDVI en estas zonas. En el contexto mexicano, el 48% de estas se encuentran en los estados del sureste; 32.11% en la península de Yucatán, de los cuales el 12.43% se ubican en el estado de Campeche (INEGI, 2018). Por su condición de zona tropical, la alta frecuencia de nubes en la península de Yucatán representa un reto mayor para los estudios a largo plazo basados en imágenes satelitales (Mitchard et al., 2012); esta característica pudiera explicar los escasos artículos científicos dedicados al monitoreo de los cambios paulatinos y abruptos de la vegetación de esta región.
OBJETIVO
Identificar las tendencias globales y locales en una serie de tiempo NDVI Landsat de 2014 a 2020 en una zona de interés en la península de Yucatán, México.
MATERIALES y MÉTODOS
Área de estudio y cubo de datos
El estudio se realizó en la región oriental de Escárcega, Campeche (Figura 1) en la península de Yucatán (90°3’47.01” y 90°25’37.88” O, 18°59’45.01” y 18°11’50.11” N). La región de interés colinda con la Reserva de la Biosfera Calakmul, un territorio con presión antropogénica por el cambio de uso de suelo, debido a la expansión de la frontera agropecuaria del sureste (Velasco y Torres, 2019).
Datos del inventario nacional de usos de suelo y tipos de vegetación del Instituto Nacional de Estadística y Geografía (INEGI, 2018; https://www.inegi.org.mx/app/biblioteca/ficha.html?upc=889463842781), han evidenciado que en el área de estudio hay una superficie importante de selva en estado conservado. Se utilizó la información de este inventario concerniente a la clasificación de tipo de suelo y vegetación como material auxiliar en esta investigación, aunque al analizar las tendencias se excluyeron del área de interés los polígonos correspondientes a asentamientos humanos y cuerpo de agua.
Los tipos de selva que se encuentran en la zona son selva baja espinosa subperennifolia y selva mediana subperennifolia. Ambos tipos de vegetación, en estado primario y secundario (arbustivo y arbóreo), tienen alturas del estrato arbóreo clasificado como bajo (4 a 14 m) y mediano (15 a 30 m). El componente arbóreo de esta vegetación pierde entre 25-50% de su follaje estacionalmente. Esta característica está determinada por la composición florística y estructura de cada subregión, variables meteorológicas o, incluso, como consecuencia de la actividad antropogénica (Van Schaik, Terborgh y Wright, 1999; Zamora-Crescencio et al., 2017). Esta vegetación se ubica a una altitud de 2 a 255 m. Según García (2001), el clima es considerado Aw1 (cálido subhúmedo), con un rango de temperatura anual de 24 °C a 26 °C. La precipitación anual es de 1000 a 1500 mm (Cuervo-Robayo et al., 2014).
La evidencia empírica sobre la proximidad de mediciones in-situ de la productividad primaria neta y el NDVI ha permitido que este índice, atribuido a Tucker (1979), haya sido empleado en una copiosa cantidad de estudios sobre monitoreo de la vegetación. Por completitud, se recuerda que
donde NIR y RED hacen referencia a los valores de reflectancia del infrarrojo cercano y reflectancia del rojo, respectivamente. Dependiendo del sensor, las bandas RED y NIR se ubican en posiciones distintas, por ejemplo, para las imágenes Landsat-7 ETM+, RED y NIR están en las bandas 3 y 4, respectivamente, mientras que para el sensor Landsat 8 OLI estas se encuentran en las bandas 4 y 5.
Este tipo de índices de vegetación espectral se consideran como una medida de abundancia de la clorofila y absorción de energía (Myneni et al., 1995). Por construcción, el NDVI toma valores en el intervalo [-1,1]: valores mayores a 0.4 son típicamente relacionados a una mayor actividad fotosintética mientras que valores menores o cercanos a cero suelen reportarse en áreas con poca o nula cobertura vegetal.
El análisis de la calidad de las imágenes, la correspondiente imputación estadística de datos faltantes, el análisis de tendencias y la generación de gráficos para resumir algunos resultados fueron efectuados con diversos paquetes del lenguaje de programación R (R Core Team, 2022).
Procesamiento de imágenes
Se utilizaron imágenes Landsat 8 OLI del path/row 20/47 durante el periodo 2014-2020. Estas imágenes fueron descargadas gratuitamente del repositorio GloVis (GloVis, 2021; https://glovis.usgs.gov/). La resolución espacial (30 m) y temporal (16 días) hacen ideal el uso de este acervo para monitoreo remoto de la vegetación.
A pesar de la conocida falla del Scan Line Corrector del sensor Landsat-7 ETM+, la cual representa la pérdida de cerca del 22% de información en cada imagen obtenida después de mayo de 2003 (Zhang, Li, y Travis, 2007), las bandas espectrales de este sensor tienen la mayor cercanía con aquellas de Landsat 8 OLI. Por ejemplo, las bandas azul, verde y roja de ETM+ coinciden en 98.83, 98.30 y 98.9% con las correspondientes bandas espectrales de Landsat 8 OLI. Además, el intervalo de longitud de onda que define a la banda NIR de 8 OLI está contenido en aquel del NIR de Landsat-7 ETM+ (Roy et al., 2016).
Debido a que, en 2014, 2016, 2018 y 2019 solo fue posible adquirir 22 imágenes Landsat 8 OLI, aunado a que algunas de las imágenes originales presentaban a lo más 15% de datos disponibles (debido a la alta frecuencia de nubes en la zona), se optó por reemplazar 11 imágenes Landsat 8 OLI por las imágenes Landsat-7 ETM+ con la fecha de adquisición más cercana. De este modo se aseguró la existencia de 23 imágenes por año -condición deseable en estudios temporales basados en imágenes Landsat- y se mejoró la calidad del acervo estudiado, véase Tabla 1. Este tipo de estrategia de consolidación del acervo a estudiar es un procedimiento estándar (Zhu et al., 2016; Novo-Fernández et al., 2018; Xu y Guo, 2014).
Fecha Landsat 8 OLI | % d.d. | Fecha Landsat-7 ETM+ | % d.d. | % d.d.d.g. |
2/22/2014 | 14.04 | 2/14/2014 | 83.42 | 99.73 |
4/11/2014 | 12.38 | 4/3/2014 | 65.21 | 65.21 |
5/29/2014 | 0.22 | 5/21/2014 | 51.54 | 75.26 |
11/1/2014 | 0 | 11/13/2014 | 78.34 | 79.46 |
11/21/2014 | 0 | 11/29/2014 | 81.7 | 82.71 |
12/7/2014 | 3.53 | 12/15/2014 | 82.38 | 82.38 |
6/1/2015 | 0.44 | 6/9/2015 | 43.5 | 84.58 |
6/30/2016 | 0 | 6/27/2016 | 60.64 | 90.73 |
6/22/2017 | 0.39 | 6/30/2017 | 58.7 | 96.69 |
8/9/2017 | 0 | 8/1/2017 | 62.43 | 80.54 |
10/12/2017 | 3.68 | 10/4/2017 | 56.25 | 91.67 |
12/15/2017 | 0 | 12/7/2017 | 46.81 | 92.51 |
1/30/2018 | 0 | 1/24/2018 | 38.03 | 83.26 |
12/21/2019 | 0 | 12/29/2019 | 62 | 91.68 |
12/7/2020 | 0 | 12/15/2020 | 38.66 | 94.56 |
Nota: las letras d.g. abrevia la expresión “después de aplicar gapfill”.
Las imágenes Landsat-7 ETM+ fueron generadas a partir del Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS), con un nivel de procesamiento L1T de reflectancia de la superficie (Vermote et al., 1997; Masek et al., 2006). La corrección radiométrica de las imágenes Landsat 8 OLI, así como el enmascaramiento de nubes y sombras de nubes de todo el acervo (161 imágenes) se efectuó con el paquete RSToolbox de Leutner, Horning, y Schwalb-Willmann (2019).
Inicialmente, las imágenes Landsat-7 ETM+ (con el bandeo conocido) y las Landsat 8 OLI (con datos faltantes) fueron rellenadas con el algoritmo gapfill, el cual es un método espacio-temporal que permite rellenar datos faltantes en imágenes satelitales siempre y cuando algunas consideraciones teóricas se satisfagan (Gerber et al., 2018). La verificación automatizada de estas condiciones junto con el correspondiente rellenado de datos faltantes, está implementada en la función Gapfill() del paquete gapfill (Gerber, 2021) la cual fue utilizada en este estudio.
En aquellas zonas del acervo con muy alta frecuencia de datos faltantes, las condiciones teóricas del algoritmo no pudieron satisfacerse y, por tanto, esas zonas no pudieron ser rellenadas con gapfill. Para obtener una base de datos sin valores faltantes, al igual que ha ocurrido con otros estudios que han usado gapfill (Weiss et al., 2014; Sun et al., 2017; Li et al., 2018; Zhu et al., 2018), fue necesario aplicar una estrategia complementaria. En este estudio, se empleó el modelo clásico de interpolación lineal implementado en la función na_interpolation() del paquete imputeTS de Moritz y Bartz-Beielstein (2017).
Clasificación de tendencias en series de tiempo de NDVI
Como primer paso se aplicó bfast01 de De Jong et al. (2013) a los píxeles del cubo de datos para determinar estadísticamente la existencia de un punto de cambio en la tendencia de la serie de tiempo -a nivel píxel- de NDVI. Posteriormente, se aplicó la función bfast01classify sobre el objeto devuelto por bfast01 para clasificar el píxel por tipo de tendencia. Cuando bfast01 determina que no existe un punto de cambio en la serie de tiempo, bfast01classify aún puede determinar si existe una tendencia global, lineal y monótona (con pendiente positiva o negativa). Tanto bfast01 como bfast01classify están implementados en el paquete de R bfast (Verbesselt et al., 2021). A continuación, explicamos los detalles técnicos alrededor de la estrategia de análisis.
bfast01: determinando estadísticamente un cambio abrupto en la tendencia
Denotemos por y 𝑡 el valor del NDVI al tiempo 𝑡 y suponemos que se satisface la siguiente representación:
donde Tt y S t denotan los componentes de tendencia y estacionalidad, respectivamente; εt denota un error de medición el cual se supone aleatorio, independiente, con media cero y varianza común σ2. La letra N denota el número de observaciones en la serie de tiempo.
Aunque la representación de arriba es idéntica a la empleada en bfast de Verbesselt et al. (2010), bfast01 considera una tendencia con representación
Es decir, bfast01 considera como hipótesis inicial la existencia de una tendencia lineal global para la serie de tiempo y como hipótesis secundaria la existencia de dos tendencias lineales locales en igual número de segmentos definidos por medio de un cambio abrupto (punto de cambio) en el valor nominal de la tendencia. En comparación, la hipótesis secundaria de bfast considera múltiples segmentos y tendencias lineales locales. Esta simplificación hace posible que la estimación de bfast01 sea muy sencilla ya que solo es necesario resolver un problema de optimización para ajustar el modelo; bfast requiere resolver múltiples problemas de optimización de forma iterativa, lo que deriva en cómputo más intenso. Al enfocarse en encontrar un punto de cambio, bfast01 determina aquel cambio abrupto con mayor influencia sobre la serie de tiempo, es decir bfast01 reporta el cambio abrupto con mayor magnitud (valor absoluto del salto de la serie de tiempo en el punto de cambio estimado).
El procedimiento estadístico usado en bfast01 para determinar un cambio abrupto en la tendencia, depende de un ancho de banda, denotado como h, un número en el intervalo 0,1 y que puede interpretarse como el porcentaje de datos empleados en distintas partes del algoritmo. Como argumentan Tecuapetla-Gómez, Villamil-Cortez y Cruz-López (2021), el ancho de banda h en esencia controla la distancia mínima necesaria entre dos puntos de cambio consecutivos para que éstos sean detectados. En el presente estudio N = 161, por lo que usar h = 0.15 impica que la distancia mínima necesaria entre cualquier punto del tiempo y un punto de cambio estadísticamente significativo es aproximadamente 24 observaciones NDVI (1 año y 16 días, de acuerdo a la resolución temporal de Landsat). Observe que esto también implica, con este valor de ancho de banda, que no es posible determinar puntos de cambio anteriores al 2015 ni posteriores a 2019. Similarmente, si h = 0.23, entonces son necesarias 37 observaciones (1 año, 6 meses y 16 días) para poder identificar un cambio abrupto.
bfast01classify: clasificando tendencias
Una vez aplicado bfast01 a un pixel, el siguiente paso consistió en clasificar su tendencia. Esta tarea fue realizada con la función bfast01classify la cual, a partir de la tendencia estimada, clasifica píxeles en ocho posibles clases:
Enverdecimiento: tendencia lineal global creciente (con pendiente positiva)
Pardeamiento: tendencia lineal global decreciente (con pendiente negativa)
Enverdecimiento sostenido: tendencia lineal local creciente antes y después del cambio abrupto
Pardeamiento sostenido: tendencia lineal local decreciente antes y después del cambio abrupto
Enverdecimiento demorado: tendencia lineal local estable (pendiente cercana a cero) o creciente antes de un cambio abrupto hacia abajo seguido de una tendencia lineal local creciente
Pardeamiento demorado: tendencia lineal local estable o decreciente antes de un cambio abrupto hacia arriba seguido de una tendencia lineal local decreciente
Enverdecimiento a pardeamiento: tendencia lineal local creciente con un cambio abrupto seguido de una tendencia lineal local decreciente
Pardeamiento a enverdecimiento: tendencia lineal local decreciente con un cambio abrupto seguido de una tendencia lineal local creciente
Más detalles sobre la definición de estas clases pueden encontrarse en De Jong et al. (2013) y la documentación de bfast01classify en Verbesselt et al. (2021).
RESULTADOS
Para ejemplificar las definiciones presentadas en la sección anterior, la Figura 2 muestra los tipos de tendencias más comunes encontradas en el área de estudio. Como se presentan en la Figura 3, las tendencias de tipo 3 (enverdecimiento sostenido) y 4 (pardeamiento sostenido) tienen una presencia marginal en la zona.
Un estudio (no incluido) de robustez de bfast01 para clasificar las tendencias en función de los valores de ancho de banda h = 0.15 y h = 0.23, arrojó resultados cuantitativamente similares, por lo que a continuación describiremos las tendencias estimadas por bfast01 con h = 0.15 (este valor hace posible la identificación de cambios abruptos en el periodo 2015-2019).
Resulta de interés categorizar los píxeles del área en tipos de tendencias y usos de suelo y vegetación. Las Figuras 4, 5 y 6 resumen la categorización que se obtiene al usar el método de interpolación lineal para rellenar datos faltantes y bfast01 con ancho de banda h = 0.15 para estimar tendencias.
La selva mediana subperennifolia cubre el porcentaje más grande del área de estudio (35.82%); de esta área, la porción más grande -25.55%- experimenta algún tipo de pardeamiento mientras que el restante 10.26% fue clasificado con algún tipo de enverdecimiento. La segunda clase más grande de vegetación es la selva baja espinosa subperennifolia -29.56%-, la cual presentó tres veces más zonas con pardeamiento (21.89%) que enverdecidas (7.67%). La vegetación secundaria arbórea de selva mediana subperennifolia cubre el 14.37% del territorio estudiado de los cuales 11.10% presentó pardeamiento y solo el 3.26% algún tipo de enverdecimiento. Similarmente, del 8.5% marcado como pastizal cultivado, 5.59% mostró pardeamiento en su NDVI contra el 2.95% de enverdecimiento. El patrón de observar más área clasificada con pardeamiento que con enverdecimiento se mantiene en el resto de los tipos de vegetación: en la vegetación secundaria arbustiva de selva mediana subperennifolia -5.36% del área-, 3.75% reporta pardeamiento contra 1.6% que presenta algún tipo de enverdecimiento; el resto de tipos de vegetación y uso de suelo cubre el 6.32% del área de estudio muestra 5.34% y 0.97% de pardeamiento y enverdecimiento, respectivamente.
Las Figuras 4, 5 y 6 revelan también que en gran parte del área de estudio encontramos tendencias globales de crecimiento y decrecimiento del NDVI. En efecto, el 47.87% del área ha sido clasificada con pardeamiento, es decir, tendencia de clase 2 (tendencia global lineal decreciente) mientras que el 15.39% ha sido clasificada con tendencia clase 1 (tendencia global lineal creciente), es decir, con enverdecimiento. La Figura 4 muestra que son prácticamente inexistentes las tendencias de las clases 3 (enverdecimiento sostenido) y 4 (pardeamiento sostenido).
Los píxeles de las clases pardeamiento demorado (20.66%), pardeamiento a enverdecimiento (6.04%), enverdecimiento demorado (5.26%) y enverdecimiento a pardeamiento (3.88%) presentan tendencias decrecientes y/o crecientes antes o después de un cambio abrupto. Esta dinámica de la tendencia puede interpretarse como signo de recuperación (pardeamiento a enverdecimiento y enverdecimiento demorado) o degradación (pardeamiento demorado y enverdecimiento a pardeamiento) del NDVI después de reportar un cambio abrupto en su valor. Con esto en mente, en los siguientes apartados se discute la distribución espacial y temporal de los tipos de tendencia mencionados en este párrafo respecto de los tipos de uso de suelo y vegetación en la zona.
Pardeamiento demorado
Estableciendo si ambos segmentos (anterior y posterior al punto de cambio) poseen tendencias significativas contribuye a entender si esta dinámica decreciente en los valores del NDVI es un fenómeno persistente en la zona. Para tal efecto, se analizaron los tipos de suelo y vegetación con pardeamiento demorado respecto a la significancia de las tendencias pre- y post-cambio abrupto en el NDVI
A partir de la Figura 7 se infiere que, en zonas con agricultura de temporal anual, bosque cultivado, sabana, tular, vegetación secundaria arbórea de selva baja y vegetación secundaria arbustiva de selva mediana subperennifolia existe poca o nula presencia de pardeamiento demorado. En particular, véase el panel B, los píxeles de estas categorías prácticamente no presentaron cambio abrupto en el NDVI de 2015 a 2019.
Se identificó que el 83% de píxeles clasificados con pardeamiento demorado muestra tendencias significativas y decrecientes antes y después de un cambio abrupto; la selva mediana subperennifolia (25.95%), la selva baja espinosa subperennifolia (24.41%) y la vegetación secundaria arbórea de selva mediana (13.99%) son los tipos de vegetación con mayor prevalencia de esta clase de tendencia (véase Figura 5). Estos resultados están en correspondencia con lo reportado por Velasco y Torres (2019).
Los años con mayor incidencia de cambio abrupto en NDVI fueron 2017 (19.15%) y 2018 (57%). Entre 2015 y 2019, los cambios abruptos se concentraron en la selva baja espinosa subperennifolia (30.92%) y la selva mediana subperennifolia (29.67%), Figura 7-B.
La distribución de los años con cambio abrupto en los diferentes tipos de vegetación pudiera estar relacionada con el crecimiento poblacional y la consecuente demanda de los habitantes por el cambio del uso de suelo forestal a campos agrícolas y ganaderos.
Pardeamiento a enverdecimiento
Este tipo de tendencia puede estar asociado a la siembra de cultivos o recrecimiento de plántulas después de haber realizado un aprovechamiento forestal. Además, es plausible que la actividad agropecuaria de la zona también contribuya a explicar este tipo de tendencia.
De los píxeles de esta clase, además de aquellos con ambos segmentos de tendencia significativa se consideraron de interés aquellos con únicamente tendencia significativa poscambio abrupto ya que este tipo de característica puede interpretarse como un signo de recuperación del NDVI.
Se encontró que cerca del 45% de los pixeles clasificados con pardeamiento a enverdecimiento reportaron tendencias significativas y crecientes poscambio abrupto. De estos, la mayoría se concentraron en la selva mediana subperennifolia (11.08%), vegetación secundaria arbórea de selva mediana (10.45%), pastizal (10.46%) y la selva baja espinosa subperennifolia (17.04%) (véase Figura 6).
Se observaron más pixeles con tendencia significativa antes del cambio abrupto que después de éste en áreas que experimentan una pérdida paulatina de vegetación. En particular, aunque marginalmente, esto es observado en el pastizal cultivado (6.81 pre- y 5.44% poscambio abrupto), la selva baja espinosa subperennifolia (5.83 pre y 4.82% poscambio abrupto) y la vegetación secundaria arbórea de selva mediana subperennifolia (6.68 pre y 5.83% poscambio abrupto).
De acuerdo con la Figura 8-B, esta clase presentó la mayor cantidad de píxeles con cambio abrupto en 2017 (34.74%), seguido de 2018 y 2016 con 22.27 y 20.66%, respectivamente. Del 2015 al 2019, las vegetaciones que reportaron la mayor cantidad de cambios abruptos fueron: selva mediana subperennifolia (25.28%), vegetación secundaria arbórea de selva mediana subperennifolia (21.32%), pastizal cultivado (21.03%) y selva baja espinosa subperennifolia (17.04%).
Enverdecimiento demorado
En analogía con la discusión presentada en la sección anterior, aquí nos enfocamos en aquellos píxeles con ambos segmentos de tendencia significativa junto con los píxeles que mostraron únicamente tendencia significativa poscambio abrupto.
De este análisis se desprende que la mayor parte de pixeles clasificados con enverdecimiento demorado (54%), mostraron tendencias significativas y crecientes post-cambio abrupto concentrándose la mayoría en la selva mediana subperennifolia (19.44%) y la selva baja espinosa subperennifolia (13.70%) (véase Figura 6).
Los años con mayor incidencia de cambio abrupto fueron 2017 y 2018 con el 43 y 20.75%, respectivamente. Estos resultados pueden estar asociados a factores como el periodo de sequía y el comienzo de las precipitaciones y su impacto en la actividad fotosintética de la vegetación.
Además de la selva mediana subperennifolia (32.7%) y la selva baja espinosa subperennifolia (28.70%), el pastizal (14.25) y la vegetación secundaria arbórea de selva mediana (13.06%) concentraron la mayor cantidad de píxeles con cambio abrupto de 2015 a 2019 (véase Figura 9-B).
Enverdecimiento a pardeamiento
En analogía con la discusión presentada en la sección anterior, aquí nos enfocamos en aquellos pixeles con ambos segmentos de tendencia significativa junto con los píxeles con únicamente tendencia significativa poscambio abrupto. En este caso, estas características pueden interpretarse como un signo de degradación del NDVI.
Más del 63% del área clasificada con enverdecimiento a pardeamiento reportó tendencias significativas y decrecientes poscambio abrupto siendo la selva baja espinosa subperennifolia (16.94%) y el tular (11.96%) las más afectadas (Figura 5). Se encontró que 2018 (31.25%) y 2019 (32%) fueron los años con más cambios abruptos. En el periodo 2015-2019, los tipos de vegetación con mayor presencia de cambio abrupto fueron la selva baja espinosa subperennifolia (27.7%), selva mediana subperennifolia (14.7%), tular (13.83%), vegetación secundaria arbórea de selva mediana subperennifolia (10.88%) y el pastizal (9.33) (Figura 10-B).
DISCUSIÓN
Los distintos trabajos enfocados en estimar tendencias de NDVI en zonas de México (Flores-Cárdenas et al., 2017; Novo-Fernández et al., 2018; Ávalos-Jiménez et al., 2018; Maldonado-Enríquez et al., 2020), se han centrado en describir tendencias lineales globales en zonas con baja densidad de nubes. A diferencia de estos artículos, el presente trabajo se enfocó en la estimación estadística de tendencias significativas, locales y globales sobre una región en la península de Yucatán.
Para rellenar los datos faltantes (inherente a la perpetuidad de nubes en la zona) de la serie de tiempo de NDVI, inicialmente se aplicó el método de predicción espacio-temporal gapfill (Gerber et al., 2018). Hasta donde sabemos, este artículo es el primero en emplear gapfill para rellenar series de tiempo de NDVI del territorio mexicano. La alta frecuencia de datos faltantes en ciertas zonas del acervo impidió que algunas consideraciones teóricas del algoritmo se cumplieran por lo que fue necesario completar el rellenado de las imágenes utilizando interpolación lineal (temporal).
Tecuapetla-Gómez, Villamil-Cortez y Cruz-López (2021) utilizaron el algoritmo bfast para determinar múltiples puntos de cambio en una serie de tiempo de NDVI Landsat-7 con el objetivo de ligar estos cambios con regímenes de severidad de quema. En el presente artículo se utilizó el algoritmo bfast01 para estimar el punto de cambio más grande en una serie de tiempo de NDVI (imágenes Landsat-7 ETM+ y Landsat 8 OLI) y clasificó posteriormente la tendencia encontrada.
La mayor parte del área estudiada presenta valores de NDVI con tendencias lineales globales. De estas tendencias, la mayor parte fue clasificada con pardeamiento (47.87%) y en menor proporción con enverdecimiento (15.39%). El pardeamiento en la vegetación de selva puede estar relacionado con actividades como la agricultura, ganadería, obras de infraestructura y la extracción selectiva de especies arbóreas. Estas actividades antropogénicas son causa principal de la transformación, degradación y destrucción de la vegetación natural (Zamora-Crescencio et al., 2017). Definir cuáles son las actividades que ocasionan cambios negativos en la cobertura vegetal requiere del establecimiento de redes de monitoreo como proponen Gaitán et al. (2021) que permitan determinar el origen de las clases de cambio abrupto de pardeamiento y la validación de las tendencias. Esto permitirá a futuro desarrollar estrategias para mantener la estructura y funcionalidad de la selva en el largo plazo.
Los resultados resaltan que el 35.84% del área estudiada presenta tendencias lineales locales: pardeamiento demorado (20.66%), pardeamiento a enverdecimiento (6.04%), enverdecimiento demorado (5.26%) y enverdecimiento a pardeamiento (3.88%).
En la vegetación subperennifolia es plausible que los cambios pardeamiento a enverdecimiento y enverdecimiento demorado sean indicadores de un comportamiento fenológico atípico como respuesta al cambio climático, un posible desfase de las precipitaciones y una sequía prolongada. Por otra parte, en la selva caducifolia, el pardeamiento demorado puede deberse a la caída anticipada de las hojas de algunas especies; se ha documentado esta caída mucho antes de producirse un estrés hídrico severo (Van Schaik, Terborgh y Wright, 1993).
En las últimas tres décadas ha habido esfuerzos plausibles por caracterizar la fenología de las selvas bajas (Reich, 1995; González-Esquinca et al., 2016), incluso la de algunas especies de la Península de Yucatán (Valdez-Hernández et al., 2010). Consideramos que el presente estudio provee una herramienta de análisis exploratorio para determinar regímenes de cambio severo en la fenología de la vegetación a nivel regional. Además de allanar el camino para estudios posteriores sobre la recuperación y degradación de la cobertura forestal, en el futuro, este tipo de estudios deberá ser complementado con trabajo de campo para delinear técnicas fiables para el cuidado y protección de la vegetación de la Península de Yucatán (De Jong et al., 2013; Gaitán et al., 2015).
Alentamos a la comunidad interesada en este tipo de análisis a visitar el repositorio https://github.com/inder-tg/cpsClass donde encontrará el código utilizado para elaborar este artículo.
RECOMENDACIONES
La alta presencia de nubes en las zonas tropicales no debe ser un factor limitante para desarrollar estrategias de monitoreo remoto de la vegetación a mediano y largo plazo. Consideramos que la problemática inherente a los datos faltantes en las imágenes de este tipo de ecosistemas puede ser apropiadamente atendida, en una primera instancia, con herramientas como gapfill. En una segunda etapa, parece esencial rellenar por completo el acervo a analizar ya que, de otro modo, aún con garantías estadísticas, métodos de clasificación de tendencias como bfast01 pueden producir resultados distintos.