Predicción de la pérdida boscosa en el Área de Conservación La Amistad-Pacífico (ACLA-P), Costa Rica
Proyecto para acceder al Título de Especialista en Explotación de Datos y Descubrimiento del Conocimiento -versión preliminar-.
Universidad de Buenos Aires
Facultad de Ciencias Exactas y Naturales
Universidad de Buenos Aires
Facultad de Ciencias Exactas y Naturales
c
Resumen
En este documento se estima la probabilidad de deforestación de los próximos cuatro años del Área de Conservación la Amistad-Pacífico (ACLAP), Costa Rica, a una resolución de 30 x 30 m con el fin de facilitar un enfoque proactivo en la gestión forestal. Para esto, se utilizó un algoritmo de aprendizaje automático -gradient boosting- que empleó como insumo principal la base de datos de Global Forest Change (GFC), la cual incluye información anual de la pérdida de cobertura arbórea del 2001 a 2019 a nivel mundial.
Los resultados de la estimación muestran un área bajo la curva ROC -receiver operating characteristic- de 0.724 en un contexto donde la deforestación es baja - 1.4 % cada 4 años- y la pérdida de cobertura boscosa tiene comportamientos erráticos y difíciles de anticipar.
Asimismo, se determina que las variables más importantes para la discriminación de una futura pérdida arbórea son: cobertura boscosa en el año 2000, elevación y tendencia en los últimos 4 años de la cobertura arbórea. Por último, el modelo muestra diferentes grados de desempeño asociados al piso de la zona de vida.
A partir de los resultados, se construyó un mapa que muestra las probabilidades de un área a ser deforestada en el 2023 y otro mapa que agrupa tales probabilidades en zonas de riesgo alto y bajo.
Palabras clave: Global Forest Change, Costa Rica, deforestacion, inteligencia artificial, gestión forestal, Área de Conservación la Amistad-Pacífico, LightGBM, gradient boosting, aprendizaje automático.
1 Introducción
Los bosques albergan la gran mayoría de la biodiversidad mundial. En ellos habitan el 80% de las especies de anfibios, 75% de las especies de pájaros y el 68% de los mamíferos (FAO and UNEP 2020). A su vez, millones de personas en el mundo dependen de los bosques para su seguridad alimentaria y sus medios de subsistencia. Sumado a esto, los bosques juegan un papel preponderante en la lucha contra el cambio climático pues absorben y almacenan cerca del 30% de las emisiones de carbono de los combustibles fósiles y de la industria en su biomasa, suelos y madera (Brack 2019). En cuanto a la adaptación al cambio climático, estabilizan y enfrían los climas locales y brindan agua a las comunidades. Por ende, detener la degradación boscosa, la deforestación, y promover su restauración es vital para proteger la biodiversidad, estabilizar el clima, reducir la carga de gases de efecto invernadero, apoyar medios de vida y contribuir con el desarrollo sostenible.
La mayoría de las investigaciones sobre deforestación se concentran principalmente en determinar sus causas e impactos. También existe un desarrollo importante en procesos de teledetección para la identificación precisa y oportuna de sitios deforestados. Sin embargo, existen pocos estudios que pronostiquen la deforestación. La presente investigación aporta al último punto al implementar un algoritmo de inteligencia artificial - Light Gradient Boosting Machine - que predice la probabilidad de que un área sea deforestada en los próximos 4 años, en el ACLAP ubicada en el Pacífico de Sur de Costa Rica.
Para la estimación del modelo se utiliza la información GFC que cuenta, entre otras cosas, con información a nivel mundial de la cobertura y pérdida arbórea desde el año 2000 hasta el año 2019. Además, se utilizan fuentes secundarias de información como la ubicación de Áreas Silvestres Protegidas, características topográficas, cercanía a carreteras, precipitación, precio de la tierra, índice de desarrollo social, entre otras.
Esta predicción permitirá focalizar esfuerzos en las zonas de deforestación esperada, optimizando varias aplicaciones como: delimitación y resguardo de zonas de protección ambiental, concientización en la población, pagos por servicios ambientales, entre otros. También aportará información para la estimación de impactos ambientales futuros y la generación de políticas públicas con un enfoque proactivo.
El documento está organizado de la siguiente manera: en la segunda sección se describe la metodología de la investigación. En la tercera sección se presentan análisis exploratorios de los datos. En la cuarta se muestran los resultados de las estimaciones del modelo y seguidamente se presentan las conclusiones del documento.
1.1 Objetivo
Predecir la probabilidad de una zona de 30 x 30 m de perder su cobertura arbórea en los próximos 4 años - hasta 2023- en el Área de Conservación La Amistad Pacífico con base en la información de pérdida boscosa anual del GFC e información complementaria de dominio público.
1.2 Área de estudio
Este trabajo se desarrollaroá en el Área de Conservación Corredor La Amistad Pacífico en Costa Rica. Costa Rica es un país centroamericano de 51 100 km \(^2\) ubicado entre Nicaragua y Panamá. A pesar de su corta extensión, alberga alrededor del 5% de la diversidad conocida en el mundo (Ministerio de Ambiente 2015) .
El Área de Conservación Corredor La Amistad Pacífico, con una extensión de 555 994 ha, se encuentra en la zona sur del país y alberga la mayor parte de la Cordillera de Talamanca en su vertiente del Pacífico. Dentro de ella existen zonas de baja altitud cercanas al nivel del mar y otras que ascienden hasta el punto más alto del país - 3821 m.s.n.m-.
Figura 1.1: Área de estudio.
Su superficie está formada por gran diversidad de ecosistemas que incluyen las zonas de vida: bosque muy húmedo premontano y tropical, páramo pluvial subalpino, bosque pluvial montano bajo, montano y premontano (Bolaños, Watson, and Tosi 2005), las precipitaciones medias van desde los 1500 mm a los 6000 mm anuales (Atlas Climatológico de Costa Rica 2011).
De acuerdo con la información del GFC, el ACLAP poseía un 81 % -450 426 ha- de su territorio con una cobertura arbórea superior al 30% en el año 2000. De 2001 a 2019, ACLAP ha perdido 30 133 ha de su cobertura arbórea. Esto equivale a un 6.7 % de la cobertura arbórea encontrada en el año 2000.
Figura 1.2: Pérdida de cobertura boscosa en ACLAP del 2000 al 2019.
Tal y como se observa en la figura 1.2, la pérdida de cobertura boscosa es variable en el tiempo. Si bien la recta de mejor ajuste tiene una pendiente negativa su valor no es estadísticamente significativo - p valor = 0.24-.
A continuación, se muestra un mapa del ACLAP en donde se puede observar la cobertura boscosa en el año 2019 según GFC. En amarillo se muestran las zonas con una cobertura no boscosa, en celeste las zonas con cobertura boscosa y en rojo aquellas zonas que fueron bosque en el año 2000 pero al año 2019 no lo son.
Figura 1.3: Estado de la cobertura boscosa en el 2019
Es importante destacar que la pérdida de cobertura boscosa de los últimos años corresponde principalmente a degradación y deforestación a pequeña escala. Deforestación sistemática por parte de la industria maderera o agrícola no se encuentra tan presente como en otras latitudes. Por estas razones, el problema de predicción de la pérdida boscosa se torna especialmente difícil.
2 Metodología
En este apartado se describirá el procedimiento realizado para entrenar y aplicar un algoritmo de aprendizaje automático supervisado que predice la deforestación en ACLAP mediante la identificación de relaciones multivariadas complejas.
Este proceso inició con la obtención y análisis de fuentes de información disponible y la creación de variables por parte del autor mediante un procedimiento conocido como “feature engineering”. Seguidamente, por tratarse de un problema de índole espacial y temporal, se dividió el set de datos en cuatro periodos de análisis para entrenar, validar y probar el modelo. Con estas particiones se especificó la estructura del algoritmo y el procedimiento utilizado para elegir el modelo que maximizó la métrica de interés en los conjuntos de validación.
El modelo seleccionado fue puesto a prueba con datos no utilizados en el proceso de entrenamiento y validación para estimar los resultados que podría obtener en el año 2023.
Finalmente, para favorecer su aplicación práctica, se generó una representación de los resultados en formato ráster en donde se muestran las probabilidades de las zonas a ser deforestadas y otro mapa que muestra el riesgo a deforestación en los próximos cuatro años.
Todos estos procedimientos fueron realizados en el lenguaje y entorno de programción R (R Core Team 2020), haciendo uso principalmente de los paquetes: lightgbm (Ke, Soukhavong, and Lamb 2020), tidymodels (Kuhn and Wickham 2020), tidyverse (Wickham et al. 2019), sf (Pebesma 2018) y raster (Hijmans 2020).
2.1 Fuentes de información
En los siguientes subapartados se describirán las fuentes de información que son geográficamente explícitas y por ende, representan polígonos, líneas o puntos en un espacio cartográfico.
2.1.1 Hansen Dataset
Corresponde a la base de datos más importante del proyecto puesto que incluye la variable a predecir -pérdida o no pérdida de bosque-. Esta información publicada por primera vez en el 2013 (Hansen et al. 2013) tiene un formato ráster, cuenta con 3 bandas principales que se basan en imágenes del satélite Landsat y tienen una resolución de 30x30 m.
2.1.1.1 Información principal
La base de datos tiene un total de 13 bandas que incluyen la información relativa a la cobertura y pérdida en la cobertura arbórea, así como composiciones de imágenes satelitales. Las 3 bandas que se utilizaron en este trabajo son las siguientes:
- Cobertura arbórea en el año 2000
La cobertura arbórea corresponde al porcentaje de vegetación con una altura mayor a 5m de altura en un píxel de 30 x 30 m en el año 2000. Es decir, no incluye cultivos de baja altura como piña, hortalizas, arroz, entre otros. Pero si puede contemplar sembradíos que cumplan esa característica como cultivos de palma maduros y teca. En este trabajo se llama bosque a todas las áreas con un porcentaje de cobertura arbórea igual o mayor al 30%.
El año 2000 es el único año que cuenta con esta información y se considera como la línea base.
- Pérdida de cobertura arbórea
Corresponde a la pérdida total de cobertura boscosa en un área de 30 x 30 m. Se debe tomar en cuenta que esta pérdida puede ser causada por fenómenos naturales o antropogénicos tales como degradación, deforestación, incendios, deslizamientos y corta de cultivos mayores a 5 m de altura. En este trabajo se define esta pérdida como deforestación.
Una vez que se identifica una pérdida, esa área queda marcada como pérdida, aunque posteriormente fuera reforestada.
- Año de pérdida
Año en que se registró la primera pérdida de cobertura arbórea.
2.1.1.2 Características de la base de datos
Esta información cuenta con varias características que la diferencian de otras bases de datos de cambio en el uso del suelo que la hacen más útil y a la vez fomentan la reproducibilidad de este algoritmo en otras latitudes:
- Libre acceso y transparencia
La información es gratuita, además la metodología y medidas de exactitud están disponibles para análisis de terceros. Es posible descargar esta información desde varios medios. Para este trabajo se utilizó la versión 1.7 2019 descargada de la plataforma Google Earth Engine.
- Comparabilidad
La metodología aplicada no depende de fronteras administrativas. Esto permite hacer comparaciones entre áreas que anteriormente no podían cotejarse, ya que utilizaban diferentes definiciones de bosque o había dudas acerca su veracidad.
- Completitud
Por temas presupuestarios, normalmente estos estudios se concentran en áreas específicas, sin embargo, esta información está disponible prácticamente a nivel global. Esta característica hace que sea posible replicar esta metodología en otras partes del mundo.
- Resolución
La resolución de estas imágenes es de 30 x 30 m por píxel lo que hace posible detectar cambios en la cobertura de la tierra con un gran nivel de detalle en comparación con otras bases de datos de escala mundial.
2.1.1.3 Exactitud
Todo análisis de esta naturaleza tiene una incertidumbre asociada. Esto puede ocurrir tanto por la metodología como por la calidad y resolución de los datos.
Los autores de esta base de datos han publicado dos estudios respecto a su exactitud -(Hansen et al. 2013) y (Tyukavina et al. 2015)-. La publicación original clasificó un total de 1500 bloques de 120 m de lado a partir de imágenes Landsat, MODIS y Google Earth para compararlas con los mapas de pérdida y ganancia globales. El segundo estudio se enfocó en los trópicos y utilizó 3 000 puntos de control con la misma resolución que la base de datos.
El cuadro 2.1 muestra la exactitud encontrada en estos estudios:
| Estudio | Bioma/zona | Falsos positivos | Falsos negaticos |
|---|---|---|---|
| Hansen, 2013 | Tropical | 13% | 16.9% |
| Tyukavina, 2015 | Latinoamérica | 4% | 17% |
Recientemente Cunningham et al (Cunningham, Cunningham, and Fagan 2019) analizaron la exactitud de la clasificación de bosque o no-bosque de esta y otras 2 bases de datos, según el umbral del porcentaje de cobertura arbórea en Costa Rica. Los resultados obtenidos para el GFC Dataset son alentadores en el sentido que son comparables con productos realizados a nivel nacional y global para ciertos umbrales. Como ejemplo, si se considera como bosque aquel píxel que tiene una cobertura arbórea superior al 30%, la exactitud es cerca del 80% y es la mejor de los tres datasets analizados.
Los resultados de ese estudio indican que existen varios sesgos que hacen que la clasificación del GFC sea errónea principalmente en precipitaciones menores a los 189 mm anuales y en elevaciones mayores a los 2 000 m.s.n.m.. En el ACLAP, los valores de exactitud son muy similares a los encontrados a nivel nacional y por ende se descarta que haya un sesgo relevante en comparación con otras zonas como los bosques tropicales secos.
En este trabajo se definieron como bosques aquellos píxeles que cuentan con un porcentaje de cobertura arbórea mayor o igual a un 30%. Consecuentemente, la deforestación registrada en zonas con un menor porcentaje de cobertura no fue tomada en cuenta.
2.1.2 Otras capas de información
A continuación se describen las otras fuentes de información utilizadas como insumo para el modelo de aprendizaje automático supervisado.
| Capa | Descripción | Formato | Fuente | Año |
|---|---|---|---|---|
| Elevación | Modelo de elevación digital con resolución 30x30m | Ráster | United States Geological Survey | 2000 |
| Pendiente | Pendiente del terreno basado en la elevación | Ráster | United States Geological Survey | 2000 |
| Calles y carreteras | Ubicación de calles y carreteras | Polilíneas | Open Street Map | 2020 |
| Poblados | Ubicación de centros de poblados | Puntos | Atlas Digital CR2004 | 2004 |
| Áreas Silvestres Protegidas | Ubicación de áreas que están protegidas por el estado. | Polígonos | SINAC | 2020 |
| Corredores Biológicos | Ubicación de corredores biológicos. | Polígonos | SINAC | 2020 |
| Capacidad de Uso del Suelo | Capacidad natural para soportar distintas formas de uso. | Polígonos | Fundación Neotrópica | 1995 |
| Precio | Valor por metro cuadrado en zonas homogéneas. | Polígonos | ONT / Ministerio de Hacienda | ND |
| Índice de Desarrollo Social | IDS a nivel distrital | Polígonos | INEC | 2013 |
| Biotemperatura | Biotemperatura media anual de la zona de vida | Polígonos | Centro Científico Tropical | 2004 |
| Precipitación | Precipitación media anual de la zona de vida | Polígonos | Centro Científico Tropical | 2004 |
| Piso | Piso altitudinal de la zona de vida | Polígonos | Centro Científico Tropical | 2004 |
| Zona de vida | Nombre de la zona de vida | Polígonos | Centro Científico Tropical | 2004 |
2.2 Feature Engineering
Se realizaron diferentes modificaciones a las capas antes descritas para ajustarlas y crear variables que fueron de utilidad en el proceso de clasificación.
2.2.1 Transformaciones básicas
La unidad de análisis de este trabajo es el píxel o área de terreno de 30 x 30 m. En concordancia, las variables a utilizar deben tener un formato tal que pueda ser asignado a la unidad de estudio. Teniendo esto en cuenta, se transformó la información de tipo “polilínea” o “punto” por medio de cálculos basados en la proximidad. Tales operaciones se aplicaron para obtener las siguientes variables:
- Cercanía a calles y carreteras.
Para cada píxel en el dataset, se calculó la distancia más cercana de su centroide a una carretera y calle. Así se obtuvieron dos variables que representan la distancia mínima a una carretera y la distancia mínima a una calle para cada píxel.
- Cercanía a centros de poblados
Del mismo modo, para cada centroide de cada píxel se calculó la distancia mínima a un poblado y se creó una variable que representa la cercanía a un poblado.
2.2.2 Operaciones con rásters
Puesto que un píxel está rodeado por otros que pueden influir en su estado futuro, se deben realizar operaciones con rásters que permitan capturar el comportamiento y estado de los vecinos. Estas operaciones son llevadas a cabo para las variables que tienen un formato ráster y que están relacionadas con la cobertura o pérdida boscosa.
Se calcularon operaciones con rásters para cada año y en radios de influencia de 50, 100, 500 y 1000 m para abarcar comportamientos en diferentes grados de lejanía.
2.2.2.1 Cobertura boscosa promedio
Calcula el promedio de la cobertura boscosa de los píxeles cercanos al píxel en estudio. Esta variable junto con “cantidad de bosques” representa la salud del bosque y da información respecto al porcentaje de cobertura boscosa de los píxeles que son deforestados.
Figura 2.1: Cálculo de cobertura boscosa promedio.
2.2.2.2 Cantidad de bosques
Suma la cantidad de bosques - mayor o igual al 30% de cobertura boscosa- que se encuentran en su proximidad. Esta variable puede ser un indicador del fraccionamiento de los bosques y es la base para determinar el avance de la pérdida boscosa cada año.
Figura 2.2: Cálculo de cantidad de bosques.
2.2.2.3 Cantidad de píxeles deforestados
Suma la cantidad de píxeles con pérdida de cobertura boscosa en la proximidad del píxel. Esta variable muestra el estado de la deforestación año tras año.
Figura 2.3: Cálculo de cantidad de píxeles deforestados.
Es importante indicar que estas variables por sí solas no necesariamente aportan información relevante, sin embargo, al interactuar entre sí, se complementan y permiten una caracterización apropiada del problema de pérdida boscosa.
Para la aplicación de tales cálculos, se realizó un “buffer” de 2 km alrededor del área de estudio para calcular las operaciones con rásters en todas las celdas dentro de ACLAP sin tener valores faltantes -NA en las figuras-.
2.2.3 Temporales
Para cada una de las variables descritas en el apartado “operaciones con rásters” se calculó el cambio en el tiempo dentro de cada periodo de estudio -4 años-. La variación temporal es representada como la pendiente de la recta de mejor ajuste - \(\beta_1\)- entre la variable en estudio y los años:
\[Píxeles\quad deforestados = \beta_0 + \beta_1 * años\]
2.3 Periodos de estudio y sets de datos
Puesto que el análisis es de carácter temporal, se agrupó la información disponible en periodos. Los periodos de estudio se ajustan a la información disponible del GFC cuya pérdida boscosa fue registrada durante el 2001-2019. Para la realización de este análisis fue necesario contar con cuatro sets de datos que funcionan como entrenamiento, validación, prueba y predicción. Considerando este requerimiento y los años disponibles en la base de datos, la información fue dividida en cuatro períodos de cuatro años cada uno.
En el cuadro 2.3 se muestran las características más importantes de estas agrupaciones:
| Periodo | Set | Variables | Etiqueta deforestación |
|---|---|---|---|
| 1 | Entrenamiento - validación | 2004-2007 | 2011 |
| 2 | Entrenamiento - validación | 2008-2011 | 2015 |
| 3 | Prueba o testeo | 2012-2015 | 2019 |
| 4 | Predicción | 2016-2019 | Variable a predecir |
Es importante notar que una vez que la base de datos del GFC cuente con la pérdida boscosa del año 2020, los periodos podrán ser de 5 años. Esto ampliará el horizonte predictivo y aportará más información para el entrenamiento del modelo.
2.4 Algoritmo predictivo
Se utilizó un ensamble de aprendizaje automático de tipo “boosting”. Este tipo de ensambles entrenan de manera secuencial múltiples algoritmos normalmente con un alto grado de sesgo y poca varianza. En cada iteración el algoritmo se concentra en las observaciones que fueron más difíciles de clasificar en la iteración anterior lo que permite reducir el sesgo de los algoritmos individuales.
Específicamente se utilizó el algoritmo desarrollado por Microsoft en el 2017 denominado “Light Gradient Boosting Machine” conocido popularmente como LightGBM. Este algoritmo es un ensamble de tipo boosting que consiste en la aplicación de gradientes descendientes en árboles de decisión.
2.5 Modelado
Una gran parte del proceso de inteligencia artificial consiste encontrar un modelo con los hiperparámetros adecuados que generen un balance óptimo entre el sesgo y la varianza para evitar el sobreajuste a los datos de entrenamiento.
De acuerdo con lo anterior, se entrenó y validó el modelo con diferentes periodos para evitar utilizar datos de más de un periodo en el conjunto de entrenamiento ya que esto puede ser una fuente de sobreajuste. El procedimiento consiste en aplicar un conjunto de combinaciones de hiperparámetros y entrenar el modelo con el “Set 1”, validarlo con el “Set 2” y viceversa. La métrica utilizada para determinar la utilidad del modelo fue el área bajo la curva ROC -Receiver Operation Characteristic- ampliamente utilizada en algoritmos de clasificación desbalanceados.
Este proceso fue realizado reiteradas ocasiones hasta alcanzar el modelo que mejor se ajustó a los datos de validación. Se optimizaron los siguientes hiperparámetros que son parte del marco del paquete parsnip:
mtry: El número de predictores que serán aleatoriamente seleccionados en cada partición cuando se crea un modelo de árbol.
trees: El número de árboles contenidos en el ensamble.
min_n: La cantidad mínima de observaciones en un nodo que son requeridos para dividir el nodo una vez más.
tree_depth: Profundidad máxima del árbol.
learn_rate: La tasa a la que el algoritmo de boosting se adapta tras cada iteración.
loss_reduction: La reducción en la función de pérdida para dividir el nodo una vez más.
sample_size: Cantidad de datos expuestos a la rutina de entrenamiento.
Una vez obtenidos los hiperparámetros óptimos, el algoritmo se entrenó con el Set 1 y Set 2, posteriormente se aplicó en el Set 3 -test set- para conocer los resultados considerando datos nunca vistos por el modelo. Finalmente, el modelo se entrenó con todos los datos disponibles para predecir la pérdida de cobertura boscosa en el año 2023.
3 Análisis exploratorios
El Área de Conservación La Amistad Pacífico tiene un total de 6.18 millones de píxeles y por ende el set de datos inicial tiene esa misma cantidad de observaciones. Esta cantidad de píxeles fue utilizada para computar las operaciones con rásters, sin embargo, las áreas que no son bosques en el periodo de estudio no son consideradas en el entrenamiento pues no aportan información al modelo. Es por esta razón que la cantidad de filas del set de datos se reduce conforme avanzan los periodos.
Por otro lado, los píxeles que experimentan una pérdida de cobertura arbórea son muy inferiores al total de píxeles con bosque. En el cuadro 3.1, se muestran los píxeles deforestados por periodo:
| Periodo | Total píxeles | Píxeles deforestados | P. deforestados (%) |
|---|---|---|---|
| 1 | 4 872 287 | 108 275 | 2.22 |
| 2 | 4 764 012 | 39 198 | 0.82 |
| 3 | 4 724 814 | 55 001 | 1.16 |
| 4 | 4 669 813 | NA | NA |
3.1 Variables utilizadas
Cada uno de los píxeles en el set de datos tiene asociado un total de 95 características. Esto hace que para cada periodo de análisis la cantidad de datos sea de aproximadamente 4.76 millones para un total de 19.03 millones en todos los periodos. Las 95 variables se describen a continuación:
- label: Indica si el píxel experimentará deforestación en los próximos 4 años. Variable binaria.
- period: Indica el periodo correspondiente al set de datos. Variable numérica.
- x, y: Coordenadas geográficas CRTM05. Variable numérica.
- tree_cover: Porcentaje de cobertura arbórea en el año 2000. Variable numérica.
- slope: Pendiente del terreno. Variable numérica.
- dem: Elevación del terreno. Variable numérica.
- asp: Nombre del Área Silvestre Protegida. Variable categórica.
- is_asp: Variable binaria que indica la presencia de un ASP. Variable binaria.
- corredor: Nombre del Corredor Biológico. Variable categórica.
- is_corredor: Variable binaria que indica la presencia de un corredor biológico. Variable binaria.
- canton: Nombre del distrito. Variable categórica.
- ids: Índice de desarrollo social distrital. Variable numérica.
- precio: Precio por metro cuadrado por zona homogénea. Variable numérica.
- z_vida: Nombre de la zona de vida de Holdridge. Variable categórica
- precip = Promedio de precipitación anual según zona de vida en mm. Variable numérica.
- piso: Piso de la zona de vida: basal, premontano, montano bajo, montano, subalpino. Enumerados del 1 al 5. Variable numérica.
- biotemp: Biotemperatura de la zona de vida en grados centígrados. Variable numérica.
- cap_uso: Capacidad de uso del suelo. Variable numérica.
- camino_dist: Distancia más corta del píxel a un camino en metros. Variable numérica.
- calle_dist: Distancia más corta del píxel a una calle en metros. Variable numérica.
- poblados_dist: Distancia más corta del píxel a un poblado en metros. Variable numérica.
- 16 variables del tipo “loss_count”: Conteo de pérdida boscosa en las cercanías del píxel. Existe una para cada año dentro de los periodos de estudio y para radios de 50, 100, 500 y 1000 m. t_0 indica el último año del periodo, t_1 representa un año antes del último año del periodo y así sucesivamente. Variable numérica.
- 4 variables tipo “sum_loss”: Suma de las pérdidas de cobertura boscosa de cada periodo para todos los radios de cercanía.
- 16 variables del tipo “forest_count”: Conteo de celdas con cobertura arbórea superiores al 30% en las cercanías del píxel. Existe una para cada año dentro de los periodos de estudio y para radios de 50, 100, 500 y 1000 m. t_0 indica el último año del periodo, t_1 representa un año antes del último año del periodo y así sucesivamente. Variable numérica.
- 4 variables tipo “mean_forest”: Promedio de los píxeles con cobertura boscosa de cada periodo para todos los radios de cercanía.
- 16 variables del tipo “cover_mean”: Promedio de cobertura boscosa en las cercanías del píxel. Existe una para cada año dentro de los periodos de estudio y para radios de 50, 100, 500 y 1000 m. t_0 indica el último año del periodo, t_1 representa un año antes del último año del periodo y así sucesivamente. Variable numérica.
- 4 variables tipo “mean_cover”: Promedio de la cobertura boscosa de cada periodo para todos los radios de cercanía.
- 4 variables del tipo loss_trend: Calcula la tendencia anual de pérdida boscosa en las cercanías del píxel. Existe una tendencia para cada uno de los 4 radios de cercanía.
- 4 variables del tipo forest_trend: Calcula la tendencia anual de la presencia de bosques en las cercanías del píxel. Existe una tendencia para cada uno de los 4 radios de cercanía.
- 4 variables del tipo cover_trend: Calcula la tendencia anual del promedio de cobertura arbórea en las cercanías del píxel. Existe una tendencia para cada uno de los 4 radios de cercanía.
3.2 Distribución de las variables
A continuación, se muestran varios gráficos de tipo “boxplot” que detallan la distribución de cada variable con respecto a su estado futuro - deforestación o no deforestación-.
Figura 3.1: Distribución de variables I
Algunas observaciones de las variables expuestas son:
- Las áreas cercanas a calles, caminos y poblados son más propensas a ser deforestados que los que se encuentran más alejados.
- La tierra con menor capacidad de uso es más propensa a deforestación probablemente por su vocación agrícola o urbana.
- Las pendientes y elevaciones bajas son más propensas a deforestación, esto es evidente puesto que los asentamientos más importantes se encuentran en la zona más plana y baja del ACLAP. Esto es captado especialmente en la variable “piso” que indica que las zonas más bajas tienen una mayor propensión a la deforestación que las zonas que se encuentran en pisos altidudinales mayores.
- Los terrenos en zonas más caras son más propensos a deforestación que los que no.
- Las áreas con una mayor cobertura boscosa son más propensos a sufrir una pérdida boscosa.
La figura 3.2 algunas variables que toman en cuenta las áreas cercanas al píxel en estudio.
Figura 3.2: Distribución de variables II
Como se puede observar en la figura anterior, cuando hay una menor cantidad de bosque en las cercanías, el píxel es más propenso a ser deforestado. En la misma línea, es más probable que un área sea deforestada si ha habido pérdidas boscosas en su cercanía.
La figura 3.3 muestra los valores agregados de las operaciones con rásters de los cuatros años de cada periodo.
Figura 3.3: Distribución de variables III
Finalmente se muestran algunas de las variables que representan cambios en el tiempo.
Figura 3.4: Distribución de variables IV
Es más probable que el píxel sea deforestado si existe una tendencia a la baja en cuanto a la cobertura boscosa y cantidad de bosques.
4 Resultados
Luego de realizar varias iteraciones con diferentes metodologías -grid search y optimización bayesiana- que buscan en el espacio de hiperparámetros la mejor combinación, se obtuvo el modelo con mejores resultados en el conjunto de validación.
El resultado de las predicciones de la probabilidad de deforestación en el año 2023 se muestra en la figura 4.1:
Figura 4.1: Probabilidad de cada píxel a ser deforestado en el 2023.
Descargue este mapa en formato ráster aquí.
En el mapa anterior se pueden observar que las zonas en rojo o rojizas son las más probables a estar deforestadas en el 2023, por el contrario, las áreas azules tienen una mayor probabilidad de mantener su cobertura boscosa. Las áreas que no son boscosas o se encuentran fueran del área de estudio son representadas en gris.
4.1 Métricas
El modelo alcanzó un valor de ROC-AUC de 0.827 en los conjuntos de validación y un 0.724 en el conjunto de prueba. La diferencia en este valor puede deberse a: variaciones en el comportamiento de la pérdida boscosa en ACLAP, al hecho de que la metodología del GFC fue actualizada a partir del 2011 -ver referencia- que coincide con el periodo de testeo, o a un sobreajuste en el conjunto de entrenamiento. Es de esperar que una vez que los datos de GFC cuenten con una metodología uniforme, los resultados experimenten una mejoría.
En la figura 4.2 se aprecia como el modelo de inteligencia artificial clasifica las áreas a ser deforestadas en el set de datos de prueba.
Figura 4.2: Distribución de las predicciones según etiqueta.
En la figura anterior es posible apreciar como los píxeles que efectivamente experimentaron una pérdida boscosa tienden a tener una mayor probabilidad a ser deforestados. Como es de esperar, la clasificación es difusa y por tanto el algoritmo incluirá una importante cantidad de píxeles sin pérdida boscosa independientemente de la probabilidad de corte elegida.
En muchos casos no es conveniente describir las predicciones de manera númerica -probabilidades- si no que es mejor especificarla categóricamente como “deforestación” o “no deforestación”. Esto se puede alcanzar al definir una probabilidad de corte que discrimine una categoría de la otra. Para elegir una probabilidad de corte adecuada es necesario tener claro los conceptos de sensibilidad y especificidad:
\[Sensibilidad = verdaderos\:positivos/(verdaderos\: positivos+falsos\: negativos)\] \[Especificidad = verdaderos\:negativos/(verdaderos\: negativos+falsos\: positivos)\]
Seguidamente, en la figura 4.3 se muestran tales conceptos para diferentes probabilidades de ser deforestado en el set de prueba.
Figura 4.3: Sensibilidad y especificidad para diversas probabilidades de ser deforestado.
Esta figura, muestra el “trade-off” entre la sensibilidad y especificidad. Si se fija una probabilidad de corte baja, se clasificará la mayoría del área de ACLAP como futura deforestación y por ende incluirá la gran mayoría de píxeles de deforestación -verdaderos positivos- pero cometerá muchos errores en los píxeles que no se deforestarán -falsos positivos-. En caso contrario, una alta probabilidad de corte implicará una baja sensibilidad pero una alta especificidad. Según la aplicación particular, la persona usuaria se decantará por alguna probabilidad de corte considerando su impacto en la clasificación.
4.2 Contribución de variables
Los modelos complejos de inteligencia artificial son difíciles de interpretar más allá de la respuesta -output- que brindan. En el año 2017, se desarrolló una metodología denominada “Shapley Additive Explanations” (Lundberg and Lee 2017) que aporta un marco que facilita el entendimiento de las predicciones del modelo.
Con esta metodología es posible cuantificar la contribución de cada variable relativa a la predicción promedio de los datos. En este caso, es interesante conocer las características que brindan más información predictiva al modelo para definir si un píxel es deforestado o no.
La figura 4.4 muestra estos valores calculados para una muestra de alrededor de 2.4 millones de píxeles del set de entrenamiento y validación.
Figura 4.4: Variables más importantes del modelo según el promedio del valor absoluto de su contribución.
De la figura anterior se puede notar que:
La cobertura boscosa en el año 2000 es la variable que aporta más información para la clasificación.
Cuatro de las seis variables más importantes tienen que ver con la tendencia del bosque y cobertura boscosa a los alrededores del píxel en cuestión.
Es notorio que los efectos locales con distancias inferiores a los 500 m no tienen un alto impacto en el modelo. Al ser una zona con muy poca deforestación, es probable que las áreas más cercanas no cuenten con una cantidad suficiente de pérdida boscosa para generar un efecto relevante.
La elevación -dem- es importante en un marco en donde la actividad antrópica se concentra en las áreas bajas.
Se debe recalcar que la contribución estimada es de una muestra que contiene todas las observaciones con etiqueta “deforestación” y 15 veces más con la etiqueta “no deforestación” -técnica de “undersampling”-. Esta condición sumada a que muchas de las variables no sean independientes hace que las estimaciones de las contribuciones no sean del todo exactas.
Ahora bien, es interesante conocer la naturaleza de los aportes predictivos de algunas variables. Esto se puede realizar al graficar las contribuciones de las variables y sus respectivos valores en cada píxel. A continuación, se presenta tal información de las tres variables con la mayores aportes -en valor absoluto- promedio al modelo.
Figura 4.5: Constribución en la probabilidad de ser deforestado de cada píxel de las tres variables más relevantes.
Las contribuciones individuales de las variables muestran como las elevaciones más altas -dem-, la tendencia a la baja de bosques a 1000 m a la redonda -forest_trend_1000- y una cobertura boscosa alta en el año 2000 -tree_cover- aportan a favor de la deforestación y viceversa.
Resulta interesante que las variables que representan la tendencia en la cobertura boscosa más cerca del píxel en cuestión - 50 m, 100 m y 500 m - fueron relegadas por aquellas calculadas a una mayor distancia -1000 m-. Esto puede deberse a que las tendencias de pérdida boscosa locales tienen una variabilidad mayor, producto de la alta dispersión de los píxeles con deforestación en el espacio. Es de esperar que en un contexto donde la pérdida boscosa sea más abundante, las variables que se encuentran más cerca de la unidad de estudio tengan un peso mayor y de esa forma permitan una clasificación más localizada.
4.3 Desagregación del modelo
Con el fin de determinar si el desempeño del modelo realizado en ACLAP tiene un comportamiento homogéneo en toda el área, se calculó la métrica ROC-AUC para tres zonas distintas asociadas a los pisos definidas en las zonas de vida (Bolaños, Watson, and Tosi 2005) de la región y también para los cantones de la zona.
Las métricas correspondientes a cada uno de los pisos de las zonas de vida son:
| Piso | Área bajo la curva ROC |
|---|---|
| Basal | 0.56 |
| Premontano | 0.64 |
| Montano bajo, montano y subalpino | 0.82 |
Las métricas para los cantones son:
| Cantones | Área bajo la curva ROC |
|---|---|
| Coto Brus | 0.74 |
| Buenos Aires | 0.68 |
| Pérez Zeledón | 0.70 |
| otros | 0.63 |
Al analizar ambos cuadros se hace evidente que el modelo tiene un rendimiento con alta dispersión según los pisos de las zonas de vida -asociados con elementos topográficos y del tipo de vegetación que los cubre-. En ACLAP, el área urbana está concentrada en los primeros dos pisos mientras que en los pisos superiores el uso de suelo predominante es el boscoso.
Los cantones cuentan con una topografía más heterogénea que atraviesa los diferentes pisos de las zonas de vida y diversos usos de la tierra. Esto resulta en métricas con menor dispersión que en el caso de los pisos de las zonas de vida.
En la figura 4.6 , se muestra la ubicación de las zonas de vidas y los cantones.
Figura 4.6: Ubicación de zonas de vidas y cantones.
4.4 Discretización de los resultados
Según la aplicación específica, se podrá fijar un nivel de sensibilidad y especificidad acorde. Para efectos de este proyecto, se elegirá una probabilidad de corte de 0.055 que permite incluir una importante cantidad de píxeles deforestados sin añadir demasiados falsos positivos.
Considerando esta probabilidad de corte, se pueden predecir correctamente un total de 22.4% de los píxeles que sufrirán deforestación y se incluyen como falsos positivos un total de 10.2% de todos los píxeles que mantendrán su cobertura boscosa. En el cuadro 4.3 y en la figura 4.7 se presenta la matriz de confusión en términos de píxeles.
| Pred./Verdad | Verdad deforestación | Verdad no-deforestación |
|---|---|---|
| Pred. deforestación | 12 314 | 473 989 |
| Pred. no-deforestación | 42 687 | 4 195 824 |
Figura 4.7: Matriz de confusión representada espacialmente.
Teniendo en cuenta la probabilidad de corte elegida, se discretizaron los resultados de la predicción de deforestación en el año 2023. Esto permite generar la figura 4.8 que muestra las sitios con un mayor riesgo -probabilidad mayor a la probabilidad de corte- y menor riesgo a ser deforestados -probabilidad menor a la probabilidad de corte- en el futuro cercano.