Datos satelitales en GRASS GIS

Los datos satelitales en general vienen en formato raster, por lo tanto aplican las mismas reglas que vimos anteriormente. Los comandos i.* se orientan explícitamente al procesamiento de datos satelitales aunque algunos puedan usarse para otros datos raster.

Para ejemplificar el flujo de trabajo para procesamiento de datos satelitales en GRASS GIS, en esta primer sesión vamos a trabajar con datos del satélite Landsat 8.

Antes de empezar y para ganar tiempo, conectamos nuestro drive e instalamos GRASS en Google Colab.

# import drive from google colab
from google.colab import drive
# mount drive
drive.mount("/content/drive")

%%bash
DEBIAN_FRONTEND=noninteractive 
sudo add-apt-repository ppa:ubuntugis/ubuntugis-unstable 
apt update 
apt install grass subversion grass-dev
apt remove libproj22

Chequeamos el path de instalación de GRASS

!grass --config path

e instalamos otras librerías de python que nos serán útiles.

!pip install pygdal
!pip install landsatxplore

Datos para esta sesión

Escenas Landsat 8 (OLI)

• Fechas: 14/01/2020 y 02/03/2020
• Path/Row: 229/082 (Córdoba)
• CRS: UTM zona 20 N (EPSG:32620)

Historia de la mision Landsat

Lanzamientos de satélites Lansat desde 1972

El sistema de escáner multiespectral (MSS) a bordo del Landsats 1-5 disponía de cuatro bandas. El Thematic Mapper (TM) a bordo de Landsats 4 y 5 tenía siete bandas. El Enhanced Thematic Mapper Plus (ETM+) del Landsat 7 tiene 8 bandas y los Landsats 8 y 9 tienen 11 bandas. Fuente: https://landsat.gsfc.nasa.gov/satellites/landsat-9/landsat-9-bands/.

Comparación entre las bandas de todos los satélites Landsat

Nota

Más detalles sobre las misiones Landsat pueden encontrarse en: https://www.usgs.gov/landsat-missions

Manos a la obra

Iniciamos GRASS GIS

Iniciamos GRASS GIS en el proyecto posgar2007_4_cba y mapset PERMANENT

import os

# data directory
homedir = "/content/drive/MyDrive/curso_grass_2023"

# GRASS GIS database variables
grassdata = os.path.join(homedir, "grassdata")
project = "posgar2007_4_cba"
mapset = "PERMANENT"

# import standard Python packages we need
import sys
import subprocess

# ask GRASS GIS where its Python packages are to be able to run it from the notebook
sys.path.append(
    subprocess.check_output(["grass", "--config", "python_path"], text=True).strip()
)

Ahora sí, estamos listos para importar los paquetes de GRASS e iniciar una sesión:

# import the GRASS GIS packages we need
import grass.script as gs
import grass.jupyter as gj

# Start the GRASS GIS Session
session = gj.init(grassdata, project, mapset)

Corroboramos la proyección

# check the CRS
print(gs.read_command("g.proj", flags="p"))

Crear un nuevo mapset

Para diferenciar sesiones dentro del curso, vamos a crear un nuevo mapset llamado landsat8:

# Create a new mapset
gs.run_command("g.mapset", 
               mapset="landsat8",
               flags="c")

# Check we are in the mapset just created
gs.run_command("g.mapset",
               flags="p")

Listamos los mapsets accesibles para corroborar que tenemos acceso a PERMANENT

# list all the mapsets in the search path
gs.mapsets(search_path=True)

Nota

Desde cualquier mapset en el que estemos trabajando, siempre vamos a tener acceso a PERMANENT.

Listemos los mapas vectoriales disponibles

# list vector maps in all mapsets in the search path
gs.list_grouped(type="vector")

Región de interés

Como primer paso para trabajar con datos raster,necesitamos establecer la región computacional. Para ello, vamos a extraer el radio urbano de Córdoba del vector de radios urbanos y luego lo vamos a usar para establecer los límites de nuestra region computacional.

# extract Cordoba urban area from `radios_urbanos`
gs.run_command("v.extract", 
               input="radios_urbanos", 
               where="nombre == 'CORDOBA'", 
               output="radio_urbano_cba")

# set the computational region to the extent of Cordoba urban area
gs.run_command("g.region", 
               flags="p", 
               vector="radio_urbano_cba")

Descargar e importar los datos Landsat 8

Entre las extensiones disponibles, hay una extensión i.landsat, que nos permite buscar, filtrar, descargar, importar y aplicar bandas de calidad a imágenes Landsat. Para la busqueda y descarga de datos, esta extensión depende de una librería de python denominada landsatxplore.

Desafortunadamente, landsatxplore fue abandonada por su autor, y es difícil de mantener. Estamos considerando utilizar otras opciones, pero es un WIP. Por este motivo, vamos a instalar la extensión igualmente porque nos sirve para la importación y enmascarado de pixeles según la calidad.

Nota

Voluntari@s para actualizar y testear los pull request de landsatxplore o actualizar i.landsat.download para que use otro paquete distinto, son mas que bienvenid@s!!

# install i.landsat toolset
gs.run_command("g.extension",
               extension="i.landsat")

La búsqueda de escenas se basa en la región computacional definida y funcionaría como se detalla a continuación.

# search for Landsat 8 scenes
# l8_credentials = os.path.join(homedir, "gisdata", "landsat")
# gs.run_command("i.landsat.download", 
#               settings=l8_credentials,
#               dataset="landsat_8_c1",
#               clouds="35",
#               start="2019-10-27",
#               end="2020-03-15",
#               flags="l")

Por defecto, usando el comando anterior sin la opción l y proveyendo una carpeta de descarga, se descargarían todas las imágenes encontradas anteriormente. Sin embargo, también es posible descargar imágenes seleccionadas via su id.

# download selected scenes
# gs.run_command("i.landsat.download", 
#               settings=l8_credentials,
#               id="LC82290822020062LGN00,LC82290822020014LGN00",
#               output=os.path.join(homedir, "gisdata", "landsat_data"))

Como el módulo de descarga no está andando y porque la descarga puede tomar mucho tiempo con conexiones lentas, vamos a usar imágenes previamente descargadas. Como primer paso vamos a imprimir todas las bandas presentes dentro de la carpeta landsat_data.

# print all landsat bands within landsat_data folder
# (1: projection match, 0: projection does not match)
gs.run_command("i.landsat.import",
               input=os.path.join(homedir, "gisdata", "landsat_data"),
               flags="p")

También podemos imprimir sólo las bandas seleccionadas con un patrón.

# print a selection of bands - might be sloooow
gs.run_command("i.landsat.import",
               input=os.path.join(homedir, "gisdata", "landsat_data"), 
               pattern="B(2|3|4|5|6|8)",
               flags="p")

Para empezar a trabajar con los datos L8, vamos entonces a importar todas las bandas, recortar a la región y reproyectar al vuelo, ya que vimos que su CRS no coincide con el del proyecto donde nos encontramos.

# import all bands, subset to region and reproject
gs.run_command("i.landsat.import", 
               input=os.path.join(homedir, "gisdata", "landsat_data"), 
               extent="region",
               flags="r")

Listamos las bandas importadas y revisamos los metadatos de alguna de ellas.

# list raster maps
gs.list_grouped(type="raster")["landsat8"]

# check metadata of some imported bands
gs.raster_info(map="LC08_L1TP_229082_20200114_20200127_01_T1_B4")

Pre-procesamiento de datos satelitales

Workflow de pre-procesamiento de datos satelitales

De número digital (ND) a reflectancia y temperatura

• Los datos L8 OLI vienen en 16-bits con rango de datos entre 0 y 65535.
• i.landsat.toar convierte los números digitales (ND) en reflectancia TOA (y temperatura de brillo) para todos los sensores Landsat. Opcionalmente proporciona reflectancia de superficie (BOA) después de la corrección DOS.
• i.atcorr proporciona un método de corrección atmosférica más complejo para gran variedad de sensores (S6).

Antes de comenzar a procesar los datos, vamos a definir la resolución de la región computacional a una banda de 30 m.

# set the region to a 30m band
gs.run_command("g.region", 
               raster="LC08_L1TP_229082_20200114_20200127_01_T1_B4", 
               flags="p")

Ahora sí, convertimos los ND a reflectancia superficial y temperatura usando el método Dark Object Substraction (DOS).

# convert from DN to surface reflectance and temperature
metadata = os.path.join(homedir, "gisdata", "landsat_data", "LC08_L1TP_229082_20200114_20200127_01_T1_MTL.txt")
gs.run_command("i.landsat.toar",
               input="LC08_L1TP_229082_20200114_20200127_01_T1_B",
               output="LC08_229082_20200114_SR_B",
               sensor="oli8",
               metfile=metadata,
               method="dos1")

Corroboremos los metadatos antes y después de la conversión ND >> RS para una banda.

# list output maps
gs.list_grouped(type="raster",
                pattern="*SR*")["landsat8"]

# check info before and after for one band
print(gs.read_command("r.info", 
                      map="LC08_L1TP_229082_20200114_20200127_01_T1_B3"))

print(gs.read_command("r.info", 
                      map="LC08_229082_20200114_SR_B3"))

# Visualize results
b3_map=gj.InteractiveMap(width = 500, tiles="OpenStreetMap")
b3_map.add_raster("LC08_L1TP_229082_20200114_20200127_01_T1_B3")
b3_map.add_raster("LC08_229082_20200114_SR_B3")
b3_map.add_layer_control(position = "bottomright")
b3_map.show()

Tarea Opcional

Seguir los mismos pasos para la escena del 02/03/2020. ¿Qué notan de diferente?

Ajuste de color y composiciones RGB

Para lograr una buena visualización en composiciones RGB, primero realizamos un ajuste de colores utilizando el módulo i.colors.enhance.

# enhance the colors
gs.run_command("i.colors.enhance",
               red="LC08_229082_20200114_SR_B4",
               green="LC08_229082_20200114_SR_B3", 
               blue="LC08_229082_20200114_SR_B2",
               strength="95")

Visualicemos la combinacion RGB color natural usando gj.Map.

# display RGB
rgb_map = gj.Map(width=450, use_region=True)
rgb_map.d_rgb(red="LC08_229082_20200114_SR_B4",
              green="LC08_229082_20200114_SR_B3", 
              blue="LC08_229082_20200114_SR_B2",)
rgb_map.show()

Tarea

Hagamos una composición falso color 543. ¿Sobre qué bandas debieran realizar el ajuste?

Enmascarado de nubes con banda QA

• Landsat 8 proporciona una banda de calidad (quality assessment, QA) con valores enteros de 16 bits que representan las combinaciones de superficie, atmósfera y condiciones del sensor que pueden afectar la utilidad general de un determinado pixel.
• La extensión i.landsat.qa reclasifica la banda QA de Landsat 8 de acuerdo a la calidad del pixel.

Nota

Más información sobre la banda QA de L8 en la guía de usuario.

Vamos a utilizar i.landsat.qa para crear las reglas necesarias para identificar las nubes y sombras de nubes en las escenas L8.

# create a rule set
gs.run_command("i.landsat.qa",
               dataset="landsat_8_c1",
               cloud_shadow_confidence="Medium,High",
               cloud_confidence="Medium,High",
               output=os.path.join(homedir, "Cloud_Mask_rules.txt"))

Con las reglas que creamos anteriormente, reclasificamos la banda QA para obtener un mapa ráster que podamos luego usar como máscara.

# reclass the BQA band based on the rule set created
gs.run_command("r.reclass",
               input="LC08_L1TP_229082_20200114_20200127_01_T1_BQA",
               output="LC08_229082_20200114_Cloud_Mask",
               rules=os.path.join(homedir, "Cloud_Mask_rules.txt"))

Para estimar la superficie cubierta por nubes y sombras de nubes, usamos el módulo r.report.

# report % of clouds and shadows
print(gs.read_command("r.report",
                      map="LC08_229082_20200114_Cloud_Mask",
                      units="p",
                      flags="e"))

Visualizamos el mapa reclasificado.

# display reclassified map over RGB
rgb_map.d_rast(map="LC08_229082_20200114_Cloud_Mask")
rgb_map.show()

Fusión de datos/Pansharpening

Vamos a usar la banda pancromática (15 m) para mejorar la definición de las bandas espectrales de 30 m, por medio de: i.fusion.hpf, que aplica un método de adición basado en un filtro de paso alto. Otros métodos de pansharpening están implementados en el módulo i.pansharpen.

Como primer paso, instalamos la extensión i.fusion.hpf.

# Install the reqquired addon
gs.run_command("g.extension",
               extension="i.fusion.hpf")

Luego, como vamos a “mejorar” la resolución espacial de las bandas multiespectrales a 15 m, necesitamos alinear la resolución de la región computacional a la banda PAN.

# Set the region to PAN band (15m)
gs.run_command("g.region",
               raster="LC08_229082_20200114_SR_B8",
               flags="p")

… y ejecutamos la fusión.

# list bands
ms_bands = gs.list_grouped(type="raster", 
                           pattern="*_SR_B[1-7]")["landsat8"]
ms_bands

# Apply the fusion based on high pass filter
gs.run_command("i.fusion.hpf",
               pan="LC08_229082_20200114_SR_B8", 
               msx=ms_bands,
               suffix="hpf", 
               center="high", 
               modulation="max", 
               trim="0.0", 
               flags="lc")

Finalmente, listamos los mapas resultantes usando un patrón de búsqueda

# list the fused maps
gs.list_grouped(type="raster", 
                pattern="*.hpf")["landsat8"]

# check info of a pansharpened band
gs.raster_info("LC08_229082_20200114_SR_B4.hpf")

y visualizamos las diferencias con gj.InteractiveMap.

# display original and fused maps
hpf_map = gj.InteractiveMap(width = 500, tiles="OpenStreetMap")
hpf_map.add_raster("LC08_229082_20200114_SR_B4")
hpf_map.add_raster("LC08_229082_20200114_SR_B4.hpf")
hpf_map.add_layer_control(position = "bottomright")
hpf_map.show()

Índices de agua y vegetación

En esta sección vamos a estimar los conocidos índices de vegetación y agua. No obstante, necesitamos primero, establecer la máscara de nubes y sombras de nubes obtenida anteriormente para evitar el cómputo de los índices sobre estas áreas.

# Set the cloud mask to avoid computing over clouds
gs.run_command("r.mask",
               raster="LC08_229082_20200114_Cloud_Mask")

Calculamos el NDVI y establecemos la paleta de colores.

# Compute NDVI
ndvi_formula = "LC08_229082_20200114_NDVI = (LC08_229082_20200114_SR_B5.hpf - LC08_229082_20200114_SR_B4.hpf) / (LC08_229082_20200114_SR_B5.hpf + LC08_229082_20200114_SR_B4.hpf) * 1.0"
gs.mapcalc(exp=ndvi_formula)

# Set the color palette
gs.run_command("r.colors",
               map="LC08_229082_20200114_NDVI",
               color="ndvi")

Calculamos el NDWI y establecemos la paleta de colores.

# Compute NDWI
ndwi_formula = "LC08_229082_20200114_NDWI = (LC08_229082_20200114_SR_B5.hpf - LC08_229082_20200114_SR_B6.hpf) / (LC08_229082_20200114_SR_B5.hpf + LC08_229082_20200114_SR_B6.hpf) * 1.0"
gs.mapcalc(exp=ndwi_formula)

# Set the color palette
gs.run_command("r.colors",
               map="LC08_229082_20200114_NDWI",
               color="ndwi")

Visualizamos los mapas resultantes.

ndi = gj.InteractiveMap(width=450, use_region=True)
ndi.add_raster("LC08_229082_20200114_NDVI")
ndi.add_raster("LC08_229082_20200114_NDWI")
ndi.add_layer_control(position = "bottomright")
ndi.show()

Tarea Opcional

Estimar NDVI y NDWI para la otra escena usando el módulo i.vi.

Clasificación No Supervisada

Los pasos para realizar una clasificación no supervisada en GRASS, implican:

• Asignar semantic labels a las bandas con r.semantic.label o r.support
• Agrupar las bandas (i.e., hacer un stack): i.group
• Generar firmas para n número de clases: i.cluster
• Clasificar usando las firmas: i.maxlik

¿Qué son los semantic labels?

Los semantic labels son etiquetas que podemos agregar a cualquier mapa ráster y que nos indican qué variable está representada en ese mapa. Estas etiquetas son especialmente relevantes para las imágenes de satélite, ya que nos permiten identificar a qué sensor y banda corresponde el mapa. Son útiles a la hora de trabajar con colecciones de imágenes de satélite y también a la hora de clasificar diferentes escenas.

Por ejemplo, si generamos firmas espectrales para unas clases usando un determinado conjunto de bandas, estas firmas pueden reutilizarse para clasificar otra escena siempre que las etiquetas semánticas sean las mismas.

¡Cuidado! Aunque es posible reutilizar las firmas espectrales para cualquier escena con las mismas bandas, los cambios temporales (estaciones, impacto meteorológico) limitan su aplicabilidad sólo a escenas obtenidas más o menos al mismo tiempo.

Probemos entonces, asignar semantic labels a las bandas de L8 y los índices de vegetación.

# list the bands needed for classification
bands = gs.list_grouped(type="raster",
                        pattern="*_SR*hpf")["landsat8"]
bands

# semantic labels list
labels = ["L8_1", "L8_2", "L8_3", "L8_4", "L8_5", "L8_6", "L8_7"]

Para bandas de satélites como Landsat y Sentinel, usamos el comando r.sematic.label que toma metadatos sobre los labels de i.band.library. Es posible agregar bandas de otros satélites siempre que se siga un determinado formato.

# add semantic labels to bands
gs.run_command("r.semantic.label",
               map=bands,
               semantic_label=labels, 
               operation="add")

Para cualquier otro ráster al que deseemos agregar una etiqueta, usamos r.support que nos permite agregar distintos metadatos a nuestros datos en GRASS.

# add semantic labels to indices
for i in ["NDVI", "NDWI"]:
    gs.run_command("r.support", 
                   map=f"LC08_229082_20200114_{i}",
                   semantic_label=i)

Imprimimos los semantic labels

gs.raster_info("LC08_229082_20200114_SR_B7.hpf")["semantic_label"]

Creamos un grupo de imágenes o stack con las bandas 1 a 7 más el NDVI y el NDWI.

# add ndvi and ndwi to the list of bands
for i in ["NDVI", "NDWI"]:
    bands.append(f"LC08_229082_20200114_{i}")

bands

# create an imagery group with the list of bands
gs.run_command("i.group",
               group="l8",
               subgroup="l8",
               input=bands)

# print elements within the group 
gs.run_command("i.group",
               group="l8",
               flags="l")

Obtenemos estadísticos -firmas- para las n clases de interés con una muestra de pixeles

# statistics for unsupervised classification
gs.run_command("i.cluster",
               group="l8",
               subgroup="l8",
               sig="l8_hpf",
               classes="7",
               separation="0.6")

Usamos las firmas espectrales para realizar la clasificación no supervisada de toda la imagen

# Maximum Likelihood unsupervised classification
gs.run_command("i.maxlik",
               group="l8",
               subgroup="l8",
               sig="l8_hpf",
               output="l8_hpf_class",
               rej="l8_hpf_rej")

# Mostrar el mapa clasificado con `InteractiveMap`
clas = gj.InteractiveMap(width=450, use_region=True)
clas.add_raster("l8_hpf_class")
clas.add_layer_control(position = "bottomright")
clas.show()

Información derivada adicional

Usualmente, para realizar clasificaciones, derivamos muchas variables a partir de las bandas de sensores satelitales multiespectrales. Otra información podría obtenerse con los siguientes módulos, entre otros:

• medidas de textura: r.texture,
• medidas de diversidad: r.diversity,
• estadísticas locales con información de contexto: r.neighbors,
• transformación tasseled cap: i.tasscap,
• etc.

Info y ejercicios de clasificación en GRASS GIS

• Topic classification en los manuales de GRASS GIS
• Image classification en la wiki
• Ejemplos de clasificación en el curso dictado en Noruega en 2018
• Detección de cambios con Landsat
• Taller GRASS para sensado remoto en FOSS4G 2022