23 Jul 2020 Teledetección y Python para analizar los efectos de la COVID-19
La reducción de la movilidad y la actividad económica durante la pandemia de la COVID-19 ha provocado también una reducción significativa de las emisiones de NO2 a la atmósfera. Utilizando técnicas de teledetección y algunas librerías en Python, podemos analizar los efectos de la COVID-19 en base a la reducción de estos gases contaminantes.
En este post mostraremos como programar algunas rutinas para llevar a cabo este tipo de análisis.
Datos de partida
El satélite Sentinel-5p de la Agencia Espacial Europea lleva a bordo el sensor TROPOMI. Este sensor permite monitorizar concentraciones de gases como el ozono, el metano, el monóxido de carbono, el dióxido de nitrógeno o el óxido de azufre.
En esta ocasión nos centraremos en el dióxido de nitrógeno (NO2). Se trata de un gas generado fundamentalmente por la actividad humana, como resultado de la combustión de los vehículos y determinados procesos industriales.
Con una resolución espacial de unos 7kmx7km, Sentinel-5p cubre diariamente toda la superfície terrestre. De este modo puede hacerse un seguimiento exhaustivo de la evolución de estos gases, y analizar episodios muy concretos como el producido durante la pandemia.
A través de esta aplicación (https://maps.s5p-pal.com/) desarrollada por la ESA y en base a los datos de Sentinel-5p puede visualizarse como evoluciona la concentración de NO2 en todo el planeta.
Acceso a los datos
En un post anterior (Copernicus, el programa de observación de la tierra) ya hicimos referencia a las diferentes formas de acceder a los datos del programa Copernicus. En esta ocasión utilizaremos la API de Copernicus Open Acces Hub junto con la librería Python Sentinelsat para obtener las imágenes que necesitamos.
El acceso a las imágenes Sentinel es gratuito, aunque es necesario disponer de una cuenta de usuario. Sin embargo –por ahora– las imágenes del Sentinel-5p están disponibles a través de una cuenta genérica, por lo que no es necesario darse de alta al portal Copernicus.
Preparación del entorno
Tal y como ya hemos comentado, para la búsqueda y obtención de las imágenes utilizaremos la librería Sentinelsat. El procesado y la visualización lo haremos utilizando las librerías HARP y VISAN respectivamente. Ambas forman parte de Atmospheric Toolbox, un paquete de librerías para facilitar a científicos el consumo, procesado y análisis de datos atmosféricos obtenidos con técnicas de teledetección.
La instalación de estas librerias la haremos en un entorno virtual de Conda. En este enlace tenéis las instrucciones para realizar la instalación de Conda en un sistema Windows y en un Linux.
Una vez instalado Conda, podéis crear un entorno virtual con la instrucción: conda create –name my_env python=3
A continuación lo podréis activar con la instrucción: conda activate my_env
Con el entorno activo, podéis instalar las librerías. Primero sentinelsat (pip install sentinelsat) y a continuación VISAN (conda install -c stcorp visan). Al instalar VISAN se realizará también la instalación de HARP.
Descarga de datos con Sentinelsat
Sentinelsat es una librería que facilita la búsqueda, filtrado y descarga de las imágenes del programa Copernicus disponibles en el repositorio de Copernicus Open Acces Hub. La conexión con el repositorio de imágenes Sentinel 5p a partir de la cuenta genérica se hará del siguiente modo:
api = SentinelAPI('s5pguest', 's5pguest', 'https://s5phub.copernicus.eu/dhus')Para la búsqueda de imágenes de una determinada zona y en un período determinado, utilizaremos el siguiente código:
products = api.query(footprint,
date=(date(2020, 2, 2), date(2020, 5, 8)),
producttype='L2__NO2___',
platformname='Sentinel-5')La variable footprint hace referencia al bbox de la zona a analizar, que podemos obtener a partir de un archivo en formato geoJSON:
# indicamos el bbox de la zona de descarga
footprint = geojson_to_wkt(read_geojson('bbox_spain.geojson'))En este caso, la variable products contendrá la referencia a todos los productos (escenas) capturadas por el Sentinel5p sobre el territorio español entre el 2 de febrero y el 8 de mayo de 2020.
La descarga de datos se ejecutará con la siguiente instrucción:
api.download_all(products)
Sentinel 5p tiene un periodo diario de revisita. Es decir, disponemos de imágenes de un mismo lugar del planeta cada dia. Y toda la superfície terrestre queda cubierta por 14 productos.
Se necesitan 2 productos para cartografiar todo el territorio español con imágenes Sentinel 5p. De este modo, y según los parámetros del código Python que hemos visto, se realizará la descarga de 192 productos (2 productos cada día, del 25 de febrero al 8 de mayo).
Procesar las imágenes Sentinel-5p con HARP
HARP es un kit de herramientas para leer, procesar y comparar imágenes de satélite.
Una vez instalado HARP en el entorno virtual Conda, podemos ejecutar sus comandos a través de la terminal, pero también mediante código Python.
En este caso utilizaremos HARP para crear un mosaico con los productos diarios que hemos descargado a través de Sentinelsat, animar temporalmente dichas imágenes, aplicar una máscara y filtrar los datos.
Crear un mosaico
Lo primero que tendremos que hacer es importar los productos de uno de los días que hemos descargado con Sentinelsat y que se encuentran por defecto en NetCDF. La importación la llevaremos a cabo con la función import_product de HARP, que además de facilitar la importación, permite añadir parámetros para acotar la extensión de la imagen, seleccionar la variable a representar, aplicar filtros, etc.
no2_escena1 = harp.import_product(S5P_OFFL_L2__NO2____20200301T122239_20200301T140410_12345_01_010302_20200304T172320.nc,
operations="latitude > 32 [degree_north]; latitude < 60 [degree_north]; tropospheric_NO2_column_number_density_validity > 75;bin_spatial(231,-55,0.5,721,-180,0.5)",
post_operations="bin(); squash(time, (latitude,longitude))")Cabe prestar especial atención a la operación bin_spatial(), que es la encargada de homogeneizar todas las variables del producto en una malla regular y facilitar de este modo su correcta visualización.
Aplicamos también un filtro a la variable tropospheric_NO2_column_number_density_validity, para importar solamente aquellos píxeles cuyo valor es superior a 75. Se recomienda aplicar este filtro para descartar la presencia de nubes, hielo y nieve.
Podemos importar tantos productos como necesitemos. Si los almacenamos primero en una lista, posteriormente los podremos juntar y exportar el resultado en formato NetCDF.
# almacenamos los productos importados en una lista # ambos productos cubren el territorio español products = [no2_escena1,no2_escena2] # juntamos los productos. product_bin = harp.execute_operations(products, "", "bin()") # exportamos el resultado de juntar ambos productos harp.export_product(product_bin, resultado_union_dia1.nc)
De este modo conseguimos crear un mosaico con los 2 productos Sentinel-5p que cubren el territorio objeto de estudio durante un día.
Este código debería repetirse para cada uno de los días a analizar.
Fusionar los mosaicos manteniendo la variable temporal
Una vez finalizado el paso anterior, dispondremos de 96 archivos, uno para cada uno de los días a analizar. Ejecutando el siguiente comando desde la terminal y con el entorno de Conda activo, los podremos juntar en uno solo manteniendo la variable temporal. De este modo, los podremos visualizar de forma animada:
harpmerge resultado_union*.nc mosaico_temporal.nc
Optimizar el tamaño de los productos
Un solo producto Sentinel-5p tiene un tamaño aproximado de 400Mb. Teniendo en cuenta que para cubrir España son necesarios 2 de estos productos, necesitaremos 800Mb de disco duro para almacenarlos. Si multiplicamos este valor por 96 (los días que queremos analizar), la cifra se nos dispara hasta los 77GB aproximadamente.
Sin embargo, al realizar la importación de los productos a HARP, filtrarlos y acotar su extensión, este tamaño disminuye significativamente.
El archivo anterior, utilizado para representar el NO2 de todo el territorio Español el día 20 de mayo de 2020 y que hemos procesado con HARP, tiene un tamaño de 123Mb. Una cifra sustancialmente inferior.
En caso de no disponer de mucho espacio de almacenamiento en el equipo de trabajo, una opción es programar el script que ejecuta el análisis para que, a medida que se vayan generando los mosaicos diarios, se eliminen los archivos originales. Así se irá liberando espacio en el sistema.
Programar una rutina con Python
El hecho de utilizar estas librerías facilita poder desarrollar scripts que automaticen flujos de trabajo.
En este enlace encontraréis un ejemplo de script en Python que realiza la descarga de 322 productos Sentinel-5p que cubren un extenso territorio (Europa y buena parte de África). Estos productos son los capturados entre el 23 de diciembre de 2019 y el 31 de mayo de 2020, y resultarán útiles para evaluar el descenso de contaminación por NO2 durante la pandemia del COVID-19.
Observaréis que el script incorpora funciones para automatizar la descarga de los datos, así como su importación y procesado con HARP. También una función para crear un mosaico con el valor promedio de 7 dias del NO2 concentrado en la zona objeto de estudio.
La construcción de estos mosaicos con el valor promedio de 7 días está programada para ser ejecutada a medida que los productos se van descargando. Y una vez creado el mosaico (en formato NetCDF) se borran los archivos originales ya utilizados. De este modo no es necesario disponer de un espacio de almacenamiento excesivamente elevado.
La ejecución de este código se demorará bastante tiempo. Y si queréis ampliar el número de días a analizar, tardará aún más. Es cuestión de tener paciencia y esperar a los resultados.
Visualizar los datos con VISAN
VISAN es una herramienta para visualizar y analizar datos atmosféricos.
Una vez instalada en el entorno de Conda, podremos arrancar la aplicación tecleando la instrucción visan en la terminal.
Con el programa arrancado, podemos importar a VISAN los ficheros NetCDF generados con HARP. La importación se hará desde el menú File > Harp Import… En file path indicaremos la ruta del archivo, y en variable name el nombre con el que queremos identificar el producto desde VISAN (por ejemplo, dia_1..nc). Finalmente, se tendrá que hacer clic en Import.
Una vez realizada la importación del archivo a VISAN, dispondremos de una consola para introducir las instrucciones en lenguaje Python.
Para visualizar la imagen, simplemente hay que utilizar:
wplot(*args, **kwargs)
A continuación se puede observar la animación de la evolución media semanal (diciembre 2019 – mayo 2020) de NO2 registrado por el sensor TROPOMI a bordo del Sentinel-5p:
Conclusiones
A través de las imágenes de satélite podemos analizar una amplia variedad de fenómenos. Incluso aquellos que a priori no se percibirían directamente desde el espacio, como una pandemia, pueden también ser analizados a partir de alguna de sus consecuencias. En este caso, la reducción de la actividad humana debido al confinamineto ha implicado una reducción de las emisiones de NO2 que tal y como hemos visto, puede observarse fácilmente.
Asimismo, la gran disponibilidad de imágenes que tenemos actualmente a nuestro alcance, junto a las herramientas adecuadas para procesarlas, nos ofrecen una amplia variedad de posibilidades que podemos explorar aplicando técnicas de teledetección. También en el ámbito de la epidemiología.
