Series de tiempo con Sentinel-2 en Google Earth Engine

intermedioLección 2 de 415 minGoogle Earth Engine

Crea una serie de medianas anuales sin nubes con Sentinel-2 y exporta cada año a Google Drive.

Introducción #

Sentinel-2 tiene una historia más corta que Landsat, pero mucho más detalle. La Agencia Espacial Europea (ESA) lanzó Sentinel-2A en 2015 y Sentinel-2B en 2017. Juntos fotografían cada punto de la Tierra cada cinco días, con 10 m de resolución en las bandas visibles y en el infrarrojo cercano. Es la mejor opción para seguir cambios recientes en áreas pequeñas, como una finca, una ciénaga o el borde de una ciudad.

En esta lección seguiremos la misma estructura de Series de tiempo con Landsat, sección por sección, para crear una mediana anual sin nubes con la colección de reflectancia de superficie.

Objetivos #

Adaptar el script de series de Landsat a Sentinel-2.
Crear una mediana anual sin nubes con la máscara SCL.
Mostrar y exportar cada año de la serie.
Entender por qué usamos la colección HARMONIZED en series de tiempo.

Recursos #

Dos detalles antes de empezar #

¿Desde qué año? ESA genera la reflectancia de superficie (nivel 2A) de forma sistemática para todo el planeta desde diciembre de 2018. En años anteriores muchas zonas tienen pocas imágenes o ninguna, por eso empezamos la serie en 2019. Si tu zona tiene datos anteriores, puedes adelantar startYear, y el script descartará los años vacíos.
¿Por qué HARMONIZED? Desde el 25 de enero de 2022, ESA suma 1000 a los valores de las imágenes nuevas por un cambio en su procesamiento. La colección S2_SR_HARMONIZED corrige ese desplazamiento para que las imágenes de antes y después de esa fecha sean comparables. En una serie de tiempo es indispensable: sin la corrección verías un salto falso en 2022.

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() divide las bandas entre 10 000 para pasarlas a reflectancia. Usamos addBands() con reemplazo, como en Landsat, porque así la imagen conserva su fecha; una división directa la perdería y después no podríamos filtrar por año. mask_clouds() usa la banda SCL, la clasificación de escena de ESA, para ocultar sombras, nubes y cirros.

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

// 0.1 → Función de factor de escala (valores digitales a reflectancia, 0 a 1)
function scaler (img) {
  var optical = img.select('B.*').divide(10000).toFloat();
  return img.addBands(optical, null, true);
}

// 0.2 → Función de máscara de nubes con la banda SCL
function mask_clouds (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);
}

Secciones 1 y 2: Bandas y colección #

Seleccionamos las diez bandas de superficie útiles y les damos nombres comunes. Dejamos fuera B1 y B9, que sirven para medir aerosoles y vapor de agua en la atmósfera y tienen 60 m de resolución. Las bandas de 20 m (borde rojo, B8A y SWIR) se remuestrean a 10 m al exportar.

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

var bn_s2 = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12', 'SCL'];
var bns = ['Blue', 'Green', 'Red', 'RE1', 'RE2', 'RE3', 'Nir', 'Nir2', 'Swir1', 'Swir2', 'SCL'];

///////////////////////////////////////////////////////////////////////////////////
// 2 → CREAR LA COLECCIÓN SENTINEL-2
///////////////////////////////////////////////////////////////////////////////////

var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED").map(scaler).select(bn_s2, bns);

Secciones 3 y 4: Parámetros y filtros #

Como Sentinel-2 pasa cada cinco días, un límite de 20 % de píxeles nubosos (CLOUDY_PIXEL_PERCENTAGE) todavía deja muchas imágenes por año.

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

var startYear = 2019;  // Año inicio
var endYear = 2025;    // Año final
var startMonth = 1;    // Mes inicial
var endMonth = 12;     // Mes final
var cloudcover = 20;   // Porcentaje máximo de píxeles nubosos
var AOI = AOI;         // Área de estudio

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

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

print("Imágenes disponibles para la zona de estudio (Sentinel-2)", filtered);

Sección 5: Calcular las medianas anuales #

Este bloque es idéntico al de Series de tiempo con Landsat: una mediana por año y eliminación de los años sin imágenes.

///////////////////////////////////////////////////////////////////////////////////
// 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: Visualizar la serie #

Agregamos una capa por año en color verdadero. Activa las capas una por una en el panel Layers para comparar los años.

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

var m_list = img_medians;
var len_w = img_medians.size();
var visArgs_S2 = {bands: 'Red,Green,Blue', min: 0.01, max: 0.22, gamma: 1};

Map.centerObject(AOI, 11);

len_w.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(m_list.get(i));
    Map.addLayer(img, visArgs_S2, '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. Aquí la escala es de 10 m.

///////////////////////////////////////////////////////////////////////////////////
// 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('Sentinel2_Median_').cat(year).getInfo(),
    folder: 'Sentinel_Median',
    fileNamePrefix: ee.String('Median_').cat(year).getInfo(),
    region: AOI,
    scale: 10,
    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 (valores digitales a reflectancia, 0 a 1)
function scaler (img) {
  var optical = img.select('B.*').divide(10000).toFloat();
  return img.addBands(optical, null, true);
}

// 0.2 → Función de máscara de nubes con la banda SCL
function mask_clouds (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);
}

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

var bn_s2 = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'B8', 'B8A', 'B11', 'B12', 'SCL'];
var bns = ['Blue', 'Green', 'Red', 'RE1', 'RE2', 'RE3', 'Nir', 'Nir2', 'Swir1', 'Swir2', 'SCL'];

///////////////////////////////////////////////////////////////////////////////////
// 2 → CREAR LA COLECCIÓN SENTINEL-2
///////////////////////////////////////////////////////////////////////////////////

var s2 = ee.ImageCollection("COPERNICUS/S2_SR_HARMONIZED").map(scaler).select(bn_s2, bns);

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

var startYear = 2019;  // Año inicio
var endYear = 2025;    // Año final
var startMonth = 1;    // Mes inicial
var endMonth = 12;     // Mes final
var cloudcover = 20;   // Porcentaje máximo de píxeles nubosos
var AOI = AOI;         // Área de estudio

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

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

print("Imágenes disponibles para la zona de estudio (Sentinel-2)", 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_S2 = {bands: 'Red,Green,Blue', min: 0.01, max: 0.22, gamma: 1};

Map.centerObject(AOI, 11);

len_w.evaluate(function(l) {
  for (var i = 0; i < l; i++) {
    var img = ee.Image(m_list.get(i));
    Map.addLayer(img, visArgs_S2, '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('Sentinel2_Median_').cat(year).getInfo(),
    folder: 'Sentinel_Median',
    fileNamePrefix: ee.String('Median_').cat(year).getInfo(),
    region: AOI,
    scale: 10,
    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 #

Con Sentinel-2 tienes menos años que con Landsat, pero cada uno con 10 m de detalle. Elige Landsat cuando necesites mirar décadas atrás y Sentinel-2 cuando te interesen cambios recientes en áreas pequeñas.

Referencias #

Drusch, M., Del Bello, U., Carlier, S., Colin, O., Fernandez, V., Gascon, F., Hoersch, B., Isola, C., Laberinti, P., Martimort, P., Meygret, A., Spoto, F., Sy, O., Marchese, F. y Bargellini, P. (2012). Sentinel-2: ESA's Optical High-Resolution Mission for GMES Operational Services. Remote Sensing of Environment, 120, 25–36.