Índices espectrales en series de tiempo con Landsat en Google Earth Engine

intermedioLección 3 de 420 minGoogle Earth Engine

Calcula el EVI de cada año de una serie de Landsat con una función cargada desde un módulo con require().

Introducción #

Un índice calculado en una sola fecha muestra cómo está la vegetación en ese momento. Calculado año tras año, muestra tendencias: una zona que se deforesta, un humedal que se seca o un cultivo que se expande.

En esta lección calcularemos el EVI de cada mediana anual de Landsat. La función del índice viene de un módulo, un script aparte que guarda funciones para reutilizarlas en otros scripts.

Objetivos #

Cargar funciones de un módulo con require().
Aplicar una función a cada imagen de una colección con map().
Calcular el EVI para cada año de la serie.
Visualizar y exportar los índices.

Recursos #

¿Por qué el EVI? #

En bosques y cultivos muy densos, el NDVI se acerca a su máximo y deja de cambiar aunque la vegetación siga creciendo. Se dice que se satura. Huete y sus colaboradores (2002) diseñaron el EVI para reducir ese efecto: corrige la influencia del suelo y de la atmósfera con las bandas roja y azul. Lo vimos en Índices de vegetación con Landsat 9.

Secciones 0 a 5: Construir la serie anual #

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.

Las secciones 0 a 5 son las mismas de Series de tiempo con Landsat: funciones, bandas, unión de Landsat 4 a 9, parámetros, filtros y medianas anuales. Esta vez la serie empieza en 2020.

///////////////////////////////////////////////////////////////////////////////////
// 0 → CREAR FUNCIONES PRINCIPALES
///////////////////////////////////////////////////////////////////////////////////

// 0.1 → Función de factor de escala
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);
}

// 0.2 → Función de máscara de nubes
function mask_clouds (img) {
  // Banda de calidad QA_PIXEL (renombrada como BQA)
  var qa = img.select('BQA');
  // Bits de QA_PIXEL en Landsat Colección 2
  var cirrus = (1 << 2);              // bit 2: cirros
  var cloudsBitMask = (1 << 3);       // bit 3: nubes
  var cloudShadowBitMask = (1 << 4);  // bit 4: sombras de nubes
  var snow = (1 << 5);                // bit 5: nieve
  // Todos los bits deben estar en cero (condiciones despejadas)
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
    .and(qa.bitwiseAnd(cirrus).eq(0))
    .and(qa.bitwiseAnd(snow).eq(0));
  return img.updateMask(mask);
}

///////////////////////////////////////////////////////////////////////////////////
// 1 → DEFINIR BANDAS Y NOMENCLATURAS
///////////////////////////////////////////////////////////////////////////////////

var bn9 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn8 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn7 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bn5 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bns = ['uBlue', 'Blue', 'Green', 'Red', 'Swir1', 'BQA', 'Nir', 'Swir2'];

///////////////////////////////////////////////////////////////////////////////////
// 2 → CREAR Y UNIFICAR COLECCIONES LANDSAT (4, 5, 7, 8 Y 9)
///////////////////////////////////////////////////////////////////////////////////

var ls4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2").map(scaler).select(bn5, bns);
var ls5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").map(scaler).select(bn5, bns);
var ls7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2").map(scaler)
                            .filterDate('1999-04-15', '2013-05-30')
                            .select(bn7, bns);
var ls8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").map(scaler).select(bn8, bns);
var ls9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2").map(scaler).select(bn9, bns);

var merged = ls4.merge(ls5).merge(ls7).merge(ls8).merge(ls9);

///////////////////////////////////////////////////////////////////////////////////
// 3 → DEFINIR PARÁMETROS DE BÚSQUEDA
///////////////////////////////////////////////////////////////////////////////////

var startYear = 2020;  // Año inicio
var endYear = 2025;    // Año final
var startMonth = 1;    // Mes inicial
var endMonth = 12;     // Mes final
var cloudcover = 30;   // Cobertura de nubes (%)
var AOI = AOI;         // Área de estudio

///////////////////////////////////////////////////////////////////////////////////
// 4 → APLICAR FILTROS A LA COLECCIÓN FUSIONADA
///////////////////////////////////////////////////////////////////////////////////

var filtered = merged.filterDate(startYear + '-01-01', endYear + '-12-31')
                     .filter(ee.Filter.lte('CLOUD_COVER_LAND', cloudcover))
                     .filter(ee.Filter.calendarRange(startMonth, endMonth, 'month'))
                     .sort('system:time_start', false)
                     .filterBounds(AOI);

print("Imágenes disponibles para la zona de estudio (Landsat 4, 5, 7, 8 y 9)", filtered);

///////////////////////////////////////////////////////////////////////////////////
// 5 → CALCULAR MEDIANAS ANUALES
///////////////////////////////////////////////////////////////////////////////////

var years = ee.List.sequence(startYear, endYear);

function yearly_medians(year) {
  var start = ee.Date.fromYMD(year, 1, 1);
  var end = start.advance(1, 'year');
  var L = filtered.map(mask_clouds);
  var image = L.filterDate(start, end).median().clip(AOI).toFloat();
  return image.set('year', year);
}

// Generar medianas por año
var img_medians = years.map(yearly_medians);

// Remover imágenes vacías
img_medians = img_medians.map(function(ele) {
  var count = ee.Image(ele).bandNames().size();
  return ee.Algorithms.If(count.gt(0), ele, 0);
}).flatten();
img_medians = img_medians.removeAll([0]);

// Colección de medianas anuales
var colmeans = ee.ImageCollection(img_medians);
var años = ee.Array(colmeans.aggregate_array('year'));

print("Mediana de imágenes por año", colmeans);
print("Lista de años", años);

Sección 6: Calcular el índice de cada año #

El EVI viene del módulo público LIST_EDUCAINDEX, que reúne trece índices de vegetación, agua y fuego. Cada índice es una función que el módulo exporta con exports, y nuestro script carga el módulo con require() en la variable x. Esta es la función que usamos:

// Módulo LIST_EDUCAINDEX (fragmento): función EVI
exports.EVI = function(image){
  return image.expression(
    '2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
      'NIR': image.select(['Nir']),
      'RED': image.select(['Red']),
      'BLUE': image.select(['Blue'])
    }).rename('evi').toFloat();
};

Otros índices con la misma línea. Cambia x.EVI(image) por x.NDVI(image), x.NDMI(image) o cualquier otra función del módulo. Con x.calcAllIndices(image) obtienes los trece índices como bandas de una sola imagen. En la página del módulo está la lista completa, con lo que recibe y devuelve cada función.

XIndex() calcula el EVI de una imagen y le copia sus propiedades con copyProperties(), para que cada resultado conserve su año. copyProperties() devuelve un objeto genérico, por eso lo envolvemos en ee.Image(). Después aplicamos la función a toda la colección con map().

///////////////////////////////////////////////////////////////////////////////////
// 6 → CÁLCULO DE ÍNDICES ESPECTRALES
///////////////////////////////////////////////////////////////////////////////////

var x = require('users/sigeduca/routines:Indexes/LIST_EDUCAINDEX'); // MÓDULO DE EVI

// Crear función EVI
function XIndex (image){
  var INDEX = x.EVI(image);
  return ee.Image(INDEX.copyProperties(image));
}

var CalIndex = colmeans.map(XIndex);
print('Índices espectrales: ', CalIndex);

Sección 7: Visualizar el índice #

Mostramos el EVI de la primera imagen de la colección, que corresponde al primer año de la serie. Las medianas siguen el orden de la lista de años, aunque en la sección 4 hayamos ordenado las imágenes de la más reciente a la más antigua.

///////////////////////////////////////////////////////////////////////////////////
// 7 → VISUALIZACIÓN EN EL MAPA
///////////////////////////////////////////////////////////////////////////////////

Map.addLayer(CalIndex.first(), {}, 'EVI', true);
Map.centerObject(AOI, 12);

Sección 8: Exportar los índices #

La exportación viene desactivada: todo el bloque está entre /* y */, así que Earth Engine lo ignora. Para usarla, borra las dos líneas que dicen Eliminar comentario y vuelve a ejecutar. En la pestaña Tasks aparecerá una tarea por año y tendrás que ejecutar cada una. crs fija la proyección de salida: EPSG:32618 es UTM zona 18 Norte, que cubre buena parte de Colombia; cámbiala por la zona UTM de tu AOI.

///////////////////////////////////////////////////////////////////////////////////
// 8 → EXPORTACIÓN OPCIONAL A GOOGLE DRIVE
///////////////////////////////////////////////////////////////////////////////////

/*// Eliminar comentario

var list = CalIndex.toList(CalIndex.size());

// Función que exporta cada imagen de la lista
var exportImages = function(index) {
  var img = ee.Image(list.get(index));
  var year = ee.Number(img.get('year')).format('%d');

  Export.image.toDrive({
    image: img,
    description: ee.String('index_Median_').cat(year).getInfo(),
    folder: 'Landsat_Index_Median',
    fileNamePrefix: ee.String('Index_Median_').cat(year).getInfo(),
    region: AOI,
    scale: 30,
    crs: 'EPSG:32618',  // UTM 18N, cámbiala por la zona de tu AOI
    maxPixels: 1e13
  });
};

// Ejecuta una exportación por año
var count = CalIndex.size();
count.evaluate(function(n) {
  for (var i = 0; i < n; i++) {
    exportImages(i);
  }
});

*/// Eliminar comentario

Código completo #

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

///////////////////////////////////////////////////////////////////////////////////
// 0 → CREAR FUNCIONES PRINCIPALES
///////////////////////////////////////////////////////////////////////////////////

// 0.1 → Función de factor de escala
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);
}

// 0.2 → Función de máscara de nubes
function mask_clouds (img) {
  // Banda de calidad QA_PIXEL (renombrada como BQA)
  var qa = img.select('BQA');
  // Bits de QA_PIXEL en Landsat Colección 2
  var cirrus = (1 << 2);              // bit 2: cirros
  var cloudsBitMask = (1 << 3);       // bit 3: nubes
  var cloudShadowBitMask = (1 << 4);  // bit 4: sombras de nubes
  var snow = (1 << 5);                // bit 5: nieve
  // Todos los bits deben estar en cero (condiciones despejadas)
  var mask = qa.bitwiseAnd(cloudShadowBitMask).eq(0)
    .and(qa.bitwiseAnd(cloudsBitMask).eq(0))
    .and(qa.bitwiseAnd(cirrus).eq(0))
    .and(qa.bitwiseAnd(snow).eq(0));
  return img.updateMask(mask);
}

///////////////////////////////////////////////////////////////////////////////////
// 1 → DEFINIR BANDAS Y NOMENCLATURAS
///////////////////////////////////////////////////////////////////////////////////

var bn9 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn8 = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B6', 'QA_PIXEL', 'SR_B5', 'SR_B7'];
var bn7 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bn5 = ['SR_B1', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B5', 'QA_PIXEL', 'SR_B4', 'SR_B7'];
var bns = ['uBlue', 'Blue', 'Green', 'Red', 'Swir1', 'BQA', 'Nir', 'Swir2'];

///////////////////////////////////////////////////////////////////////////////////
// 2 → CREAR Y UNIFICAR COLECCIONES LANDSAT (4, 5, 7, 8 Y 9)
///////////////////////////////////////////////////////////////////////////////////

var ls4 = ee.ImageCollection("LANDSAT/LT04/C02/T1_L2").map(scaler).select(bn5, bns);
var ls5 = ee.ImageCollection("LANDSAT/LT05/C02/T1_L2").map(scaler).select(bn5, bns);
var ls7 = ee.ImageCollection("LANDSAT/LE07/C02/T1_L2").map(scaler)
                            .filterDate('1999-04-15', '2013-05-30')
                            .select(bn7, bns);
var ls8 = ee.ImageCollection("LANDSAT/LC08/C02/T1_L2").map(scaler).select(bn8, bns);
var ls9 = ee.ImageCollection("LANDSAT/LC09/C02/T1_L2").map(scaler).select(bn9, bns);

var merged = ls4.merge(ls5).merge(ls7).merge(ls8).merge(ls9);

///////////////////////////////////////////////////////////////////////////////////
// 3 → DEFINIR PARÁMETROS DE BÚSQUEDA
///////////////////////////////////////////////////////////////////////////////////

var startYear = 2020;  // Año inicio
var endYear = 2025;    // Año final
var startMonth = 1;    // Mes inicial
var endMonth = 12;     // Mes final
var cloudcover = 30;   // Cobertura de nubes (%)
var AOI = AOI;         // Área de estudio

///////////////////////////////////////////////////////////////////////////////////
// 4 → APLICAR FILTROS A LA COLECCIÓN FUSIONADA
///////////////////////////////////////////////////////////////////////////////////

var filtered = merged.filterDate(startYear + '-01-01', endYear + '-12-31')
                     .filter(ee.Filter.lte('CLOUD_COVER_LAND', cloudcover))
                     .filter(ee.Filter.calendarRange(startMonth, endMonth, 'month'))
                     .sort('system:time_start', false)
                     .filterBounds(AOI);

print("Imágenes disponibles para la zona de estudio (Landsat 4, 5, 7, 8 y 9)", filtered);

///////////////////////////////////////////////////////////////////////////////////
// 5 → CALCULAR MEDIANAS ANUALES
///////////////////////////////////////////////////////////////////////////////////

var years = ee.List.sequence(startYear, endYear);

function yearly_medians(year) {
  var start = ee.Date.fromYMD(year, 1, 1);
  var end = start.advance(1, 'year');
  var L = filtered.map(mask_clouds);
  var image = L.filterDate(start, end).median().clip(AOI).toFloat();
  return image.set('year', year);
}

// Generar medianas por año
var img_medians = years.map(yearly_medians);

// Remover imágenes vacías
img_medians = img_medians.map(function(ele) {
  var count = ee.Image(ele).bandNames().size();
  return ee.Algorithms.If(count.gt(0), ele, 0);
}).flatten();
img_medians = img_medians.removeAll([0]);

// Colección de medianas anuales
var colmeans = ee.ImageCollection(img_medians);
var años = ee.Array(colmeans.aggregate_array('year'));

print("Mediana de imágenes por año", colmeans);
print("Lista de años", años);

///////////////////////////////////////////////////////////////////////////////////
// 6 → CÁLCULO DE ÍNDICES ESPECTRALES
///////////////////////////////////////////////////////////////////////////////////

var x = require('users/sigeduca/routines:Indexes/LIST_EDUCAINDEX'); // MÓDULO DE EVI

// Crear función EVI
function XIndex (image){
  var INDEX = x.EVI(image);
  return ee.Image(INDEX.copyProperties(image));
}

var CalIndex = colmeans.map(XIndex);
print('Índices espectrales: ', CalIndex);

///////////////////////////////////////////////////////////////////////////////////
// 7 → VISUALIZACIÓN EN EL MAPA
///////////////////////////////////////////////////////////////////////////////////

Map.addLayer(CalIndex.first(), {}, 'EVI', true);
Map.centerObject(AOI, 12);

///////////////////////////////////////////////////////////////////////////////////
// 8 → EXPORTACIÓN OPCIONAL A GOOGLE DRIVE
///////////////////////////////////////////////////////////////////////////////////

/*// Eliminar comentario

var list = CalIndex.toList(CalIndex.size());

// Función que exporta cada imagen de la lista
var exportImages = function(index) {
  var img = ee.Image(list.get(index));
  var year = ee.Number(img.get('year')).format('%d');

  Export.image.toDrive({
    image: img,
    description: ee.String('index_Median_').cat(year).getInfo(),
    folder: 'Landsat_Index_Median',
    fileNamePrefix: ee.String('Index_Median_').cat(year).getInfo(),
    region: AOI,
    scale: 30,
    crs: 'EPSG:32618',  // UTM 18N, cámbiala por la zona de tu AOI
    maxPixels: 1e13
  });
};

// Ejecuta una exportación por año
var count = CalIndex.size();
count.evaluate(function(n) {
  for (var i = 0; i < n; i++) {
    exportImages(i);
  }
});

*/// Eliminar comentario

Conclusiones #

Con map() y una función, cualquier cálculo que hiciste sobre una imagen sirve para toda una serie. Y con un módulo, esa función se escribe una sola vez y se reutiliza en todos tus scripts. En la siguiente lección usaremos la misma idea para clasificar el agua de cada año con el método DSWE.

Referencias #

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.