jueves, 28 de junio de 2018

Herramienta para visualizar datos ráster en el tiempo - EO Time Series Viewer



En esta oportunidad se muestra una herramienta disponible desde el QGIS 3 como complemento que permite visualizar imágenes de distintos sensores remotos con diferentes fechas de adquisición, nos referimos al EO Time Series Viewer, desarrollado por Benjamin Jakimow. En general el objetivo de esta entrada será demostrar el uso de la herramienta que nos ayuda en la visualización de un grupo de imágenes provenientes de los sensores del Sentinel2 y el Landsat 8, los cuales han sido adquiridos durante los últimos meses.

Datos de Ejemplo


Para el desarrollo de esta demostración se descargó de las fuentes disponibles las siguientes imágenes del Sentinel 2:
  1. S2A_MSIL1C_20180201T153611_N0206_R068_T17MPN_20180201T190246
  2. S2A_MSIL1C_20180502T153621_N0206_R068_T17MPN_20180502T221752
  3. S2A_MSIL1C_20180522T153621_N0206_R068_T17MPN_20180522T203523
  4. S2B_MSIL1C_20180127T153609_N0206_R068_T17MPN_20180127T184619
  5. S2B_MSIL1C_20171218T153609_N0206_R068_T17MPN_20171218T184243
Para el caso del Landsat 8 se uso : LC08_L1TP_010065_20180227_20180308_01_T1

Con todas las 6 imágenes se realizó el preprocesamiento necesario, luego fueron recortados en una cuadrícula de 20x20 km de una zona de interés, ubicada en el Departamento de Lambayeque en Perú. Luego de recortar todas las bandas, las mismas fueron apiladas (stack) para contar con las imágenes compuestas. Todo esto luego fue ordenado en grupos dentro de nuestro panel de capas, los cuales se pueden apreciar en la siguiente figura.


Primeras Experiencias


Antes de encontrar a EO Time Series Viewer, cuando intentaba evaluar diversas imágenes satelitales una herramienta que me ayuda mucho fue Temporal/Spectral Profile Tool, a través del cual se puede realizar el trazado interactivo de información temporal o espectral almacenada en imágenes multibandas. Cuando lo activamos podemos visualizar el trazado de los datos que se extraen de las imágenes como multibandas, desde donde podemos ir incorporando más de una imagen.


Incluso se cuenta con la ventaja de que podías extraer lo valores para ser copiados en el portapapeles y tener nuestro propio cuadro con los datos espectrales extraídos.



Pero cuando se incorpora la imagen Landsat 8, podemos apreciar que el gráfico se distorsiona un poco para lograr ajustar al número de bandas y en este caso como este último tiene solo 6 en lugar de los 10 del Sentinel 2.



Conociendo a EO Time Series Viewer


Luego de instalar el complemento junto con sus dependencias (paquetes de python: pyqtgraph y pyopengl), para lo cual se recomienda revisar aquí los pasos a seguir. Cuando activamos el complemento nos aparece la siguiente ventana con la interface de usuario.


Al igual que el QGIS, los paneles que lo componen pueden ser ajustados para que se ubiquen en distintos lugares, del mismo modo pueden ser activados y desactivados.

Sin pretender ser una guía detallada del complemento, una información más detallada lo pueden obtener desde aquí. Lo que se busca es resaltar las principales opciones y ventajas que se tiene al emplear este complemento.

Primeros Pasos


Lo primero que debemos hacer es incorporar nuestras imágenes, para ello empleamos la opción denominada "Add images to the time series". Luego dentro de "Sensors/Products", vamos a tener que editar el cuadro inicial indicando el satélite desde donde provienen las imágenes.


Lo interesante de esta herramienta es para visualizar las imágenes, a las mismas lo ordena por fecha de adquisición, sin importar si se trata de una imagen Landsat o de Sentinel. Si tenemos activado la pestaña "Time Series", se puede apreciar en la parte inferior una tabla que nos brinda información básica como el sensor, el número de bandas (nb), los píxeles en dirección x (ns) y en dirección y (nl), el sistemas de referencia (CRS) y al final la denominación de nuestra imagen.



Visualización de Mapas


Quizás el panel denominado "Maps" es uno de los más importantes por la variedad de opciones disponibles. Lo primero que vamos a resaltar es el hecho que podemos incorporar distintas vistas, todos ellos pueden diferenciarse por la combinación de bandas a emplearse. En este caso generamos tres vistas que nos permiten ver nuestra imagen con el zoom más adecuado con tres tipos de combinaciones. Es importante resaltar que para cada tipo de sensor se activa la opción para seleccionar las bandas disponibles.
En la imagen de abajo se muestra el detalle de una vista denominada "Agricultura", en donde para Landsat 8 (6,5,2), mientras que para el Sentinel 2 (11,8A, 2). Para otras combinaciones se puede revisar desde aquí. Tener en cuenta que la combinación de las bandas que se muestra está en función al orden de apilado realizado para cada sensor.



Otra ventaja que tenemos es la posibilidad de adjuntar una capa vectorial, sobre todo si tenemos identificado una zona de interés que nos permitirá realizar una mejor evaluación. En este caso a modo de ejemplo se adjunta una distribución de parcelas. Lo interesante es que podemos ajustar un estilo tal como lo haríamos en el QGIS.


Otro detalle que podemos resaltar es que nos permite contar con un punto de mira detallado y que puede ser personalizado con resaltado de píxeles y barra de escala. Se puede apreciar que se diferencia en función a la resolución de la imagen.



Existe un panel dedicado a identificar los valores en una posición específica, nos referimos a "Cursor Location Values", con el mismo concepto de QGIS cuando usamos una herramienta para identificador de valores ráster (por ejemplo "Value tool"). Nos brinda la información del píxel en las coordenadas que se indica.




Librería Espectral


Finalmente vamos a comentar sobre el panel denominado "Spectral Library", desde donde podemos generar la firma espectral en función a las bandas disponibles. Por defecto se encuentra activa la opción denominada "Add new profiles automatically to opened Spectral Library", con el cual podemos generar los perfiles en los puntos seleccionados sobre las imágenes, cada vez que hacemos un clic sobre la imagen se irán creando los valores, los mismo que luego pueden ser exportados en formatos de texto.


También podemos emplear la opción denominada "Select new profile from map", para visualizar una gráfica que irá variando según el punto seleccionado en nuestras imágenes. Por ejemplo si nos ubicamos dentro de una zona en donde existe presencia de vegetación, la curva tendrá la siguiente forma.





Ahora si ubicamos en una zona sin vegetación tendremos la siguiente gráfica.



Perfil Temporal


Para el caso del panel "Temporal Profiles", es interesante tener la opción de generar un gráfica que nos permite tener los valores (DN) en el tiempo (Date y DOY) para una determinada ubicación, pero sobre todo poder combinar nuestras bandas empleando para este caso el índice NDVI (para el sensor Sentinel). También se cuenta con la opción de establecer un estilo a la gráfica de perfil que se muestra.





Bueno, espero que puedan probar este complemento, hasta el momento puedo decir que me está ayudando mucho cuando requiero analizar varias imágenes al mismo tiempo, sobre todo sería muy útil para el monitoreo de coberturas en el tiempo.




martes, 2 de enero de 2018

Aplicación de máscaras para nubes en imágenes Landsat 8 con QGIS









En esta oportunidad vamos a presentar a Cloud Masking,  un plugin de QGIS con el que podemos enmascarar nubes, sombras de nubes, cirros (nubes grises), aerosoles, hielo/nieve y agua para productos Landsat (4, 5, 7 y 8), utilizando diferentes procesos y filtros como Fmask, Blue Band, Cloud QA , Aerosol y Pixel QA.

Fuente de datos



Para este procedimiento vamos a considerar una escena del Landsat 8, las cuales pueden ser descargadas desde distintas fuentes, en este caso se realizó desde: https://earthexplorer.usgs.gov/. La zona de interés corresponde al Path/Row 8/64.




Descarga de las bandas de Landsat con ESPA



Con el objetivo de conseguir mejores resultados para el procesamiento de enmascaramiento y aplicación de los filtros, se requiere otros productos o bandas del Landsat, los cuales deberás ser solicitados a través de la interface del Earth Resources Observation and Science (Eros) Center Science Processig Architecture (ESPA),  para lo cual es necesario tener un registro.




Con nuestros datos de registro, accederemos a una página que nos solicita que como primer paso le indiquemos la lista de escenas, en este caso del Landsat 8, para ello debemos generar un archivo de texto (*.txt), con los códigos que identifican las escenas; para nuestro caso sería algo como la siguiente figura.


Ahora cuando nos soliciten la lista de las escenas, simplemente tenemos que indicar la ubicación de nuestro archivo de texto. Para las otras opciones, considerar las que se indican en la figura de abajo, en donde para este caso también se solicitó en la sección de configuración de salida, que nuestras escenas se guarden como GeoTiff y tengan una proyección UTM Zona 18S.



Luego de realizar el pedido, nos llegará un correo anunciando que el enlace para la descarga está listo y nosotros ya podemos descargar todo nuestro pedido.




Plugin Cloud Masking



Una vez que ya contamos con nuestras bandas requeridas, nuestro siguiente paso es instalar y luego activar el plugin Cloud Masking. Tal como se aprecia en la figura de abajo, se genera un panel del plugin.




El plugin presente esta dividido y ordenado en tres secciones; (1) Open and Load, desde donde podremos abrir el archivo de los metadatos (MTL) y las bandas apiladas, (2) Filters and Mask, que nos permite habilitar y configurar los filtros para aplicar las máscara de nubes y (3) Apply and Save, para guardar la máscara creada y aplicarla sobre las bandas apiladas.


Abriendo los Metadatos 


Lo primero que tenemos que hacer con el plugin es ubicar el archivo de los metadatos de nuestra escena (MTL), para abrirlo e identificar nuestras bandas.



Abriendo un Apilado de bandas (Stack)



No es un paso obligatorio, pero se recomienda hacer una combinación de bandas de tipo RGB, con el objetivo de poder visualizar y revisar la escena. En el lado derecho podemos apreciar que existen botones con las combinaciones más comunes. Para este caso vamos a elegir "False Color", luego hacemos clic en "Load stack", con ello podremos apreciar en nuestra vista de mapa, la combinación de bandas establecida.




Seleccionando los filtros a aplicarse



Ahora en la siguiente sección debemos seleccionar nuestros filtros. El plugin tiene cuatro filtros diferentes para aplicar, Fmask y Blue Band están disponibles para cualquier escena Landsat, pero Aerosol y Pixel QA (para el caso del Landsat 8), están disponibles solo para las escenas que son solicitadas a través de ESPA, tal como lo vimos anteriormente. Se puede activar más de un filtro a la vez, el complemento acumula los filtros en el mismo orden (de abajo hacia arriba).



Filtro FMask


El primer filtro que activaremos será el FMask, dejando las opciones por defecto. El proceso de Fmask utiliza una implementación de un algoritmo de python llamado fmask, descrito  en http://pythonfmask.org/en/latest/, actuando como una biblioteca interna en el plugin.




Filtro de la Banda Azul


Este filtro utiliza la banda azul del Landsat (Banda 2 para Landsat 8) para enmascarar todos los píxeles con valores inferiores al umbral establecido, este último también depende de la versión empleada, por lo tanto, para Landsat 8 el umbral es 0-65534 (16 bits).




Filtro Aerosol



Este tipo de filtro solo se puede ejecutar en escenas del Landsat 8, para lo cual se emplea una banda (Aerosol) de 8 bits. Cuando revisamos los archivos solicitados a través de ESPA, veremos que existe una banda que tiene la terminación  *_sr_aerosol.tif. Estos 8 bits son:



En el complemento se implementa este filtro bit a bit (solo los bits útiles) y se puede habilitar uno o más bits.



Para obtener más información, es importante consultar esta guía del producto.


Filtro Pixel QA


Para emplear este filtro, debemos contar con una banda de 16 bits, la misma que presenta la terminación *_pixel_qa.tif. Este archivo solo estará disponibles cuando seleccionamos en nuestro pedido a través de ESPA, la opción de "Surface Reflectance (SR)", lo que representa una colección especial de Landsat. Estos 16 bits para Landsat 8 son:



Las opciones del Pixel QA a emplear son:





Generando la máscara de nubes


A partir de aquí ya podemos generar las máscaras en función a los filtros previamente habilitados. Existen tres opciones para su generación: (1) La imagen completa (Whole Image), la cual viene por defecto, esto ocurre si las otras opciones están desactivadas; (2) Solo en el área seleccionada (In only select area), cuando se selecciona esta opción, podemos procesar la máscara en un área específica de la imagen, que podrá ser seleccionada en la vista del mapa (solo un área rectangular). Si la opción de Mantener el tamaño de la imagen original (Keep the original image size) está habilitada, la máscara de resultado tendrá la misma extensión, de lo contrario, la máscara de resultado se recortará en el área seleccionada; (3) La opción In only shape area, con esta opción, se puede procesar la máscara solo dentro de un Shape particular, el mismo que debe ser creado antes de aplicarse. Si la opción Eliminar datos fuera de la forma está habilitada, la máscara resultante tendrá la misma extensión pero los datos fuera del shape se enmascararán, de lo contrario la máscara resultante tendrá la misma y los datos fuera de la forma serán válidos (solo enmascarados dentro del shape).

En nuestro caso vamos a emplear la opción 2 y luego haremos clic en "Generate mask".




Si vemos nuestro panel de capas, veremos que se generó una capa temporal con el resultado de la máscara generada aplicando los filtros habilitados.




Aplicando la máscara


El plugin nos ayuda a crear, aplicar y guardar las  bandas apiladas con el Landsat cargado, teniendo tres opciones para aplicar la máscara sobre las bandas apiladas. (1) Raw Bands, que son las bandas originales de los productos landsat, uno mismo puede seleccionar el orden de apilado, por defecto son las bandas de reflectancia; (2) SR Bands, estos son los productos de Landsat obtenidos por ESPA, Surface Reflectancia (SR), esta opción solo está habilitada si está procesando este tipo de producto, y también se puede seleccionar el orden de la bandas para ser el apilado, por defecto son las bandas de reflectancia; (3) Particular file, esta última opción se aplica a una apilado o archivo de imagen en particular.

Para nuestro caso emplearemos la opción 2 y luego le indicamos la ubicación en donde se guardará nuestra imagen apilada con la máscara aplicada.





Finalmente cuando se ejecuta la máscara, podemos hacer una combinación de bandas y para verificar nuestro resultado podemos emplear el plugin "Value Tool", y si nos ubicamos en una zona donde existe nubes y se aplicó una máscara, comprobaremos que los valores son "cero".




Se recomienda grabar la capa temporal inicial con los filtros aplicados, para ello es importante salvar el estilo de la capa para que lo podamos agregar en otro proyecto.

Bueno, por ahora es todo lo que quería mostrarles, les invito a probar el plugin.


* Versión en Inglés elaborado por geogeek.xyz.

Referencias Consultadas: