Índices de vegetación con Sentinel-2 en Google Earth Engine

básicoLección 12 de 1815 minGoogle Earth Engine

Calcula, visualiza y exporta siete índices de vegetación con Sentinel-2 en Google Earth Engine.

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 Sentinel-2. Gracias a sus 10 m de resolución podremos ver diferencias entre lotes, parques o bordes de bosque que en Landsat se pierden.

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 #

Sentinel-2 Surface Reflectance (COPERNICUS/S2_SR_HARMONIZED)

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	Sentinel-2	Resolución
Azul (B)	B2	10 m
Rojo (R)	B4	10 m
NIR	B8	10 m
SWIR1	B11	20 m
SWIR2	B12	20 m

Paso 1: Preparar la imagen #

Aplicamos la máscara de nubes, que además convierte los valores en reflectancia, y tomamos la primera escena con poca nubosidad sobre el AOI.

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.

// Máscara de nubes con la banda SCL y paso a reflectancia (0 a 1)
function maskS2clouds(img) {
  var scl = img.select('SCL');
  var mask = scl.neq(3)    // sombra de nubes
    .and(scl.neq(8))       // nubes, probabilidad media
    .and(scl.neq(9))       // nubes, probabilidad alta
    .and(scl.neq(10));     // cirros
  return img.updateMask(mask).select('B.*').divide(10000);
}

// Primera imagen de Sentinel-2 con poca nubosidad sobre el AOI
var s2_image = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterDate("2023-01-01", "2023-12-01")
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 10))
  .filterBounds(AOI)
  .map(maskS2clouds)
  .first()
  .clip(AOI);
print(s2_image);

// Imagen en color verdadero
Map.addLayer(s2_image, {bands: "B4,B3,B2", min: 0, max: 0.3, gamma: 1}, "Sentinel-2");
Map.centerObject(AOI, 12);

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.

// Índices con las bandas de Sentinel-2 (B2 azul, B4 rojo, B8 NIR, B11 SWIR1, B12 SWIR2)
var NDVI = s2_image.normalizedDifference(['B8', 'B4']).rename('NDVI');
var EVI = s2_image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4'), 'BLUE': s2_image.select('B2')
}).rename('EVI');
var SAVI = s2_image.expression('((NIR - RED) / (NIR + RED + 0.5)) * 1.5', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4')
}).rename('SAVI');
var MSAVI = s2_image.expression('(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED))) / 2', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4')
}).rename('MSAVI');
var NDMI = s2_image.normalizedDifference(['B8', 'B11']).rename('NDMI');
var NBR = s2_image.normalizedDifference(['B8', 'B12']).rename('NBR');
var NBR2 = s2_image.normalizedDifference(['B11', 'B12']).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 en una sola imagen y la exportamos a 10 m.

var index = NDVI.addBands([EVI, SAVI, MSAVI, NDMI, NBR, NBR2]).toFloat();

Export.image.toDrive({
  image: index,
  description: "sentinel_2_vegetation_index",
  folder: "index",
  region: AOI,
  scale: 10,
  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.

// Máscara de nubes con la banda SCL y paso a reflectancia (0 a 1)
function maskS2clouds(img) {
  var scl = img.select('SCL');
  var mask = scl.neq(3)    // sombra de nubes
    .and(scl.neq(8))       // nubes, probabilidad media
    .and(scl.neq(9))       // nubes, probabilidad alta
    .and(scl.neq(10));     // cirros
  return img.updateMask(mask).select('B.*').divide(10000);
}

// Primera imagen de Sentinel-2 con poca nubosidad sobre el AOI
var s2_image = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED")
  .filterDate("2023-01-01", "2023-12-01")
  .filter(ee.Filter.lte('CLOUDY_PIXEL_PERCENTAGE', 10))
  .filterBounds(AOI)
  .map(maskS2clouds)
  .first()
  .clip(AOI);
print(s2_image);

// Imagen en color verdadero
Map.addLayer(s2_image, {bands: "B4,B3,B2", min: 0, max: 0.3, gamma: 1}, "Sentinel-2");
Map.centerObject(AOI, 12);

// Índices con las bandas de Sentinel-2 (B2 azul, B4 rojo, B8 NIR, B11 SWIR1, B12 SWIR2)
var NDVI = s2_image.normalizedDifference(['B8', 'B4']).rename('NDVI');
var EVI = s2_image.expression('2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4'), 'BLUE': s2_image.select('B2')
}).rename('EVI');
var SAVI = s2_image.expression('((NIR - RED) / (NIR + RED + 0.5)) * 1.5', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4')
}).rename('SAVI');
var MSAVI = s2_image.expression('(2 * NIR + 1 - sqrt(pow((2 * NIR + 1), 2) - 8 * (NIR - RED))) / 2', {
  'NIR': s2_image.select('B8'), 'RED': s2_image.select('B4')
}).rename('MSAVI');
var NDMI = s2_image.normalizedDifference(['B8', 'B11']).rename('NDMI');
var NBR = s2_image.normalizedDifference(['B8', 'B12']).rename('NBR');
var NBR2 = s2_image.normalizedDifference(['B11', 'B12']).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');

var index = NDVI.addBands([EVI, SAVI, MSAVI, NDMI, NBR, NBR2]).toFloat();

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

Conclusiones #

Sentinel-2 es ideal para seguir la vegetación a escala de finca o de ciudad. Recuerda que NDMI, NBR y NBR2 usan bandas SWIR de 20 m, así que en esos índices el detalle es algo menor que en NDVI.

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.