Índices de vegetación con Landsat 9 en Google Earth Engine

básicoLección 11 de 1815 minGoogle Earth Engine

Calcula, visualiza y exporta siete índices de vegetación con Landsat 9, y conoce quién propuso cada uno y para qué.

Introducción #

Las plantas sanas absorben casi toda la luz roja para hacer fotosíntesis y reflejan con fuerza el infrarrojo cercano. Los índices espectrales de vegetación aprovechan ese contraste para resumir en un solo número qué tan vigorosa, densa o húmeda está la cobertura de cada píxel.

En esta lección calcularemos siete índices con Landsat 9, los visualizaremos con paletas adecuadas y los exportaremos juntos en un solo GeoTIFF.

Objetivos #

Calcular NDVI, EVI, SAVI, MSAVI, NDMI, NBR y NBR2.
Usar normalizedDifference() y expression() para escribir fórmulas.
Visualizar cada índice con su propia simbología.
Exportar todos los índices como bandas de una misma imagen.

Recursos #

Los índices #

Índice	Fórmula	Para qué sirve
NDVI	(NIR − R) / (NIR + R)	Vigor y densidad de la vegetación
EVI	2,5 · (NIR − R) / (NIR + 6R − 7,5B + 1)	Vegetación densa, corrige efectos de la atmósfera
SAVI	(NIR − R) / (NIR + R + 0,5) · 1,5	Zonas con suelo expuesto
MSAVI	(2NIR + 1 − √((2NIR + 1)² − 8(NIR − R))) / 2	Vegetación escasa sin fijar el factor L
NDMI	(NIR − SWIR1) / (NIR + SWIR1)	Contenido de humedad de la vegetación
NBR	(NIR − SWIR2) / (NIR + SWIR2)	Áreas quemadas
NBR2	(SWIR1 − SWIR2) / (SWIR1 + SWIR2)	Severidad de quemas y humedad del suelo

De dónde viene cada índice #

Cada índice nació para resolver un problema concreto. Conocer esa historia ayuda a elegir el adecuado para tu caso.

NDVI. Rouse, Haas, Schell y Deering lo presentaron en 1974 para seguir la vegetación de las Grandes Llanuras de Estados Unidos con el satélite ERTS-1, el que después se llamó Landsat 1. Unos años más tarde, Tucker (1979) mostró su relación con la biomasa verde y lo convirtió en el índice más usado del mundo.
SAVI. Huete (1988) notó que, donde la vegetación es escasa, el NDVI cambia con el color y el brillo del suelo. Para reducir ese efecto añadió un factor de ajuste L, que suele fijarse en 0,5.
MSAVI. Qi y colaboradores (1994) propusieron que ese factor L se calculara a partir de la propia imagen en lugar de fijarlo a mano. La fórmula que usamos corresponde a la segunda versión de su trabajo, conocida también como MSAVI2.
EVI. Huete y colaboradores (2002) lo afinaron para el sensor MODIS de la NASA. Usa la banda azul para corregir el efecto de la atmósfera y no se satura tan rápido como el NDVI en bosques muy densos.
NDMI. Cambia el rojo por el infrarrojo de onda corta (SWIR1), una banda sensible al agua que contienen las hojas. Hardisky, Klemas y Smart lo propusieron en 1983 con el nombre de NDII, y Wilson y Sader (2002) lo llamaron NDMI al usarlo con Landsat para detectar talas de bosque.
NBR. El fuego reduce el reflejo en el infrarrojo cercano, porque desaparece la vegetación sana, y lo aumenta en el SWIR2 por el suelo y la ceniza expuestos. López García y Caselles (1991) aprovecharon ese contraste para cartografiar quemas con Landsat, y Key y Benson lo estandarizaron en el protocolo FIREMON del Servicio Forestal de Estados Unidos, donde se compara antes y después del incendio (dNBR).
NBR2. Es una variante con las dos bandas SWIR que el USGS incluye en sus productos de índices de Landsat. Es sensible a la humedad de la vegetación y se usa para seguir la recuperación después de un incendio.

Banda	Landsat 8 y 9
Azul (B)	SR_B2
Rojo (R)	SR_B4
NIR	SR_B5
SWIR1	SR_B6
SWIR2	SR_B7

Paso 1: Preparar la imagen #

Reutilizamos las funciones de factor de escala y máscara de nubes. Filtramos Landsat 9 por fecha, región y nubosidad, y tomamos la primera escena con first(). Los índices deben calcularse siempre sobre reflectancia, nunca sobre los valores digitales.

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, 10);

Paso 2: Calcular los índices #

Para los índices que son una diferencia normalizada entre dos bandas, normalizedDifference() es lo más corto. Para fórmulas con constantes o raíces usamos expression(), que asigna un nombre a cada banda dentro de la ecuación.

// NDVI: (NIR - Rojo) / (NIR + Rojo)
var NDVI = ls9_image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

// EVI: 2.5 * (NIR - Rojo) / (NIR + 6 * Rojo - 7.5 * Azul + 1)
var EVI = ls9_image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4'),
  'BLUE': ls9_image.select('SR_B2')
}).rename('EVI');

// SAVI: (NIR - Rojo) / (NIR + Rojo + 0.5) * 1.5
var SAVI = ls9_image.expression('((NIR - RED) / (NIR + RED + 0.5)) * 1.5', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4')
}).rename('SAVI');

// MSAVI: (2 * NIR + 1 - raíz((2 * NIR + 1)^2 - 8 * (NIR - Rojo))) / 2
var MSAVI = ls9_image.expression('(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED))) / 2', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4')
}).rename('MSAVI');

// NDMI: (NIR - SWIR1) / (NIR + SWIR1)
var NDMI = ls9_image.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDMI');

// NBR: (NIR - SWIR2) / (NIR + SWIR2)
var NBR = ls9_image.normalizedDifference(['SR_B5', 'SR_B7']).rename('NBR');

// NBR2: (SWIR1 - SWIR2) / (SWIR1 + SWIR2)
var NBR2 = ls9_image.normalizedDifference(['SR_B6', 'SR_B7']).rename('NBR2');

Paso 3: Visualizar los índices #

Los índices de vegetación van de 0 (sin vegetación) a 1 (vegetación muy densa), así que comparten una paleta que pasa del café al verde oscuro. NDMI y NBR tienen paletas propias.

// Simbologías
var vegPalette = {min: 0, max: 1, palette: ['0000ff', 'DF923D', 'F1B555', 'FCD163', '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01',
  '011D01', '011301']};
var ndmiParams = {min: -0.5, max: 0.5, palette: ["ff0000", "ffa500", "ffff00", "00ff00", "187a07"]};
var nbrParams = {min: 0, max: 0.5, palette: ["e47f44", "ffe483", "5d61ed", "2902ff"]};

Map.addLayer(NDVI, vegPalette, 'NDVI');
Map.addLayer(EVI, vegPalette, 'EVI');
Map.addLayer(SAVI, vegPalette, 'SAVI');
Map.addLayer(MSAVI, vegPalette, 'MSAVI');
Map.addLayer(NDMI, ndmiParams, 'NDMI');
Map.addLayer(NBR, nbrParams, 'NBR');
Map.addLayer(NBR2, nbrParams, 'NBR2');

Paso 4: Exportar los índices #

Unimos los siete índices con addBands() para que el GeoTIFF tenga una banda por índice, y lo enviamos a Google Drive.

// Unir los índices en una sola imagen y exportarla
var index = NDVI.addBands([EVI, SAVI, MSAVI, NDMI, NBR, NBR2]).toFloat();

Export.image.toDrive({
  image: index,
  description: "landsat_9_vegetation_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, 10);

// NDVI: (NIR - Rojo) / (NIR + Rojo)
var NDVI = ls9_image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');

// EVI: 2.5 * (NIR - Rojo) / (NIR + 6 * Rojo - 7.5 * Azul + 1)
var EVI = ls9_image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4'),
  'BLUE': ls9_image.select('SR_B2')
}).rename('EVI');

// SAVI: (NIR - Rojo) / (NIR + Rojo + 0.5) * 1.5
var SAVI = ls9_image.expression('((NIR - RED) / (NIR + RED + 0.5)) * 1.5', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4')
}).rename('SAVI');

// MSAVI: (2 * NIR + 1 - raíz((2 * NIR + 1)^2 - 8 * (NIR - Rojo))) / 2
var MSAVI = ls9_image.expression('(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED))) / 2', {
  'NIR': ls9_image.select('SR_B5'),
  'RED': ls9_image.select('SR_B4')
}).rename('MSAVI');

// NDMI: (NIR - SWIR1) / (NIR + SWIR1)
var NDMI = ls9_image.normalizedDifference(['SR_B5', 'SR_B6']).rename('NDMI');

// NBR: (NIR - SWIR2) / (NIR + SWIR2)
var NBR = ls9_image.normalizedDifference(['SR_B5', 'SR_B7']).rename('NBR');

// NBR2: (SWIR1 - SWIR2) / (SWIR1 + SWIR2)
var NBR2 = ls9_image.normalizedDifference(['SR_B6', 'SR_B7']).rename('NBR2');

// Simbologías
var vegPalette = {min: 0, max: 1, palette: ['0000ff', 'DF923D', 'F1B555', 'FCD163', '99B718',
  '74A901', '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01', '012E01',
  '011D01', '011301']};
var ndmiParams = {min: -0.5, max: 0.5, palette: ["ff0000", "ffa500", "ffff00", "00ff00", "187a07"]};
var nbrParams = {min: 0, max: 0.5, palette: ["e47f44", "ffe483", "5d61ed", "2902ff"]};

Map.addLayer(NDVI, vegPalette, 'NDVI');
Map.addLayer(EVI, vegPalette, 'EVI');
Map.addLayer(SAVI, vegPalette, 'SAVI');
Map.addLayer(MSAVI, vegPalette, 'MSAVI');
Map.addLayer(NDMI, ndmiParams, 'NDMI');
Map.addLayer(NBR, nbrParams, 'NBR');
Map.addLayer(NBR2, nbrParams, 'NBR2');

// Unir los índices en una sola imagen y exportarla
var index = NDVI.addBands([EVI, SAVI, MSAVI, NDMI, NBR, NBR2]).toFloat();

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

Conclusiones #

Con unas pocas líneas tienes siete índices listos para analizar el estado de la vegetación, la humedad o el efecto de un incendio. Prueba a calcularlos sobre una composición sin nubes en lugar de una sola escena, y compara dos fechas para ver cómo cambia tu zona.

Referencias #

Rouse, J. W., Haas, R. H., Schell, J. A. y Deering, D. W. (1974). Monitoring vegetation systems in the Great Plains with ERTS. Third ERTS Symposium, NASA SP-351, 309–317.
Tucker, C. J. (1979). Red and photographic infrared linear combinations for monitoring vegetation. Remote Sensing of Environment, 8(2), 127–150.
Hardisky, M. A., Klemas, V. y Smart, R. M. (1983). The influence of soil salinity, growth form, and leaf moisture on the spectral radiance of Spartina alterniflora canopies. Photogrammetric Engineering and Remote Sensing, 49(1), 77–83.
Huete, A. R. (1988). A soil-adjusted vegetation index (SAVI). Remote Sensing of Environment, 25(3), 295–309.
López García, M. J. y Caselles, V. (1991). Mapping burns and natural reforestation using Thematic Mapper data. Geocarto International, 6(1), 31–37.
Qi, J., Chehbouni, A., Huete, A. R., Kerr, Y. H. y Sorooshian, S. (1994). A modified soil adjusted vegetation index. Remote Sensing of Environment, 48(2), 119–126.
Huete, A., Didan, K., Miura, T., Rodriguez, E. P., Gao, X. y Ferreira, L. G. (2002). Overview of the radiometric and biophysical performance of the MODIS vegetation indices. Remote Sensing of Environment, 83(1–2), 195–213.
Wilson, E. H. y Sader, S. A. (2002). Detection of forest harvest type using multiple dates of Landsat TM imagery. Remote Sensing of Environment, 80(3), 385–396.
Key, C. H. y Benson, N. C. (2006). Landscape assessment: ground measure of severity, the Composite Burn Index; and remote sensing of severity, the Normalized Burn Ratio. USDA Forest Service, RMRS-GTR-164-CD.