SIG EDUCA
YouTube Tutoriales
cursos / google earth engine intermedio

Agua superficial con DSWE (Jones, 2019) en Google Earth Engine

intermedioLección 4 de 430 minGoogle Earth Engine

Clasifica el agua de cada año con el método DSWE, calcula su área, grafica cómo cambia y crea una animación GIF de la serie de Landsat.

Introducción #

En el curso básico detectamos agua con un índice y un umbral. Eso funciona bien en lagos y ríos abiertos, pero falla en las orillas, en los humedales con vegetación y en los píxeles que solo están cubiertos de agua en parte. DSWE (Dynamic Surface Water Extent) aborda ese problema combinando cinco pruebas espectrales en lugar de una.

En esta lección aplicaremos DSWE a cada año de la serie de Landsat. Con la máscara de agua de cada año mediremos cuántos kilómetros cuadrados de agua hay, graficaremos cómo cambian y crearemos una animación GIF de la serie. Es un solo script, de la sección 0 a la 11.

Objetivos #

Recursos #

De dónde viene DSWE #

DSWE fue desarrollado por John W. Jones, del Servicio Geológico de Estados Unidos (USGS). La primera versión (Jones, 2015) nació para seguir el agua en los humedales de los Everglades, en Florida, donde buena parte de la inundación ocurre bajo la vegetación y es difícil de ver desde el satélite. Jones comparó sus resultados con los niveles de agua que mide en campo la red EDEN (Everglades Depth Estimation Network). En 2019 publicó una revisión de las dos pruebas de agua parcial, que mejoró la detección de píxeles cubiertos solo en parte por agua (Jones, 2019).

El USGS distribuye DSWE como producto oficial de Landsat, pero solo para Estados Unidos. Para cualquier otro lugar, como Colombia, hay que calcularlo. Nuestro módulo parte de la implementación de RivWidthCloud, un algoritmo para medir el ancho de los ríos publicado por Yang, Pavelsky, Allen y Donchyts (2020), adaptada a Landsat Colección 2.

Las cinco pruebas #

DSWE trabaja con reflectancia de superficie, de 0 a 1, y usa cinco valores: MNDWI, NDVI, MBSRV (verde + rojo), MBSRN (NIR + SWIR1) y AWEsh. AWEsh es la versión para zonas con sombras del índice de Feyisa et al. (2014) que vimos en el curso básico.

PruebaCondiciónQué busca
1MNDWI > 0,124Agua abierta según el MNDWI
2MBSRV > MBSRNEl visible (verde + rojo) supera al infrarrojo (NIR + SWIR1), algo típico del agua
3AWEsh > 0Agua según el AWEsh
4MNDWI > −0,44, SWIR1 < 0,09, NIR < 0,15 y NDVI < 0,7Agua parcial, con criterio conservador
5MNDWI > −0,5, Azul < 0,1, SWIR1 < 0,3, SWIR2 < 0,1 y NIR < 0,25Agua parcial, con criterio agresivo

Las tres primeras pruebas buscan agua abierta. Las dos últimas buscan agua parcial y son las que revisó Jones en 2019: la prueba 4 es estricta para no confundir el agua con bosques oscuros, y la 5 es más permisiva para no perder el agua que hay bajo pastos y juncos.

De las pruebas a las clases #

Cada prueba aporta un dígito a un código de cinco cifras, ordenado de la prueba 5 a la 1. Por ejemplo, 00111 significa que se cumplieron las pruebas 1, 2 y 3. El USGS asigna cada uno de los 32 códigos posibles a una clase:

ClaseNombreCuándo se asigna
0No aguaSe cumple una prueba o ninguna, salvo si la única es la 5
1Agua, confianza altaSe cumplen cuatro o cinco pruebas
2Agua, confianza moderadaSe cumplen tres pruebas
3Agua parcial, conservadoraSolo se cumplen las pruebas 4 y 5 (código 11000)
4Agua parcial, agresivaSe cumplen dos pruebas, salvo 4 y 5 juntas, o solo la 5

En la práctica, las clases 1 y 2 son agua abierta, y las 3 y 4 son zonas que pueden estar inundadas en parte, como humedales u orillas. La función ClassifyWaterJones2019() del módulo devuelve una máscara de agua (waterMask) con las clases 1 y 2.

DSWE se diseñó para escenas individuales. Al aplicarlo a una mediana anual obtenemos el agua típica del año, no la máxima: un desborde que duró dos semanas puede no aparecer. Si te interesa la temporada de lluvias, ajusta startMonth y endMonth para que la mediana use solo esos meses.

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 2000.

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

El módulo Jones_2019R #

Las pruebas y la tabla de clases están en el módulo público users/sigeduca/routines:Indexes/Jones_2019R, que el script carga con require(). No tienes que copiarlo; lo mostramos para que veas cómo funciona, y también está en la página Jones_2019R. La R indica la versión revisada para Landsat Colección 2, nivel 2. Cada prueba devuelve una imagen de ceros y unos. Al multiplicarlas por 1, 10, 100, 1000 y 10 000 y sumarlas obtenemos el código de cinco cifras de cada píxel; los ceros a la izquierda desaparecen, así que 00111 queda como 111. Después, cada grupo de códigos se convierte en su clase.

Agua oscura y valores negativos. El factor de escala de la Colección 2 resta 0,2, así que en el mar o en lagos profundos algunas bandas quedan un poco por debajo de 0. normalizedDifference() enmascara cualquier píxel con un valor negativo, y sin corrección esas zonas quedarían sin clasificar, sobre todo en las costas. Por eso, al inicio de Dswe(), llevamos a 0 las bandas negativas con max(0).

///////////////////////////////////////////////////////////////////////////////////
// MÓDULO JONES_2019R: AGUA SUPERFICIAL CON DSWE PARA LANDSAT COLECCIÓN 2, NIVEL 2
///////////////////////////////////////////////////////////////////////////////////

// Referencia principal:
// Jones, J. W. (2019). Improved automated detection of subpixel-scale inundation—Revised
// dynamic surface water extent (DSWE) partial surface water tests. Remote Sensing, 11(4), 374.

/* Código tomado de RivWidthCloud (Yang et al., 2020):
github.com/seanyx/RivWidthCloudPaper/blob/master/RivWidthCloud_JavaScript/functions_Landsat578/functions_waterClassification_Jones2019.js

ADAPTACIÓN A LANDSAT COLECCIÓN 2:
El original usaba Landsat Colección 1 SR, guardada x 10 000, por eso los umbrales de las
pruebas 4 y 5 eran 900, 1500, 1000, 3000 y 2500. Aquí la imagen llega escalada con scaler()
(reflectancia de 0 a 1), así que esos umbrales pasan a 0.09, 0.15, 0.1, 0.3 y 0.25,
los mismos que publica el USGS para la Colección 2.
Además, las bandas negativas se llevan a 0 al inicio de Dswe(), porque normalizedDifference()
enmascara los píxeles negativos y el agua oscura de las costas quedaba sin clasificar.
*/

// Funciones de las pruebas espectrales de DSWE
var Mndwi = function(image) {
  return(image.normalizedDifference(['Green', 'Swir1']).rename('mndwi'));
};
var Mbsrv = function(image) {
  return(image.select(['Green']).add(image.select(['Red'])).rename('mbsrv'));
};
var Mbsrn = function(image) {
  return(image.select(['Nir']).add(image.select(['Swir1'])).rename('mbsrn'));
};
var Ndvi = function(image) {
  return(image.normalizedDifference(['Nir', 'Red']).rename('ndvi'));
};
var Awesh = function(image) {
  return(image
  .expression('Blue + 2.5 * Green + (-1.5) * mbsrn + (-0.25) * Swir2', {
    'Blue': image.select(['Blue']),
    'Green': image.select(['Green']),
    'mbsrn': Mbsrn(image).select(['mbsrn']),
    'Swir2': image.select(['Swir2'])
  }));
};

var Dswe = function(i) {
  // En agua oscura (mar, lagos profundos) la reflectancia escalada puede ser negativa,
  // y normalizedDifference() enmascara esos píxeles. Se usa 0 como mínimo.
  i = i.select(['Blue', 'Green', 'Red', 'Nir', 'Swir1', 'Swir2']).max(0);
  var mndwi = Mndwi(i);
  var mbsrv = Mbsrv(i);
  var mbsrn = Mbsrn(i);
  var awesh = Awesh(i);
  var swir1 = i.select(['Swir1']);
  var nir = i.select(['Nir']);
  var ndvi = Ndvi(i);
  var blue = i.select(['Blue']);
  var swir2 = i.select(['Swir2']);
  
  
  var t1 = mndwi.gt(0.124);
  var t2 = mbsrv.gt(mbsrn);
  var t3 = awesh.gt(0);
  var t4 = mndwi.gt(-0.44)
  .and(swir1.lt(0.09))   // Colección 1: 900
  .and(nir.lt(0.15))     // Colección 1: 1500
  .and(ndvi.lt(0.7));
  var t5 = mndwi.gt(-0.5)
  .and(blue.lt(0.1))     // Colección 1: 1000
  .and(swir1.lt(0.3))    // Colección 1: 3000
  .and(swir2.lt(0.1))    // Colección 1: 1000
  .and(nir.lt(0.25));    // Colección 1: 2500
  
  var t = t1.add(t2.multiply(10)).add(t3.multiply(100)).add(t4.multiply(1000)).add(t5.multiply(10000));
  
  var noWater = t.eq(0)
  .or(t.eq(1))
  .or(t.eq(10))
  .or(t.eq(100))
  .or(t.eq(1000));
  var hWater = t.eq(1111)
  .or(t.eq(10111))
  .or(t.eq(11011))
  .or(t.eq(11101))
  .or(t.eq(11110))
  .or(t.eq(11111));
  var mWater = t.eq(111)
  .or(t.eq(1011))
  .or(t.eq(1101))
  .or(t.eq(1110))
  .or(t.eq(10011))
  .or(t.eq(10101))
  .or(t.eq(10110))
  .or(t.eq(11001))
  .or(t.eq(11010))
  .or(t.eq(11100));
  var pWetland = t.eq(11000);
  var lWater = t.eq(11)
  .or(t.eq(101))
  .or(t.eq(110))
  .or(t.eq(1001))
  .or(t.eq(1010))
  .or(t.eq(1100))
  .or(t.eq(10000))
  .or(t.eq(10001))
  .or(t.eq(10010))
  .or(t.eq(10100));
  
  var iDswe = noWater.multiply(0)
  .add(hWater.multiply(1))
  .add(mWater.multiply(2))
  .add(pWetland.multiply(3))
  .add(lWater.multiply(4));
  
  return(iDswe.rename('dswe'));
};

// Máscara de agua: clases 1 (confianza alta) y 2 (confianza moderada)
exports.ClassifyWaterJones2019 = function(img) {
  var dswe = Dswe(img);
  var waterMask = dswe.eq(1).or(dswe.eq(2)).rename(['waterMask']);
  return(waterMask);
};

Ojo con los umbrales si usas otro código. En RivWidthCloud, las pruebas 4 y 5 usan valores como 900, 1500 o 3000, porque se escribieron para Landsat Colección 1, que guardaba la reflectancia multiplicada por 10 000. Nuestras imágenes ya pasaron por scaler() y están entre 0 y 1, así que esos umbrales equivalen a 0,09, 0,15 y 0,3, los mismos que publica el USGS para la Colección 2. Si se dejan los valores de la Colección 1 con datos escalados, las condiciones de reflectancia se cumplen siempre, y DSWE marca como agua más píxeles de los que debería.

Sección 6: Aplicar DSWE a la serie #

Jones() agrega la máscara de agua como banda nueva, así que cada año conserva su reflectancia y su año. Después seleccionamos solo waterMask y guardamos esas máscaras en masked.

///////////////////////////////////////////////////////////////////////////////////
// 6 → APLICAR DSWE JONES (2019)
///////////////////////////////////////////////////////////////////////////////////

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

// Crear función Jones (DSWE)
function Jones (image){
  var river = x.ClassifyWaterJones2019(image);
  return image.addBands(river);
}

var water_jones_2019 = colmeans.map(Jones);
var masked = water_jones_2019.select("waterMask"); // colección de máscaras de agua

print("Máscara de agua, Jones (2019)", masked);

Sección 7: Calcular el área de agua #

Un mapa de agua por año es útil, pero la pregunta que suele importar es más directa: ¿hay más o menos agua que antes? ee.Image.pixelArea() da el área de cada píxel en metros cuadrados, y la dividimos entre un millón para pasarla a km². Multiplicada por la máscara de agua, solo suman los píxeles de agua. reduceRegion() suma ese valor dentro del AOI, y el resultado se guarda como propiedad water_area de cada imagen. Con bestEffort: true, Earth Engine aumenta el tamaño de píxel si el AOI es muy grande, para que el cálculo no falle.

///////////////////////////////////////////////////////////////////////////////////
// 7 → CALCULAR ÁREA DE CUERPOS DE AGUA
///////////////////////////////////////////////////////////////////////////////////

// Función para calcular el área de agua de cada imagen
function water_area (img){
  var water = img.eq(1);
  var area = ee.Image.pixelArea();
  var waterArea = water.multiply(area).divide(1000000).rename('water_area');
  var stats = waterArea.reduceRegion({reducer: ee.Reducer.sum(), geometry: AOI, scale: 30,
                                      bestEffort: true, maxPixels: 10000000000});
  return img.set(stats).copyProperties(img);
}

var water_area_c = masked.map(water_area);
var year_si = ee.Array(masked.aggregate_array('year'));
var wa_y = ee.Array(water_area_c.aggregate_array('water_area'));

Sección 8: Graficar la evolución #

aggregate_array() reúne los años y las áreas en dos listas, y ui.Chart.array.values() las dibuja con los años en el eje horizontal.

///////////////////////////////////////////////////////////////////////////////////
// 8 → VISUALIZACIÓN DE RESULTADOS EN LA CONSOLA
///////////////////////////////////////////////////////////////////////////////////

// Gráfico del cambio de área en la consola
var chart_W = ui.Chart.array.values({array: wa_y, axis: 0, xLabels: year_si})
  .setOptions({
    title: 'EVOLUCIÓN DEL ÁREA DE AGUA',
    colors: ['cf513e'],
    hAxis: {title: 'Años', titleTextStyle: {italic: false, bold: true}},
    vAxis: {title: 'Área (km²)', titleTextStyle: {italic: false, bold: true}},
    lineSize: 3,
    pointSize: 3,
    dataOpacity: 0.4,
    legend: {position: 'none'}
  });

print(chart_W);

Cuidado con los huecos. Las nubes que la mediana no logró quitar dejan píxeles sin datos, y en los años de Landsat 7 aparecen además sus franjas vacías. Esos píxeles no cuentan ni como agua ni como tierra, así que un año con muchos huecos mostrará menos agua de la real. Antes de interpretar una caída en el gráfico, activa la máscara de ese año en el mapa y revisa si tiene zonas vacías.

Sección 9: Visualizar la serie #

Para cada año agregamos dos capas apagadas: la imagen en color verdadero y la máscara de agua, en la que el agua aparece en blanco y lo demás en negro.

///////////////////////////////////////////////////////////////////////////////////
// 9 → VISUALIZACIÓN EN EL MAPA
///////////////////////////////////////////////////////////////////////////////////
// Visualización RGB

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

// Visualización de la máscara de agua de Jones
var list_masked = masked.toList(masked.size());
var len_masked = list_masked.size();

len_masked.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(list_masked.get(i));
    Map.addLayer(img, {}, 'Máscara Jones - ' + i.toString(), false);
  }
});

Sección 10: Crear la animación GIF #

Un gráfico resume el cambio en números, pero una animación lo muestra. Cada cuadro del GIF debe ser una imagen RGB de 8 bits, y visualize() entrega justo eso. El año lo escribe el paquete text de Gennadii Donchyts, muy usado en la comunidad de Earth Engine y disponible para todos.

gifParams fija el tamaño máximo en píxeles (dimensions) y la velocidad (framesPerSecond). En annotations indicamos que el texto vaya arriba y que lo tome de la propiedad año. El tamaño de la letra depende de la escala del mapa cuando ejecutas el script; si sale muy pequeña o muy grande, cambia el 2 de Map.getScale() * 2.

makeGifAnimation() recibe una colección, unos parámetros de visualización y un nombre. Convierte cada imagen en un cuadro, escribe el año y une la serie con su versión invertida para que la animación vaya y vuelva sin cortes. La llamamos dos veces: el falso color SWIR1-NIR-Rojo muestra el agua casi negra, la vegetación verde y el suelo desnudo en tonos rosados; la máscara de agua muestra solo el agua, en azul.

///////////////////////////////////////////////////////////////////////////////////
// 10 → CREACIÓN DE ANIMACIÓN CON FORMATO GIF
///////////////////////////////////////////////////////////////////////////////////

// Parámetros:
//   collection → colección de imágenes a animar
//   visArgs    → argumentos de visualización (bandas, min, max, palette, etc.)
//   label      → nombre para identificar la animación en la consola

var text = require('users/gena/packages:text');

var gifParams = {
  'region': AOI,
  'dimensions': 800,
  'crs': 'EPSG:4326',
  'framesPerSecond': 6
};

var annotations = [{
  position: 'top',
  offset: '5%',
  margin: '5%',
  property: 'año',
  scale: Map.getScale() * 2
}];

function makeGifAnimation(collection, visArgs, label) {

  // Convertir cada imagen a visualización RGB y añadir la etiqueta del año
  var visFun = function(img) {
    var y = ee.Number(img.get('year')).toInt();
    return img.visualize(visArgs).set('año', ee.String(y));
  };

  // Escribir el año sobre cada imagen
  var text_f = function(img) {
    return text.annotateImage(img, {}, AOI, annotations)
               .copyProperties(img, img.propertyNames());
  };

  // Aplicar visualización y anotación, y construir la animación (ida y vuelta)
  var animated = collection.map(visFun).map(text_f);
  var basicAnimation = animated.merge(animated.sort('año', false));

  // Mostrar en la consola
  print('Enlace del GIF → ' + label, basicAnimation.getVideoThumbURL(gifParams));
  print('Animación → ' + label, ui.Thumbnail(basicAnimation, gifParams));

  return basicAnimation;
}

// Animación en falso color (Swir1, Nir, Red)
var visArgs_RGB = {bands: 'Swir1,Nir,Red', min: 0.0125, max: 0.41448, gamma: 1};
var anim_RGB = makeGifAnimation(ee.ImageCollection(img_medians), visArgs_RGB, 'Landsat falso color');

// Animación de la máscara de agua (Jones, 2019)
var visArgs_WATER = {bands: 'waterMask', opacity: 1, palette: ['0c0c0c', '0400ff']};
var anim_WATER = makeGifAnimation(masked, visArgs_WATER, 'Máscara de agua Jones (2019)');

Para descargar el GIF, haz clic en el enlace que aparece en la consola y guárdalo desde el navegador. Si Earth Engine responde que la animación es demasiado grande, reduce dimensions o el tamaño del AOI.

Sección 11: Exportar las máscaras #

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.

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

/*// Eliminar comentario

var list = masked.toList(masked.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('Jones_Median_').cat(year).getInfo(),
    folder: 'Jones_Median',
    fileNamePrefix: ee.String('Jones_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 = masked.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 = 2000;  // 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 → APLICAR DSWE JONES (2019)
///////////////////////////////////////////////////////////////////////////////////

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

// Crear función Jones (DSWE)
function Jones (image){
  var river = x.ClassifyWaterJones2019(image);
  return image.addBands(river);
}

var water_jones_2019 = colmeans.map(Jones);
var masked = water_jones_2019.select("waterMask"); // colección de máscaras de agua

print("Máscara de agua, Jones (2019)", masked);

///////////////////////////////////////////////////////////////////////////////////
// 7 → CALCULAR ÁREA DE CUERPOS DE AGUA
///////////////////////////////////////////////////////////////////////////////////

// Función para calcular el área de agua de cada imagen
function water_area (img){
  var water = img.eq(1);
  var area = ee.Image.pixelArea();
  var waterArea = water.multiply(area).divide(1000000).rename('water_area');
  var stats = waterArea.reduceRegion({reducer: ee.Reducer.sum(), geometry: AOI, scale: 30,
                                      bestEffort: true, maxPixels: 10000000000});
  return img.set(stats).copyProperties(img);
}

var water_area_c = masked.map(water_area);
var year_si = ee.Array(masked.aggregate_array('year'));
var wa_y = ee.Array(water_area_c.aggregate_array('water_area'));

///////////////////////////////////////////////////////////////////////////////////
// 8 → VISUALIZACIÓN DE RESULTADOS EN LA CONSOLA
///////////////////////////////////////////////////////////////////////////////////

// Gráfico del cambio de área en la consola
var chart_W = ui.Chart.array.values({array: wa_y, axis: 0, xLabels: year_si})
  .setOptions({
    title: 'EVOLUCIÓN DEL ÁREA DE AGUA',
    colors: ['cf513e'],
    hAxis: {title: 'Años', titleTextStyle: {italic: false, bold: true}},
    vAxis: {title: 'Área (km²)', titleTextStyle: {italic: false, bold: true}},
    lineSize: 3,
    pointSize: 3,
    dataOpacity: 0.4,
    legend: {position: 'none'}
  });

print(chart_W);

///////////////////////////////////////////////////////////////////////////////////
// 9 → VISUALIZACIÓN EN EL MAPA
///////////////////////////////////////////////////////////////////////////////////
// Visualización RGB

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

// Visualización de la máscara de agua de Jones
var list_masked = masked.toList(masked.size());
var len_masked = list_masked.size();

len_masked.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(list_masked.get(i));
    Map.addLayer(img, {}, 'Máscara Jones - ' + i.toString(), false);
  }
});

///////////////////////////////////////////////////////////////////////////////////
// 10 → CREACIÓN DE ANIMACIÓN CON FORMATO GIF
///////////////////////////////////////////////////////////////////////////////////

// Parámetros:
//   collection → colección de imágenes a animar
//   visArgs    → argumentos de visualización (bandas, min, max, palette, etc.)
//   label      → nombre para identificar la animación en la consola

var text = require('users/gena/packages:text');

var gifParams = {
  'region': AOI,
  'dimensions': 800,
  'crs': 'EPSG:4326',
  'framesPerSecond': 6
};

var annotations = [{
  position: 'top',
  offset: '5%',
  margin: '5%',
  property: 'año',
  scale: Map.getScale() * 2
}];

function makeGifAnimation(collection, visArgs, label) {

  // Convertir cada imagen a visualización RGB y añadir la etiqueta del año
  var visFun = function(img) {
    var y = ee.Number(img.get('year')).toInt();
    return img.visualize(visArgs).set('año', ee.String(y));
  };

  // Escribir el año sobre cada imagen
  var text_f = function(img) {
    return text.annotateImage(img, {}, AOI, annotations)
               .copyProperties(img, img.propertyNames());
  };

  // Aplicar visualización y anotación, y construir la animación (ida y vuelta)
  var animated = collection.map(visFun).map(text_f);
  var basicAnimation = animated.merge(animated.sort('año', false));

  // Mostrar en la consola
  print('Enlace del GIF → ' + label, basicAnimation.getVideoThumbURL(gifParams));
  print('Animación → ' + label, ui.Thumbnail(basicAnimation, gifParams));

  return basicAnimation;
}

// Animación en falso color (Swir1, Nir, Red)
var visArgs_RGB = {bands: 'Swir1,Nir,Red', min: 0.0125, max: 0.41448, gamma: 1};
var anim_RGB = makeGifAnimation(ee.ImageCollection(img_medians), visArgs_RGB, 'Landsat falso color');

// Animación de la máscara de agua (Jones, 2019)
var visArgs_WATER = {bands: 'waterMask', opacity: 1, palette: ['0c0c0c', '0400ff']};
var anim_WATER = makeGifAnimation(masked, visArgs_WATER, 'Máscara de agua Jones (2019)');

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

/*// Eliminar comentario

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

*/// Eliminar comentario

Conclusiones #

Con un solo script recorriste todo el análisis: construiste la serie, clasificaste el agua de cada año con DSWE, mediste su área y la mostraste en una animación. Antes de sacar conclusiones, compara siempre la misma temporada y revisa los años con huecos. La función makeGifAnimation() sirve para cualquier colección, así que puedes usarla también con la serie de Sentinel-2 o con los índices de vegetación.

Si trabajas con Sentinel-2, el proyecto OPERA del JPL de la NASA publica DSWx-HLS, un producto basado en DSWE que combina imágenes de Landsat y Sentinel-2. También está en el catálogo de Earth Engine como OPERA/DSWX/L3_V1/HLS.

Referencias #

Volver al curso← Google Earth Engine intermedio SiguienteTodos los cursos →
© 2026 SIG EDUCA