Series de tiempo con Landsat 4 a 9 en Google Earth Engine

intermedioLección 1 de 420 minGoogle Earth Engine

Une las colecciones de Landsat 4, 5, 7, 8 y 9 y crea una mediana anual sin nubes para analizar cuatro décadas de cambios.

Introducción #

Hasta ahora trabajamos con una imagen o con la mediana de unos pocos meses. En el nivel intermedio damos un paso más: usaremos colecciones completas para ver cómo cambia un lugar a lo largo de los años.

El programa Landsat es ideal para eso. Observa la Tierra desde 1972, y desde 1982, con Landsat 4, sus sensores miden la superficie con 30 m de resolución y bandas muy parecidas entre sí. Eso permite comparar una misma zona durante más de cuatro décadas. En esta lección uniremos Landsat 4, 5, 7, 8 y 9 en una sola colección y crearemos una mediana sin nubes para cada año.

Objetivos #

Unir las colecciones de Landsat 4, 5, 7, 8 y 9 con nombres de bandas comunes.
Recorrer una lista de años y crear una mediana anual sin nubes.
Eliminar los años que no tienen imágenes.
Mostrar y exportar cada año de la serie.

Recursos #

Cinco satélites, una serie #

Satélite	Sensor	Datos en la serie	Detalle
Landsat 4	TM	1982 a 1993	Primer sensor TM, con 30 m y bandas SWIR.
Landsat 5	TM	1984 a 2012	Operó casi 29 años, un récord Guinness.
Landsat 7	ETM+	1999 a 2013	Desde 2003 sus escenas tienen franjas sin datos.
Landsat 8	OLI	2013 en adelante	Más sensibilidad radiométrica que TM y ETM+.
Landsat 9	OLI-2	2021 en adelante	Gemelo de Landsat 8. Juntos pasan cada 8 días.

Landsat 7 merece un comentario. El 31 de mayo de 2003 falló su corrector de líneas de escaneo (SLC), y desde entonces cada escena tiene franjas sin datos que suman cerca del 22 % de la imagen. En una mediana anual esas franjas se rellenan con otras fechas, así que el satélite sigue siendo útil. En el script lo usamos hasta mayo de 2013, cuando Landsat 8 empezó a entregar datos. Así cubrimos 2012, cuando Landsat 5 ya casi no tomaba imágenes y Landsat 8 todavía no estaba en órbita.

Sección 0: Crear las funciones principales #

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.

scaler() convierte los valores digitales en reflectancia de superficie, de 0 a 1, y la banda térmica en grados Kelvin. mask_clouds() oculta los píxeles marcados en la banda de calidad QA_PIXEL, que en el script se llama BQA. Cada condición atmosférica ocupa un bit de esa banda, y en la Colección 2 son estos:

Bit	Qué marca	Variable
2	Cirros	`cirrus`
3	Nubes	`cloudsBitMask`
4	Sombras de nubes	`cloudShadowBitMask`
5	Nieve	`snow`

Ojo con los bits si vienes de la Colección 1. En la antigua banda pixel_qa de la Colección 1, el bit 3 era la sombra y el 5 la nube. En QA_PIXEL de la Colección 2 el 3 es la nube, el 4 la sombra y el 5 la nieve. Si se usan los bits viejos, las sombras de nubes no se enmascaran y quedan como manchas oscuras en la mediana.

///////////////////////////////////////////////////////////////////////////////////
// 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);
}

Sección 1: Definir bandas y nombres #

Landsat 8 y 9 agregaron una banda azul costera al inicio, y por eso las demás se corrieron un número. El azul es SR_B1 en Landsat 4, 5 y 7, pero SR_B2 en Landsat 8 y 9. Para unir las colecciones renombramos las bandas con select(), que recibe una lista con las bandas de entrada y otra con los nombres nuevos. Landsat 4, 5 y 7 no tienen azul costero, así que repetimos SR_B1 como uBlue para que todas las colecciones tengan las mismas ocho bandas.

Nombre común	Landsat 4, 5 y 7	Landsat 8 y 9
uBlue	SR_B1 (repetida)	SR_B1, azul costero
Blue	SR_B1	SR_B2
Green	SR_B2	SR_B3
Red	SR_B3	SR_B4
Swir1	SR_B5	SR_B6
BQA	QA_PIXEL	QA_PIXEL
Nir	SR_B4	SR_B5
Swir2	SR_B7	SR_B7

Los sensores no son idénticos. Las bandas de OLI son más estrechas que las de TM y ETM+, así que un mismo lugar puede dar valores un poco distintos con Landsat 7 y con Landsat 8. Para ver tendencias generales no es un problema. Para análisis cuantitativos finos, Roy et al. (2016) publicaron coeficientes para armonizar ambos sensores.

///////////////////////////////////////////////////////////////////////////////////
// 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'];

Sección 2: Crear y unir las colecciones #

A cada colección le aplicamos el factor de escala y los nombres comunes, y después las unimos con merge(). Landsat 7 se limita a sus fechas útiles, como explicamos arriba.

///////////////////////////////////////////////////////////////////////////////////
// 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);

Sección 3: Definir los parámetros #

Reunimos al inicio los valores que vas a cambiar. Con startMonth y endMonth puedes quedarte con una sola temporada, por ejemplo 1 y 3 para enero a marzo, la época seca en buena parte del Caribe colombiano. Comparar siempre la misma temporada evita confundir un cambio real con la diferencia entre la época seca y la de lluvias.

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

var startYear = 1980;  // 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

Sección 4: Filtrar la colección #

Filtramos la colección unida por fechas, nubosidad sobre tierra (CLOUD_COVER_LAND), meses y área de estudio, y la ordenamos por fecha. Earth Engine solo procesa las imágenes que pasan los filtros, aunque scaler() se haya aplicado antes.

///////////////////////////////////////////////////////////////////////////////////
// 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);

Sección 5: Calcular las medianas anuales #

yearly_medians() recibe un año, toma las imágenes sin nubes de ese año y calcula su mediana. Al aplicarla a la lista de años obtenemos una mediana por año. En los primeros años del archivo puede no haber ninguna escena de tu zona, y esa mediana queda sin bandas. ee.Algorithms.If() cambia esas imágenes vacías por un 0, y removeAll([0]) las elimina de la lista. Al final convertimos la lista en una colección, colmeans.

///////////////////////////////////////////////////////////////////////////////////
// 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);

¿Y la reproyección? No reproyectamos las medianas dentro del script. Las operaciones de esta serie trabajan píxel a píxel, y Earth Engine calcula cada resultado en la escala que se le pide. La proyección y el tamaño de píxel finales se fijan al exportar, con crs y scale. Usar reproject() antes obliga a calcular todo a esa escala, incluso lo que ves en el mapa, y suele causar errores de memoria en series largas.

Sección 6: Visualizar la serie #

evaluate() trae al navegador el número de medianas, y con un ciclo normal de JavaScript agregamos una capa por año. Las capas se cargan apagadas y se numeran desde 0, en orden: actívalas una por una en el panel Layers para ver cómo cambia tu zona.

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

var m_list = img_medians;
var len_w = img_medians.size();
var visArgs_LAND = {bands: 'Red,Green,Blue', min: 0.002, max: 0.25, gamma: 1};

Map.centerObject(AOI, 10);

len_w.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(m_list.get(i));
    Map.addLayer(img, visArgs_LAND, 'RGB-' + i.toString(), false);
  }
});

Sección 7: Exportar cada año #

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.

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

/*// Eliminar comentario

var list = colmeans.toList(colmeans.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('Landsat_Median_').cat(year).getInfo(),
    folder: 'Landsat_Median',
    fileNamePrefix: ee.String('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 = colmeans.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 = 1980;  // 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 → VISUALIZACIÓN EN EL MAPA
///////////////////////////////////////////////////////////////////////////////////

var m_list = img_medians;
var len_w = img_medians.size();
var visArgs_LAND = {bands: 'Red,Green,Blue', min: 0.002, max: 0.25, gamma: 1};

Map.centerObject(AOI, 10);

len_w.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(m_list.get(i));
    Map.addLayer(img, visArgs_LAND, 'RGB-' + i.toString(), false);
  }
});

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

/*// Eliminar comentario

var list = colmeans.toList(colmeans.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('Landsat_Median_').cat(year).getInfo(),
    folder: 'Landsat_Median',
    fileNamePrefix: ee.String('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 = colmeans.size();
count.evaluate(function(n) {
  for (var i = 0; i < n; i++) {
    exportImages(i);
  }
});

*/// Eliminar comentario

Conclusiones #

Ya tienes una serie de cuatro décadas lista para analizar. En la siguiente lección construiremos la misma serie con Sentinel-2, más corta pero con 10 m de resolución. Después usaremos estas medianas anuales para calcular índices y detectar agua año a año.

Referencias #

Roy, D. P., Kovalskyy, V., Zhang, H. K., Vermote, E. F., Yan, L., Kumar, S. S. y Egorov, A. (2016). Characterization of Landsat-7 to Landsat-8 reflective wavelength and normalized difference vegetation index continuity. Remote Sensing of Environment, 185, 57–70.