Predicción de zonas de desertificación usando IA.

Zaragoza. Primera Edición. 2022

La desertificación se produce por la degradación de la tierra que se da en zonas secas: áridas, semiáridas y subhúmedas secas. La actividad humana y las variaciones climáticas están entre las causas de esta degradación del suelo, que impacta en los ecosistemas y en los recursos y modos de vida de los habitantes de las zonas afectadas. Estas zonas secas son aquellas con un índice de aridez (IA=P/ETP) inferior a 0,65.

España es uno de los territorios que sufre esta problemática. En Aragón el 75% del territorio está en riesgo de desertización. El cambio climático y el incremento de las temperaturas amenazan las tierras semiáridas que predominan en el Valle del Ebro, desde la ciudad de Huesca hacia el sur. Identificar a tiempo estas zonas será de gran ayuda para el apoyo a la toma de acciones que permitirán mitigar el impacto que la desertificación trae sobre la seguridad alimentaria así como la reducción de la pérdida de hábitats naturales y biodiversidad, todo esto, alineado con los Objetivos de Desarrollo Sostenible de Naciones Unidas.

Enfoque sobre los objetivos de Desarrollo Sostenibles 13 y 15. Accion por el clima y biodiversidad

Enfoque sobre los objetivos de Desarrollo Sostenibles 13 y 15. Fuente: https://sdgs.un.org/goals

El proyecto explora el uso de imágenes satelitales para observar cambios en la vegetación de Aragón y estimar valores de NDVI. Es un primer acercamiento al estudio de la degradación del suelo, con posibles aplicaciones en reforestación y gestión de cultivos; no equivale por sí solo a predecir la desertificación. Con el apoyo de Saturdays.AI Zaragoza, cuyo objetivo es acercar la IA a todos los ciudadanos y aplicarla a problemas sociales, queremos dar solución o una aproximación a resolver este escenario que nos hemos planteado como equipo de trabajo.

Nos planteamos tres preguntas: ¿cómo podemos observar cambios en la vegetación?, ¿qué imágenes sirven para estudiarlos?, ¿cómo puede ayudarnos la IA a estimar sus valores? A lo largo del artículo explicamos nuestro enfoque.

Imágenes Satelitales

Las herramientas de teledetección hacen posible la visualización y cálculo de cambios que ocurren sobre la superficie terrestre, aquellas que nos interesan, los que corresponden a la teledetección pasiva (Definición de Teledetección), en nuestro caso concreto, captan la radiación solar reflejada por la superficie. Es bien sabido que organismos vivos reflejan y absorben radiación solar en forma diferente a como lo hacen los objetos inertes como calles o edificios. Aprovechando este principio, los satélites generan imágenes que contienen estas diferencias en forma de Niveles Digitales y bandas, que se traducen en unidades numéricas que se corresponden con los valores de longitudes de onda de las diferentes secciones del espectro electromagnético (Espectro electromagnético y teledetección). En resumen, miden la radiación reflejada por objetos como edificios o autopistas o entes vivos como bosques o vegetación en general.

Principales Componentes de un sistema de Teledetección

Principales Componentes de un sistema de Teledetección. Fuente : Curso de teledección https://static.uvq.edu.ar/mdm/teledeteccion/unidad-1.html

Índice de vegetación NDVI

Los cambios en la vegetación con el soporte de las imágenes satelitales se pueden medir de diferentes formas, una de ellas es extrayendo índices de vegetación, “Los Índices de Vegetación son combinaciones de las bandas espectrales registradas por los satélites de Teledetección” (¿Cuál es el índice de Vegetación NDVI?) que permiten identificar vegetación sobre la superficie terrestre. Es allí,con ayuda de estos índices en donde se resaltan y se pueden analizar numéricamente las propiedades de la vegetación como tal.

Interacción de los elementos de la superficie terrestre con la radiación

Interacción de los elementos de la superficie terrestre con la radiación. Fuente: Universidad de Murcia. Apuntes de teledetección. <https://www.um.es/geograf/sigmur/teledet/tema02.pdf>

NDVI significa índice de vegetación de diferencia normalizada. Se calcula a partir de la reflectancia registrada en las bandas roja e infrarroja cercana y permite estudiar la presencia y el estado de la vegetación en imágenes satelitales.

Los patrones del espectro electromagnético para identificar la vegetación muestran que la vegetación absorbe la sección de radiación solar que corresponde al rojo. A su vez se reflejan valores elevados en el infrarrojo cercano y podremos identificar valores multiespectrales brillantes en esta banda ( Valores Espectrales del NDVI).

Representación de cobertura vegetal en RED y NIR

Representación de cobertura vegetal en RED y NIR. Fuente: Calculo de NDVI http://www.gisandbeers.com/calculo-del-indice-ndvi/

Así podemos diferenciar cubiertas vegetales en diferentes estados. Nuestro índice podrá ser obtenido mediante la siguiente relación:

NDVI = (Banda infrarroja cercana — Banda roja) / (Banda infrarroja cercana + Banda roja)

O lo que es lo mismo, y en términos de nomenclatura cuando trabajamos las bandas multiespectrales:

NDVI = (NIR — RED) / (NIR + RED)

Una vez se ha calculado el índice para cada uno de los píxeles de la imagen, los valores posibles oscilan entre -1 y 1.

Valores de NDVI para diferentes zonas

Valores de NDVI para diferentes zonas. Fuente: propia

En general, los valores positivos más altos se asocian con vegetación más densa o vigorosa. Su interpretación depende también del tipo de cubierta y de las condiciones de la imagen.

Manipulando la imagen

Para comprender mejor cómo podemos usar estos valores entre bandas y niveles digitales, es importante identificar las propiedades de una imagen en general y cómo se ven estos valores a nivel de computación.

Derecha: como un ser humano ve una imagen. Izquierda: como la procesa un ordenador (matrices)

Derecha: como un ser humano ve una imagen. Izquierda: como la procesa un ordenador. Fuente: Material Saturdays.AI, Zaragoza. Edición I. 2021–2022. Tratamiento de Imágenes.

Un ordenador representa una imagen mediante una matriz. Cada posición (i, j) identifica un píxel y contiene uno o varios valores que describen su intensidad o color.

En las imágenes satelitales, cada píxel contiene medidas relacionadas con la radiación registrada por el sensor en las distintas bandas espectrales.

Firmas espectrales en imágenes.

Firmas espectrales en imágenes. Fuente: Imagen 45. Utilización en imágenes multiespectrales. https://static.uvq.edu.ar/mdm/teledeteccion/unidad-1.html

Obtención y preparación de datos

Obtención

Para estudiar los cambios en la superficie terrestre utilizamos imágenes del satélite Sentinel-2A. En el flujo de trabajo original empleamos la biblioteca de Python SentinelSat para la descarga de imágenes en el área que nos interesa. Es de anotar que las imágenes se han tomado con los siguientes condicionantes: para un área específica detallada en un fichero .json, para la misma tesela (zona de captura de imagen del satélite) y sección de órbita para garantizar que se está tomando la misma zona de estudio para diferentes momentos de tiempo.

Fuente: propia (Codigo completo en https://github.com/desertificacion-AI/desertificacion-AI)

Fragmento de código para la descarga de imágenes sentinel del sitio scihub.copernicus.eu Fuente: Propia

NOTA: Para la obtención y captura de los datos se han descargado las imágenes de Sentinel con un porcentaje de nubes máximo de 10%.

Una vez ejecutado el script para la descarga de imágenes del satélite, se descarga una serie de ficheros zip que contienen cada uno la siguiente estructura al momento de descomprimirlos:

Estructura de fichero descargado del sitio scihub.copernicus.eu Fuente: propia

Cada una de las imágenes Sentinel vienen dadas con un consecutivo de nombre de imagen con la siguiente estructura:

MS2_MSILLL_YYYYMMDD

Donde:

MS2: identifica la misión de Sentinel 2, pudiendo encontrar S2A (para Sentinel 2A) o S2B (para Sentinel 2B)
MSI: indica el instrumento de operación (MultiSpectral instrument)
LLL: indica el nivel de procesado del producto pudiendo encontrar los niveles L0, L1C, L1B o L2A
YYYY: designa el momento temporal UTC (año) en el que fue tomada la imagen
MM: designa el momento temporal UTC (mes) en el que fue tomada la imagen
DD: designa el momento temporal UTC (día) en el que fue tomada la imagen

Las descargas de Sentinel para el proyecto fueron:

Relación de imágenes. Fuente: Propia

Preparación

Una vez obtenidas las imágenes, sobre ellas calcularemos los índices NDVI mediante Script en Python.

Fuente: propia (Codigo completo en https://github.com/desertificacion-AI/desertificacion-AI)

Fragmento de código para el cálculo del NDVI para las imágenes descargadas

Para entender la salida de este script, las imágenes con el NDVI calculado (nuevos ficheros tif), se debe tener en cuenta que esta transformación imágenes nos da lugar a una imagen compuesta por una matriz de valores donde cada valor corresponde a un pixel (representación en la imagen) y este valor oscila entre -1 y 1. Recordando el apartado de Índice de Vegetación, a partir de 0.3 entendemos que hay vegetación para ese píxel (valor de matriz).

Captura de pantalla de las imágenes NDVI generadas por el script. Fuente: Propia

Detalle de una de las imágenes NDVI generadas por el script de Python. Fuente: Propia

¿Y la predicción de los valores de NDVI?

Ya estamos en el punto en donde, hemos identificado la zona sobre la que trabajaremos, hemos hallado los índices de vegetación para estas zonas en diferentes momentos de tiempo, ahora continuaremos con el establecimiento de una relación entre los valores medidos y cómo estos cambian en el tiempo (nuestra variable continua) para anticipar un conjunto de nuevos valores.

Deseamos construir un modelo que describa la relación entre las variables de estudio, con el fin, principalmente, de predecir los valores de una variable a partir de los valores de la otra. Elegiremos Regresor de bosque aleatorio porque queremos estimar valores numéricos continuos de NDVI, no categorías. Este algoritmo combina las predicciones de varios árboles de decisión y promedia sus resultados (material de Saturdays.AI Zaragoza, edición I, 2021–2022, sesión de regresores).

Fuente: Definición de Random Forest regressor. Material Saturdays AI. Zaragoza. Edición I. Sesión Regresores. 2021–2022

Breve resumen del funcionamiento de Random Forest Regressor

Para este modelo supervisado reservamos una zona de las imágenes para el entrenamiento y otra zona distinta para evaluar las predicciones de NDVI.

Evaluación

Entrenamiento

Para el entrenamiento del algoritmo, se tomaron 5 imágenes con NDVI calculado y una imagen adicional como salida conocida. De cada una de las imágenes se ha elegido un subset de igual tamaño para añadirlos a un vector en Python. Lo anterior se realiza por capacidad de procesamiento de las máquinas disponibles, las imágenes originales tienen dimensiones de 10.000×10.000 píxeles que no es posible procesar de forma ágil y eficiente con los recursos disponibles en el grupo de trabajo.

Imágenes con NDVI calculado. Fuente: propia.

Nota Aclaratoria: Cada imagen de NDVI calculado tiene un tamaño de 470MB aproximadamente, por lo que 6 de ellas nos da un estimado de 2,8GB para hacer el procesamiento de este dataset.

Fragmento de código para la elección del subset de datos para el Random Forest de entrenamiento. Fuente: propia

Con ayuda de la biblioteca sklearn preparamos los datos de entrenamiento y ejecutamos después la fase de prueba.

Prueba

Para probar el algoritmo, se tomaron las mismas 5 imágenes con el NDVI calculado pero esta vez, el subset elegido corresponde a otra porción de la imagen. En el entrenamiento hemos elegido una sección, para la prueba (test) elegimos una sección distinta.

Fragmento de código para la elección de los datos de test para el algoritmo. Fuente: Propia

El margen de error para la imagen predicha es el siguiente:

Margen de error resultante de la predicción. Fuente Propia.

Margen de error resultante de la predicción. Fuente: Propia.

Es necesario saber si el algoritmo es capaz de predecir correctamente una imagen. Utilizaremos una imagen que ya se tiene para validar la efectividad del algoritmo a la hora de predecir el resultado.

Derecha: Imagen original. Izquierda: Imagen predicha. Fuente: Propia.

Visualización en Streamlit

En la página desertIAragon se encuentra desplegada la aplicación de Streamlit para visualizar los cambios en NDVI calculados

Captura de pantalla del sitio web de Streamlit para mostrar imágenes calculadas. Fuente : Propia

En la sección Imágenes Satelitales, encontramos un subset de 2000×2000 píxeles de las imágenes extraídas del sitio web de Sentinel:

Fragmentos de imágenes satelitales. Fuente: elaboración propia.

Imágenes satelitales de 2000 x 2000 px. Fuente: propia

La sección de cuadrículas de estudio, nos permite elegir la coordenada deseada para mostrar los cambios en el NDVI calculado y sus valores. Elegir la coordenada en las opciones disponibles en el panel izquierdo:

Elección de coordenada para la cuadrícula de estudio. Fuente: propia

Valor de NDVI calculado para la imagen 1. Fuente: Propia.

Predicción

El sistema se entrena con una porción igual de cada una de las imágenes 1 a 5. Esta porción es un recuadro de coordenadas aleatorias y dimensiones igual al ancho/alto de la cuadrícula indicado en la banda lateral (en Streamlit).

La razón de elegir un recuadro reducido es el elevado coste computacional que tiene el entrenamiento del sistema. Este coste puede estar rondando aproximadamente 5 minutos de procesamiento para un fragmento inferior a 250 píxeles de dimensión, en un ordenador con una RAM de 8GB y un procesador: Intel Core I5 2,67 Ghz.

De esta forma la aplicación puede mostrar unos resultados de una forma relativamente ágil. Para el test se ha elegido las imágenes de las cuadrículas 1 a 5 con las que se ha obtenido el índice NDVI.

La predicción se hace con la cuadrícula 6. Se compara la imagen original con la predicha por la IA.

Comparación fragmento de la imagen original vs. imagen predicha por la IA

Datos generales de las imágenes:

Coordenadas del fragmento mostrado: desde (0, 0) hasta (5, 5) píxeles.

Error cuadrático promedio: 0.03719

Conclusiones

Existe software propietario como ERDAS para facilitar el cálculo de NDVI y algoritmos de regresión para predecir comportamientos, sin embargo, alineados con los principios de Saturdays.AI: acercar la IA a todos los ciudadanos y aplicarla a problemas sociales, hemos optado por una solución basada en herramientas de código abierto y disponibles para todo público.
Este proyecto tiene la potencialidad de explorar como una fase siguiente la capacidad de procesamiento de datos para imágenes de dimensiones más grandes, ya sea en resolución espacial o temporal, de forma más eficiente en la nube. Tomando ventajas como procesamiento en paralelo o distribuido se podrán añadir más entradas al algoritmo y procesamientos más complejos que no se pueden suplir con los recursos actuales
En esta primera evaluación, el modelo estimó los valores de NDVI de la zona de prueba con un error cuadrático medio de 0,03719. El resultado es útil para explorar cambios en la vegetación, pero no demuestra por sí solo que podamos predecir o evitar la desertificación. Harían falta más datos, evaluación temporal y otros indicadores de degradación del suelo antes de utilizarlo para tomar decisiones sobre el territorio.