Detección de agua con Landsat 9 en Google Earth Engine

básicoLección 13 de 1810 minGoogle Earth Engine

Calcula NDWI, MNDWI y AWEInsh con Landsat 9 para detectar agua, y conoce el origen de cada índice.

Introducción #

El agua refleja poca luz en general. Devuelve algo en el azul y el verde, pero en el infrarrojo cercano y en el de onda corta absorbe casi todo. El suelo y la vegetación hacen lo contrario y reflejan mucho más en el infrarrojo. Si comparamos el verde con una banda infrarroja obtenemos un número que tiende a ser positivo sobre el agua y negativo sobre todo lo demás, y en esa idea se basan los índices espectrales para detectar agua.

En esta lección calcularemos tres de estos índices, NDWI, MNDWI y AWEInsh, con Landsat 9. Son la base para mapear inundaciones y extraer cuerpos de agua más adelante en el curso.

Objetivos #

Comprender en qué bandas se diferencia el agua del resto de coberturas.
Calcular NDWI, MNDWI y AWEInsh en Earth Engine.
Visualizar los índices con una paleta que resalte el agua.
Exportar los tres índices a Google Drive.

Recursos #

Los índices #

Índice	Fórmula	Landsat 8 y 9
NDWI	(Verde − NIR) / (Verde + NIR)	(B3 − B5) / (B3 + B5)
MNDWI	(Verde − SWIR1) / (Verde + SWIR1)	(B3 − B6) / (B3 + B6)
AWEInsh	4(Verde − SWIR1) − (0,25 NIR + 2,75 SWIR2)	B3, B6, B5, B7

NDWI es el más sencillo y funciona muy bien en lagos y embalses.
MNDWI cambia el NIR por el SWIR1 y separa mejor el agua de las zonas urbanas.
AWEInsh combina cuatro bandas con coeficientes calculados para que el umbral 0 separe bien el agua, incluso frente a superficies oscuras como techos o asfalto.

De dónde viene cada índice #

NDWI. McFeeters (1996) lo propuso para delimitar cuerpos de agua abiertos con Landsat, comparando el verde con el infrarrojo cercano. Ese mismo año, Gao publicó otro índice también llamado NDWI, pero pensado para medir el agua dentro de las hojas y calculado con dos bandas infrarrojas. Cuando encuentres "NDWI" en un artículo, revisa qué bandas usa.
MNDWI. Xu (2006) observó que con el NDWI las zonas urbanas también daban valores positivos y se confundían con agua. Reemplazó el infrarrojo cercano por el SWIR1, donde el agua sigue absorbiendo casi todo pero las superficies construidas reflejan mucho más.
AWEI. Feyisa, Meilby, Fensholt y Proud (2014) buscaban extraer agua de forma automática con un umbral fijo de 0. Propusieron dos versiones, AWEInsh para zonas sin sombras importantes y AWEIsh para paisajes con sombras de montañas o edificios. En el curso usamos AWEInsh.

Paso 1: Preparar la imagen #

Partimos de una escena de Landsat 9 con máscara de nubes y factor de escala.

El script usa una geometría llamada AOI. Dibújala con las herramientas de geometría del mapa en el Code Editor y renómbrala como AOI en el panel de importaciones.

// Factor de escala: convierte los valores digitales en reflectancia
function scaler(img) {
  var optical = img.select('SR_B.').multiply(0.0000275).add(-0.2).toFloat();
  var thermal = img.select('ST_B.*').multiply(0.00341802).add(149).toFloat();
  return img.addBands(optical, null, true).addBands(thermal, null, true);
}

// Máscara de nubes con la banda QA_PIXEL
function mask_clouds(img) {
  var qa = img.select('QA_PIXEL');
  var cirrus = 1 << 2;       // bit 2, cirros
  var cloud = 1 << 3;        // bit 3, nubes
  var cloudShadow = 1 << 4;  // bit 4, sombra de nubes
  var snow = 1 << 5;         // bit 5, nieve
  var mask = qa.bitwiseAnd(cirrus).eq(0)
    .and(qa.bitwiseAnd(cloud).eq(0))
    .and(qa.bitwiseAnd(cloudShadow).eq(0))
    .and(qa.bitwiseAnd(snow).eq(0));
  return img.updateMask(mask);
}

// Primera imagen de Landsat 9 con poca nubosidad sobre el AOI
var ls9_image = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
  .filterDate('2023-01-01', '2023-12-31')
  .filter(ee.Filter.lt("CLOUD_COVER_LAND", 10))
  .filterBounds(AOI)
  .map(mask_clouds)
  .map(scaler)
  .first()
  .clip(AOI);
print(ls9_image);

// Imagen en color verdadero
Map.addLayer(ls9_image, {bands: "SR_B4,SR_B3,SR_B2", min: 0.003, max: 0.19, gamma: 1}, "Landsat 9");
Map.centerObject(AOI, 8);

Paso 2: Calcular los índices para detectar agua #

NDWI y MNDWI son diferencias normalizadas entre dos bandas, así que basta con normalizedDifference(), que recibe primero la banda verde y después la infrarroja. AWEInsh combina cuatro bandas con coeficientes fijos, de modo que lo escribimos con expression() y asignamos un nombre a cada banda dentro de la fórmula.

// NDWI: (Verde - NIR) / (Verde + NIR)
var NDWI = ls9_image.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI');

// MNDWI: (Verde - SWIR1) / (Verde + SWIR1)
var MNDWI = ls9_image.normalizedDifference(['SR_B3', 'SR_B6']).rename('MNDWI');

// AWEInsh: 4 * (Verde - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)
var AWEInsh = ls9_image.expression('4 * (GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)', {
  'GREEN': ls9_image.select('SR_B3'),
  'SWIR1': ls9_image.select('SR_B6'),
  'SWIR2': ls9_image.select('SR_B7'),
  'NIR': ls9_image.select('SR_B5')
}).rename('AWEInsh');

Paso 3: Visualizar los índices #

Usamos una paleta de −1 a 1 que va del rojo al azul. Los cuerpos de agua quedan en tonos cian y azul, y el suelo y la vegetación en rojo, amarillo o verde. NDWI y MNDWI siempre caen entre −1 y 1. AWEInsh no está acotado a ese rango, pero su umbral también es 0, así que la misma paleta separa bien el agua de lo demás.

// Simbología: rojo para lo seco, azul para el agua
var waterPalette = {min: -1, max: 1, palette: ['red', 'yellow', 'green', 'cyan', 'blue']};

Map.addLayer(NDWI, waterPalette, 'NDWI');
Map.addLayer(MNDWI, waterPalette, 'MNDWI');
Map.addLayer(AWEInsh, waterPalette, 'AWEInsh');

Paso 4: Exportar los índices #

Unimos los tres índices con addBands() para que el GeoTIFF tenga una banda por índice, en el mismo orden en que aparecen en el mapa, y lo enviamos a Google Drive desde la pestaña Tasks.

// Unir los índices y exportarlos
var index = NDWI.addBands([MNDWI, AWEInsh]).toFloat();

Export.image.toDrive({
  image: index,
  description: "landsat_9_water_index",
  folder: "index",
  region: AOI,
  scale: 30,
  maxPixels: 1e13,
  fileFormat: "TIFF",
  formatOptions: {cloudOptimized: true}
});

Código completo #

Este es el script completo. Cópialo, pégalo en el Code Editor, dibuja tu AOI y presiona Run.

// Factor de escala: convierte los valores digitales en reflectancia
function scaler(img) {
  var optical = img.select('SR_B.').multiply(0.0000275).add(-0.2).toFloat();
  var thermal = img.select('ST_B.*').multiply(0.00341802).add(149).toFloat();
  return img.addBands(optical, null, true).addBands(thermal, null, true);
}

// Máscara de nubes con la banda QA_PIXEL
function mask_clouds(img) {
  var qa = img.select('QA_PIXEL');
  var cirrus = 1 << 2;       // bit 2, cirros
  var cloud = 1 << 3;        // bit 3, nubes
  var cloudShadow = 1 << 4;  // bit 4, sombra de nubes
  var snow = 1 << 5;         // bit 5, nieve
  var mask = qa.bitwiseAnd(cirrus).eq(0)
    .and(qa.bitwiseAnd(cloud).eq(0))
    .and(qa.bitwiseAnd(cloudShadow).eq(0))
    .and(qa.bitwiseAnd(snow).eq(0));
  return img.updateMask(mask);
}

// Primera imagen de Landsat 9 con poca nubosidad sobre el AOI
var ls9_image = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2")
  .filterDate('2023-01-01', '2023-12-31')
  .filter(ee.Filter.lt("CLOUD_COVER_LAND", 10))
  .filterBounds(AOI)
  .map(mask_clouds)
  .map(scaler)
  .first()
  .clip(AOI);
print(ls9_image);

// Imagen en color verdadero
Map.addLayer(ls9_image, {bands: "SR_B4,SR_B3,SR_B2", min: 0.003, max: 0.19, gamma: 1}, "Landsat 9");
Map.centerObject(AOI, 8);

// NDWI: (Verde - NIR) / (Verde + NIR)
var NDWI = ls9_image.normalizedDifference(['SR_B3', 'SR_B5']).rename('NDWI');

// MNDWI: (Verde - SWIR1) / (Verde + SWIR1)
var MNDWI = ls9_image.normalizedDifference(['SR_B3', 'SR_B6']).rename('MNDWI');

// AWEInsh: 4 * (Verde - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)
var AWEInsh = ls9_image.expression('4 * (GREEN - SWIR1) - (0.25 * NIR + 2.75 * SWIR2)', {
  'GREEN': ls9_image.select('SR_B3'),
  'SWIR1': ls9_image.select('SR_B6'),
  'SWIR2': ls9_image.select('SR_B7'),
  'NIR': ls9_image.select('SR_B5')
}).rename('AWEInsh');

// Simbología: rojo para lo seco, azul para el agua
var waterPalette = {min: -1, max: 1, palette: ['red', 'yellow', 'green', 'cyan', 'blue']};

Map.addLayer(NDWI, waterPalette, 'NDWI');
Map.addLayer(MNDWI, waterPalette, 'MNDWI');
Map.addLayer(AWEInsh, waterPalette, 'AWEInsh');

// Unir los índices y exportarlos
var index = NDWI.addBands([MNDWI, AWEInsh]).toFloat();

Export.image.toDrive({
  image: index,
  description: "landsat_9_water_index",
  folder: "index",
  region: AOI,
  scale: 30,
  maxPixels: 1e13,
  fileFormat: "TIFF",
  formatOptions: {cloudOptimized: true}
});

Conclusiones #

Los índices para detectar agua son una herramienta rápida para ubicar ríos, ciénagas y embalses. Ningún índice es perfecto en todos los paisajes, así que vale la pena comparar los tres en tu zona antes de elegir uno.

Referencias #

McFeeters, S. K. (1996). The use of the Normalized Difference Water Index (NDWI) in the delineation of open water features. International Journal of Remote Sensing, 17(7), 1425–1432.
Gao, B.-C. (1996). NDWI—A normalized difference water index for remote sensing of vegetation liquid water from space. Remote Sensing of Environment, 58(3), 257–266.
Xu, H. (2006). Modification of normalised difference water index (NDWI) to enhance open water features in remotely sensed imagery. International Journal of Remote Sensing, 27(14), 3025–3033.
Feyisa, G. L., Meilby, H., Fensholt, R. y Proud, S. R. (2014). Automated Water Extraction Index: A new technique for surface water mapping using Landsat imagery. Remote Sensing of Environment, 140, 23–35.