Aplicaciones de Teledetección con el QGIS y el plugin Semi-Automatic Classification

En esta oportunidad y motivado por una conversación con alguien sobre éstos temas, voy a dedicarme a mostrarles lo que se puede hacer con el plugin Semi-Automatic Classification (SCP), creado por Luca Congedo, sobre todo en su última versión 4.3, el cual desde su propia página podemos encontrar una guía de usuario bastante completa.

Nuestro objetivo va ser además de conocer sobre las funcionalidades del plugin SCP, poder realizar la descarga directa de imágenes del Landsat y luego estimar la temperatura de la superficie de la tierra empleando la banda TIRS del Landsat 8.

Conociendo el entorno del Semi-Automatic Classification Plugin (SCP)

En resumen este complemento está compuesto por una barra de herramientas (The Toolbar), dos paneles móviles integrados en QGIS (The ROI Creation dock y Classification dock) y una ventana de la interfaz principal, juntos permiten realizar una clasificación supervisada de imágenes de manera fácil y rápida, con funciones de pre-procesamiento y post-procesamiento para los datos. Por otra parte, los gráficos de firmas espectrales y de dispersión se muestran en ventanas separadas.

1. Una barra de herramientas (The Toolbar)

El cual permite la selección de la imagen de entrada, e incluye varios botones para la apertura de las principales funciones de la ventana de la interfaz principal.

Figura 1: Barra del SCP

2. Panel para creación de Regiones de Interés (The ROI Creation dock)

Permite la definición de un archivo Shapefile de entrenamiento, y para la creación de los ROIs usando un algoritmo crecimiento de regiones (1) o por dibujo manual.

Figura 2: El panel para crear ROIs

3. Panel de Clasificación (Classification dock)

Está diseñado para gestionar las firmas espectrales, y clasificar las imágenes del satélite (o un conjunto de bandas). El proceso de clasificación requiere de las firmas espectrales que definen las clases de cobertura del suelo. Adicionalmente, es posible importar las firmas espectrales de otros proyectos o desde la pestaña USGS Spectral Library.

Figura 3: El panel de Clasificación

4.  La Ventana de Interfaz Principal

Se compone de varias pestañas agrupadas en secciones. Cada sección contiene varias funciones que son útiles principalmente para el proceso de clasificación. En primer lugar se cuenta herramientas (Tools) específicas para la teledetección, en la última versión se cuenta con la facilidad para descargar directamente imágenes Landsat,  del mismo modo existen opciones para realizar trabajo de pre- procesamiento de imágenes, de procesamiento posterior, para realizar operaciones con las bandas (Band calc), para lograr definir el grupo o conjunto de bandas como valores de entrada y el centro de las longitudes de onda de las bandas esto último requerido para calcular mejor las firmas espectrales (Band set).

Figura 4: Interfaz principal del plugin SCP

A. Descargando imágenes del Landsat 8 con SCP.

Ahora vamos a realizar un ejercicio, para ello necesitamos descargar una imagen proveniente del satélite Landsat 8, además de poder hacerlo por EarthExplorer, lo vamos a realizar desde el SCP para demostrar que es muy sencillo y práctico.

Paso 1: Instalar el complemento

Para ello  verificamos que tenemos instalado el complemento Semi-Automatic Classificaction Plugin (SCP), si no lo tenemos, tenemos que instalarlo desde ComplementosàAdministrar e instalar complementos, realizamos la búsqueda por nombre y lo instalamos.

Figura 5: Zona de instalación del plugin SCP

  

Paso 2: Realizar la búsqueda de las imágenes del Landsat

Una vez que Ahora ya estamos listos para iniciar con el ejercicio, para ello nos vamos a la barra de menú del QGIS y en SCP buscamos Tools y luego Download Landsat.

Figura 6: Herramienta del SCP para descargar imágenes del Landsat.

Una vez realizado esto, nos aparecerá un entorno en donde vamos a poder realizar la búsqueda de nuestra imagen, teniendo en cuenta que para nuestro caso solamente buscaremos una imagen del Landsat 8 y debido a que vamos a realizar un ejercicio empleando la banda termal (B10), la fecha de adquisición no debe ser actual debido a los problemas con dichas bandas, según la USGS.

Figura 7: Ajustes antes de realizar la búsqueda de la imagen.

De acuerdo a la Figura 7, debemos primero marcar la casilla de "only Landsat 8", luego generar un directorio para la base de datos, que permite guardar las listas de las búsquedas realizadas, luego determinar los límites de la zona que nos interesa, para ello haremos clic tanto para determinar las coordenadas del límite superior izquierdo (UL) como las del límite inferior derecho (UR). Para nuestro caso estoy seleccionando una zona del territorio peruano y como ayuda superpuse una capa vectorial de los límites de departamento.

Figura 8: Ajustando nuestra zona de interés 

Ahora nos queda elegir el rango de fechas y también podemos indicarle un valor al valor máximo de cobertura de nubes (Max cloud cover), una vez hecho esto hacemos clic en "Find images". Como producto de ello nos van a salir una lista de imágenes, para lo cual escogemos la que nos interesa, inclusive podemos hacer una vista preliminar. Luego se debe remover de la lista las imágenes que no vamos a emplear.

Figura 9: Eligiendo la imagen y mostrando una pre-visualización

Paso 3: Realizar la descarga de la imagen del Landsat seleccionada

Como paso siguiente debemos realizar la descarga de la imagen (LC80080672014353LGN00), para nuestro caso solamente necesitaremos las bandas 2, 3, 4, 5, 6, 7 y la banda termal 10. 

Figura 10: Eligiendo las bandas que vamos a descargar

Existe también la opción de realizar el pre-procesamiento de las mismas activando la casilla respectiva, como resultado de ello vamos a tener la conversión automática de los DN (es decir, los números digitales) de las bandas 2 al 7 del Landsat para la medida física de la reflectancia de la parte superior de la atmósfera (TOA) y para las bandas del sensor TIRS (Thermal Infrared Sensor) serán convertidas desde la radiancia espectral a la temperatura de brillo usando las constantes térmicas previstas en el archivo de metadatos (LC80080672014193LGN00_MTL). Además, el SCP implementa un simple corrección atmosférica basada en imágenes utilizando el método DOS1 (Resta del Objeto Oscuro-1). Para mayor detalle de éstos temas revisar aquí. 

Una vez realizado este paso tendremos el siguiente resultado.

Figura 11: Resultados de la descarga

Como se activó la casilla para abrir las imágenes en el QGIS se obtendrá la siguiente figura.

Figura 12: Visualización de la imagen descargada y luego del pre-procesamiento

Finalmente los resultados obtenidos del pre-procesamiento se van a guardar en una carpeta distinta.

Figura 13: Carpeta generada luego del pre-procesamiento

Para ver una demostración recomiendo ver el vídeo realizado por el mismo creador del plugin: 

Download Landsat 8 Images with Semi-Automatic Classification Plugin v.4 for QGIS

B. Estimación de la temperatura de la superficie de la tierra empleando la banda TIRS del Landsat.

Para lograr realizar este ejercicio tendremos que realizar primero una clasificación supervisada de las coberturas presentes en la imagen, esto es para definir la emisividad de la superficie, requerido para el cálculo de la temperatura de la superficie del terreno.

Paso 1: Preparar nuestros datos de entrada

Debido a que no se requiere emplear toda la imagen vamos a recortar nuestra zona de interés, para ello tenemos varias opciones, para nuestro caso vemos que dentro de la pestaña del SCP dedicado al pre-procesamiento encontramos la opción que dice "Clip multiple rasters". Pero antes de realizar el corte podemos fijarnos los metadatos de las imágenes y comprobamos que tienen como SRC (Sistema de referencia de coordenadas) WGS 84 / UTM zone 18N (EPSG: 32618), por lo tanto debemos realizar una reproyección de coordenadas para pasarlo a un SRC: WGS 84 / UTM zone 18S (EPSG: 32718).

Figura 14: Realizando una reproyección de nuestra imagen

Ahora vamos a realizar el corte de la zona de interés, en esta oportunidad se generó un polígono que abarca la zona en donde se ubica la ciudad de "Carhuaz". Solo debemos seleccionar las bandas que van a ser recortadas y el polígono que se preparó, también lo pueden hacer generando sus propias coordenadas.

Figura 15: Realizando un recorte de la imagen.

Para quienes desean seguir la metodología, comparto los datos obtenidos después de realizar el recorte junto con la metadata de la imagen y el polígono de recorte.

Datos de entrada. 

La imagen resultante del recorte se debe apreciar como la siguiente figura.

Figura 16: Resultado del recorte de la imagen

Paso 2: Clasificación de las coberturas del área de estudio.

Con nuestros datos de entrada tenemos que generar nuestro conjunto de bandas, para lo cual nos vamos a la pestaña "Band set", antes que todo vamos a refrescar la lista (Refresh list), hay que tener en cuenta que son imágenes recortadas y eso se debe mostrar en la lista, luego de ello seleccionamos todos y hacemos clic en "Add rasters to set", pero no hay olvidarse de remover la banda 10. Luego ordenamos las bandas de menor a mayor, a continuación seleccionamos el sensor Landsat 8 OLI del cuadro combinado “Quick wavelength settings", con el fin de ajustar automáticamente la longitud de onda central de las bandas. 

Figura 17: Definición de las bandas de entrada para la clasificación supervisada

Para apoyarnos en el desarrollo de la clasificación vamos a generar una imagen virtual empleando tres bandas para crear una composición, para lograr esto nos vamos a la barra de menú dentro de Ráster--> Miscelánea--> Construir ráster virtual (catálogo), luego seleccionamos las bandas 3, 4 y 5 e indicamos el nombre de salida, en este caso es "rgb.vrt".

Figura 18: Generando una imagen virtual componiendo 3 bandas

Con el resultado obtenido se ajustaran las bandas para resaltar sobre todo la vegetación, realizando un RGB 5-4-3, también si deseamos podemos superponer una imagen del google earth.

Figura 19: Imagen virtual compuesta RGB: 543

Con esta imagen ya podemos iniciar la creación del archivo Shapefile de las zonas de entrenamiento. Dentro del panel de creación ROI hacer clic en el botón “New Shp”, y seleccionar dónde desea guardar el archivo Shapefile (ROI_CARH.shp).

Figura 20: Generando el polígono para la creación de ROIs

Dentro del panel de Clasificación “SCP: Classification” hacer clic en el botón Guardar “Save" con el fin de crear un archivo de lista de firmas espectrales (sig1.xml).

Figura 21: Generando el archivo para guardar las firmas espectrales

Una vez realizado estos pasos podemos iniciar la recolección de las áreas de interés. Los ROI (regiones de interés) son polígonos que se pueden crear de forma automática tal como se indicó anteriormente mediante un algoritmo de crecimiento de regiones (es decir, la imagen es segmentada en torno a una semilla de píxeles), haciendo clic en un píxel de la imagen. De esta manera, los ROI creados incluyen píxeles espectralmente homogéneos.  

El SCP permite la definición de un identificador de Macroclases (es decir, MC ID) y un identificador de Clase (es decir, C ID) para cada ROI o firma espectral, que son los códigos de identificación de las clases de cobertura del terreno. Las Macroclases permiten la clasificación de los materiales que tienen diferentes firmas espectrales (por lo tanto, son procesados de forma individual), pero pertenecen a la misma clase de cobertura del terreno (por lo tanto el mismo MC ID se asigna a estos píxeles). Por ejemplo podríamos clasificar pastos asignándole C ID = 1 y MC ID = 1 y árboles asignándole C ID = 2 y MC ID = 1, todo esto como una Macroclase de vegetación (es decir, MC ID = 1). Cada ROI (o firma espectral) deben tener un C ID único, mientras que el MC ID puede ser compartido con otros ROIs. En el panel de Clasificación es posible elegir entre la clasificación MC ID y C ID. 

Con el fin de crear un ROI, en la panel de creación ROI hacer clic en el botón “+” junto a “ROI creation”  y, a continuación, hacer clic en cualquier píxel de la imagen; acercar el mapa lo suficiente para ayudarse en la identificación de las coberturas.

En “ROI signature definition” escriba una breve descripción del ROI dentro del campo de C ID (Información de Clase) y también dentro de la Información de Macroclases (MC ID),  es decir asignar un ID Macroclass y un  ID Clase (es posible cambiar el nombre y cambiar estos códigos después desde “ROI list”, el cual  también cambia los valores en el archivo Shapefile de entrenamiento).

Figura 22: Generando  las Clases de los ROIs

Con el fin de salvar el ROI del archivo shapefile de entrenamiento, hacer clic en el botón “Save ROI”; si la casilla de verificación “Add sig. list” está marcada, entonces las firmas espectrales son calculadas (la media de los valores de los píxeles del ROI para todas las bandas) y se añade a la tabla de la lista de firmas; también, es posible añadir firmas más tarde, resaltando los ROIs de la lista y haciendo clic en “Add to Signature” (si dos o más regiones de interés (ROI)  tienen el mismo MC ID y C ID, una firma espectral única se calcula considerando todos los píxeles que están alrededor de esas regiones de interés).

Figura 23: Grabando la lista de los ROIs generados por clases

Dentro del panel “SCP: Classification” podemos definir el color de las clases que se utilizarán en la clasificación, para ello hacer doble clic sobre la columna de color en la lista de firmas “Signature list” (la lista de firmas se guarda automáticamente cuando se guarda el proyecto de QGIS, o al hacer clic en el botón “Save” ).

Figura 24: La lista de las firmas espectrales generados ajustando los colores

Debemos repetir los pasos anteriores para cada clase de cobertura del terreno y asignar a cada ROI una nueva ID de clase de manera incremental, y de las siguientes ID de Macroclases:

• Suelo (por ejemplo, desnudo, rocoso): MC ID = 1
• Vegetación (por ejemplo, árboles, cultivos): MC ID = 2
• Agua (por ejemplo, río): MC ID = 3
• Construcción (por ejemplo, viviendas, asfalto): MC ID = 4

El SCP  permite realizar vistas previas de clasificación, a fin de evaluar muy rápidamente los resultados de la clasificación. Estas vistas previas de Clasificaciones son útiles durante la recolección de ROIs, y para la selección de las firmas espectrales más precisas. Si los resultados de las vistas previas se consideran buenas (es decir, las clases se identifican correctamente), la clasificación de la totalidad de la imagen se puede realizar. Por otra parte, es posible eliminar una o más firmas espectrales, o añadir nuevas firmas espectrales creando otros ROIs tal como se describe en los pasos anteriores.

En el panel de Clasificación, dentro “Classification Preview” puedes colocar por ejemplo Set size = 200 (es decir, al lado de la vista previa de clasificación en unidades de píxeles), y luego seleccionar uno de los algoritmos disponibles, para nuestro caso consideramos que “Minimum Distance” obtuvo los mejores resultados; hacer clic en el botón +, a continuación, hacer clic en una zona de la imagen; después de unos segundos, se mostrará la vista previa de clasificación. 

Figura 25: Pre-visualización del resultado de clasificación

Para llevar a cabo la clasificación final, hacer clic en el botón “Perform Classification” y seleccionar dónde desea guardar la salida. Para nuestro caso se optó que el resultado se mostrara a nivel de Macroclases, por lo tanto la pestaña que indica "Macroclass ID", estaba activada. Finalmente se obtuvo lo siguiente:

Figura 26: Resultado final de la Clasificación de coberturas

Paso 3: Reclasificación de los valores de clasificación de coberturas por sus valores de emisividad.

Los valores de emisividad (e) se presentan en la siguiente tabla (éstos valores son indicativos ya que se deben obtener de estudios de campo):

Superficie de la tierra	Emisividad (e)
Suelo	0.93
Vegetación	0.98
Construcción (Brick: masonry)	0.94
Agua	0.98

Fuente: http://www.infrared-thermography.com/material-1.htm

Con estos datos se requiere realizar la reclasificación, para lo cual vamos a trabajar con la "Caja de herramientas de procesado", empleando uno de los algoritmos de SAGA "Reclassify grid values".

Figura 27: Ubicación del algoritmo de SAGA para realizar una Reclasificación

Una vez que definimos el ráster de entrada que sería el resultado de nuestra clasificación de coberturas, dentro de las opciones sobre los métodos que nos presenta vamos a elegir "Simple table", para lo cual vamos a requerir rellenar una tabla (5x3) con los valores de la clasificación para ser cambiados con los valores de emisividad.

Figura 28: Generando una tabla 5x3 para realizar la reclasificación

Con todo ello obtendremos el siguiente resultado en el QGIS:

Figura 29: Resultado final de la reclasificación.

Paso 4: Convertir la temperatura de brillo del satélite (brightness temperature) a la temperatura de la superficie de la tierra.

La emisividad corregida por las temperaturas de la superficie del suelo  (The emissivity-corrected land surface temperatures - LST), puede ser calculado con la siguiente formula desarrollado por: Artis and Carnahan 1982 (2), Weng et al. 2004 :

Donde:

T_B = The radiant surface temperature (K)

λ = Wavelength of emitted radiance (ver tabla)

ρ = h x c/ σ (1.4386 10^-2 mK)

σ = Boltzmann’s constant (1.38 x 10^-23 J K^-1)

h = Planck’s constant (6.626 x 10^-34 Js)

c = the velocity of light (2.998 x 108 m s-1)

ε = Emisividad espectral

Los valores de λ para las bandas termales del satélite Landsat están listados en la siguiente tabla:

Satellite	Band	Center wavelength (µm)
Landsat 4, 5, and 7	6	11.45
Landsat 8	10	10.8
Landsat 8	11	12

Para profundizar en el tema recomiendo revisar la siguiente investigación:

Land surface temprature estimation from thermal bandof landsat sensor,case study: alashtar city

Para nuestro caso vamos a trabajar con la Calculadora ráster del QGIS, en donde introduciremos la siguiente formula:

b / (1 + (10.8 * b / 14380) * ln(a))

Donde a  es el ráster obtenido de la reclasificación y b es el ráster de la temperatura de brillo.

Figura 30: Empleo de la calculadora ráster del QGIS

Finalmente obtendremos el siguiente mapa como resultado de este cálculo, el cual para una mejor visualización se hizo algunos ajustes al estilo.

Figura 31: Resultado final para el cálculo de la temperatura de la superficie de la tierra.

Bueno espero que puedan desarrollar la metodología con los datos presentados o con sus propios datos y de ser posible comenten al respecto. A continuación voy a indicar las referencias que se han revisado y son la fuente principal de la presente entrada y en donde seguro encontraran mayores detalles.

Referencias:

• SCP 4.3 allows for downloading from the entire archive ofLandsat 4, 5, 7 and 8
• Estimation of Land Surface Temperature with Landsat Thermal InfraredBand: a Tutorial Using the Semi-Automatic Classification Plugin for QGIS
• Surface Temperature from Landsat Data: a New Lab Using the Semi-AutomaticClassification Plugin
• Lab 4:remote sensing land surface temperature

---

(1) El crecimiento de regiones es un proceso de segmentación contextual tendiente a separar objetos específicos de la imagen.

(2) Artis, D. A., & Carnahan, W. H. (1982). Survey of emissivity variability in thermography of urban areas. Remote Sensing of Environment, 12, 313– 329