ADVERTENCIA: Los procedimientos descritos en este instructivo para armonizar los datos de Landsat ETM+ y OLI están desactualizados y no se recomiendan ni son necesarios cuando se trabaja con datos de reflectancia de la superficie de la Colección 2 de Landsat. Para obtener más información, consulta la entrada de Preguntas frecuentes sobre la armonización de Landsat. Los datos de la colección 1 dejaron de estar disponibles, y el código de este instructivo ya no funcionará ni se actualizará.
En este instructivo, se explica cómo armonizar la reflectancia de la superficie de Landsat ETM+ con la reflectancia de la superficie de Landsat OLI. En ella encontrarás las siguientes herramientas:
- Una función de transformación espectral
- Funciones para crear datos listos para el análisis
- Ejemplo de visualización de una serie temporal
Esta guía integral tiene como objetivo armonizar y visualizar una serie temporal regional de más de 35 años de datos de Landsat que se puede aplicar de inmediato a tus propias regiones de interés.
Ten en cuenta que los coeficientes para transformar ETM+ en OLI también se aplican a TM. Por lo tanto, en este instructivo, la referencia a ETM+ es sinónimo de TM.
Acerca de Landsat
Landsat es un programa de imágenes satelitales que recopila imágenes de la Tierra con resolución moderada desde 1972. Como el programa de observación de la Tierra más largo basado en el espacio, proporciona un valioso registro temporal para identificar las tendencias espacio-temporales en el cambio del paisaje. En este instructivo, se usan datos de los instrumentos Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) y Operational Land Imager (OLI). Están estrechamente relacionados y son relativamente fáciles de consolidar en una serie temporal coherente que produce un registro continuo desde 1984 hasta el presente, con una cadencia de 16 días por sensor y una resolución espacial de 30 metros. El instrumento Multispectral Scanner extiende el registro de Landsat hasta 1972, pero no se usa en este instructivo. Sus datos son bastante diferentes, lo que dificulta la integración con sensores posteriores. Consulta Savage et al. (2018) y Vogeler et al. (2018) para ver ejemplos de armonización en todos los sensores.
Por qué es importante la armonización
Roy et al. (2016) demuestran que existen diferencias pequeñas, pero potencialmente significativas, entre las características espectrales de Landsat ETM+ y OLI, según la aplicación. Entre los motivos por los que podrías querer armonizar conjuntos de datos, se incluyen los siguientes: producir una serie temporal extensa que abarque Landsat TM, ETM+ y OLI, generar composiciones intraanuales cercanas a la fecha para reducir los efectos de las observaciones faltantes de las brechas de SLC desactivado de ETM+ y el enmascaramiento de nubes o sombras, o aumentar la frecuencia de observación dentro de una serie temporal. Consulta el manuscrito vinculado más arriba para obtener más información.
Instrucciones
Funciones
A continuación, se muestra una serie de funciones necesarias para armonizar los datos de ETM+ con los de OLI y generar datos listos para el análisis que se usarán en la sección de la aplicación de este instructivo para visualizar el historial espectrotemporal de un píxel.
Armonización
La armonización se logra a través de la transformación lineal del espacio espectral de ETM+ al espacio espectral de OLI según los coeficientes que se presentan en la tabla 2 de los coeficientes de regresión de OLS de Roy et al. (2016).
Los coeficientes específicos de cada banda se definen en el siguiente diccionario con constantes de imagen de pendiente (slopes) y de intersección (itcps). Ten en cuenta que los valores de la intersección con el eje Y se multiplican por 10,000 para que coincidan con la escala de los datos de reflectancia de la superficie de USGS Landsat.
var coefficients = {
itcps: ee.Image.constant([0.0003, 0.0088, 0.0061, 0.0412, 0.0254, 0.0172])
.multiply(10000),
slopes: ee.Image.constant([0.8474, 0.8483, 0.9047, 0.8462, 0.8937, 0.9071])
};
Los nombres de las bandas del ETM+ y el OLI para la misma ventana de respuesta espectral no son iguales y deben estandarizarse. Las siguientes funciones seleccionan solo las bandas de reflectancia y la banda pixel_qa de cada conjunto de datos, y les cambian el nombre según el rango de longitud de onda que representan.
// Function to get and rename bands of interest from OLI.
function renameOli(img) {
return img.select(
['B2', 'B3', 'B4', 'B5', 'B6', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
// Function to get and rename bands of interest from ETM+.
function renameEtm(img) {
return img.select(
['B1', 'B2', 'B3', 'B4', 'B5', 'B7', 'pixel_qa'],
['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2', 'pixel_qa']);
}
Por último, define la función de transformación, que aplica el modelo lineal a los datos de ETM+, convierte el tipo de datos en Int16 para que sea coherente con OLI y vuelve a adjuntar la banda pixel_qa para usarla más adelante en el enmascaramiento de nubes y sombras.
function etmToOli(img) {
return img.select(['Blue', 'Green', 'Red', 'NIR', 'SWIR1', 'SWIR2'])
.multiply(coefficients.slopes)
.add(coefficients.itcps)
.round()
.toShort()
.addBands(img.select('pixel_qa'));
}
Enmascaramiento de nubes y sombras
Los datos listos para el análisis deben tener las nubes y las sombras de las nubes enmascaradas. La siguiente función usa la máscara de CF (Zhu et al., 2015). La banda pixel_qa se incluye con cada imagen de reflectancia de la superficie de Landsat USGS para establecer en nulo los píxeles identificados como nubes y sombras de nubes.
function fmask(img) {
var cloudShadowBitMask = 1 << 3;
var cloudsBitMask = 1 << 5;
var qa = img.select('pixel_qa');
var mask = qa.bitwiseAnd(cloudShadowBitMask)
.eq(0)
.and(qa.bitwiseAnd(cloudsBitMask).eq(0));
return img.updateMask(mask);
}
Cálculo del índice espectral
La próxima aplicación usa el índice espectral de la relación de quemado normalizada (NBR) para representar el historial espectral de un píxel boscoso afectado por un incendio forestal. Se usa el NBR porque Cohen et al. (2018) descubrieron que, entre 13 índices o bandas espectrales, el NBR tenía la mayor relación señal-ruido con respecto a las perturbaciones forestales (señal) en todo EE.UU. Se calcula con la siguiente función.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
En tu propia aplicación, puedes optar por usar un índice espectral diferente. Aquí es donde se puede modificar o agregar un índice adicional.
Cómo combinar funciones
Define una función de wrapper para cada conjunto de datos que consolide todas las funciones anteriores para facilitar su aplicación a sus respectivas colecciones de imágenes.
// Define function to prepare OLI images.
function prepOli(img) {
var orig = img;
img = renameOli(img);
img = fmask(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
// Define function to prepare ETM+ images.
function prepEtm(img) {
var orig = img;
img = renameEtm(img);
img = fmask(img);
img = etmToOli(img);
img = calcNbr(img);
return ee.Image(img.copyProperties(orig, orig.propertyNames()));
}
En tu aplicación, puedes decidir incluir o excluir funciones. Modifica estas funciones según sea necesario.
Ejemplo de serie temporal
Las funciones de wrapper prepOli y prepEtm anteriores se pueden asignar a las colecciones de reflectancia de superficie de Landsat para crear datos aptos para el análisis entre sensores y visualizar la cronología espectral de un píxel o una región de píxeles. En este ejemplo, crearás una serie temporal de más de 35 años y mostrarás el historial espectral de un solo píxel. Este píxel en particular relaciona el historial reciente de un parche de bosque de coníferas maduro del noroeste del Pacífico (figura 1) que experimentó algunas perturbaciones en las décadas de 1980 y 1990, y una quema de gran magnitud en 2011.
Figura 1: Ubicación y carácter del sitio para el ejemplo de área de interés. Bosque de coníferas maduro del noroeste del Pacífico en la ladera norte del monte Hood, Oregón, EE.UU. Imágenes cortesía de: Google Earth, USDA Forest Service, Landsat y Copernicus.
Cómo definir un área de interés (AOI)
El resultado de esta aplicación es una serie temporal de observaciones de Landsat para un píxel. Se usa un objeto ee.Geometry.Point para definir la ubicación del píxel.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
En tu aplicación, es posible que selecciones un píxel diferente. Para ello, reemplaza las coordenadas de longitud y latitud anteriores. Como alternativa, puedes resumir el historial espectral de un grupo de píxeles con otras definiciones de objetos ee.Geometry, como ee.Geometry.Polygon() y ee.Geometry.Rectangle().
Consulta la sección Geometries de la Guía para desarrolladores para obtener más información.
Recupera colecciones de sensores de Landsat
Obtén colecciones de reflectancia de superficie de Landsat USGS para OLI, ETM+ y TM.
Visita los vínculos para obtener más información sobre cada conjunto de datos.
var oliCol = ee.ImageCollection('LANDSAT/LC08/C01/T1_SR');
var etmCol = ee.ImageCollection('LANDSAT/LE07/C01/T1_SR');
var tmCol = ee.ImageCollection('LANDSAT/LT05/C01/T1_SR');
Define un filtro de colección de imágenes
Define un filtro para restringir las colecciones de imágenes según el límite espacial del área de interés, la temporada de mayor fotosíntesis y la calidad.
var colFilter = ee.Filter.and(
ee.Filter.bounds(aoi), ee.Filter.calendarRange(182, 244, 'day_of_year'),
ee.Filter.lt('CLOUD_COVER', 50), ee.Filter.lt('GEOMETRIC_RMSE_MODEL', 10),
ee.Filter.or(
ee.Filter.eq('IMAGE_QUALITY', 9),
ee.Filter.eq('IMAGE_QUALITY_OLI', 9)));
Modifica estos criterios de filtrado como desees. Identifica otras propiedades por las que puedes filtrar en la pestaña "Propiedades de las imágenes" de las páginas de descripción de datos vinculadas anteriormente.
Prepara las colecciones
Combina las colecciones, aplica el filtro y asigna la función prepImg a todas las imágenes. El resultado del siguiente fragmento es un solo ee.ImageCollection que contiene imágenes de las colecciones de sensores OLI, ETM+ y TM que cumplen con los criterios de filtro y se procesan para obtener el NBR listo para el análisis.
// Filter collections and prepare them for merging.
oliCol = oliCol.filter(colFilter).map(prepOli);
etmCol = etmCol.filter(colFilter).map(prepEtm);
tmCol = tmCol.filter(colFilter).map(prepEtm);
// Merge the collections.
var col = oliCol.merge(etmCol).merge(tmCol);
Crea un gráfico de serie temporal que muestre todas las observaciones
La armonía entre los sensores se puede evaluar de forma cualitativa trazando todas las observaciones como un diagrama de dispersión en el que el sensor se distingue por el color. La función ui.Chart.feature.groups de Earth Engine proporciona esta utilidad, pero primero se debe agregar el valor del NBR observado a cada imagen como una propiedad. Aplica una función de reducción de mapa a la colección de imágenes que calcula la mediana de todos los píxeles que se intersecan con el AOI y establece el resultado como una propiedad de la imagen.
var allObs = col.map(function(img) {
var obs = img.reduceRegion(
{geometry: aoi, reducer: ee.Reducer.median(), scale: 30});
return img.set('NBR', obs.get('NBR'));
});
En tu propia aplicación, si tu AOI es un área grande, puedes considerar aumentar el valor de scale y especificar los argumentos bestEffort, maxPixels y tileScale para la función reduceRegion, de modo que la operación no supere los límites máximos de píxeles, memoria o tiempo de espera. También puedes reemplazar el reductor de mediana por la estadística que prefieras. Consulta la sección Estadísticas de una región de imagen de la Guía para desarrolladores para obtener más información.
Ahora, la función ui.Chart.feature.groups puede aceptar la colección.
La función espera una colección y los nombres de las propiedades del objeto de colección para asignar al gráfico. Aquí, la propiedad "system:time_start" funciona como la variable del eje X y "NBR" como la variable del eje Y. Ten en cuenta que la función reduceRegion anterior asignó el nombre de la propiedad NBR según el nombre de la banda en la imagen que se reduce. Si usas una banda diferente (que no sea "NBR"), cambia el argumento del nombre de la propiedad del eje Y según corresponda. Por último, establece la variable de agrupación (serie) en "SATELLITE", a la que se asignan colores distintos.
var chartAllObs =
ui.Chart.feature.groups(allObs, 'system:time_start', 'NBR', 'SATELLITE')
.setChartType('ScatterChart')
.setSeriesNames(['TM', 'ETM+', 'OLI'])
.setOptions({
title: 'All Observations',
colors: ['f8766d', '00ba38', '619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
pointSize: 6,
dataOpacity: 0.5
});
print(chartAllObs);
Después de un tiempo de procesamiento, aparecerá en la consola un gráfico similar al de la Figura 2. Ten en cuenta los aspectos siguientes:
- Hay varias observaciones por año.
- La serie temporal se compone de tres sensores diferentes.
- No hay sesgos discernibles en los sensores.
- La frecuencia de observación varía anualmente.
- Existe cierta varianza intraanual, que es constante en todos los años y relativamente ajustada. Si se prueba con el NDVI (img.normalizedDifference(["SWIR1", "NIR"]).rename("NDVI");), se observa una variabilidad de respuesta intraanual e interanual mucho mayor.
- El NBR sigue siendo alto para la mayoría de las series temporales, con algunas perturbaciones menores.
- Una gran y rápida reducción en la respuesta del NBR se debe a un incendio forestal evidente en
- Sin embargo, ten en cuenta que el incendio (Dollar Lake) ocurrió en septiembre de 2011. El error en el año de detección se debe a que el período compuesto anual abarca de julio a agosto. Los cambios que se producen después de este período no se detectan hasta la siguiente observación disponible que no sea nula, que podría ser el año siguiente o más adelante.
- La recuperación de la vegetación (aumento de la respuesta del NBR) comienza dos años después de que se detecta la pérdida importante del NBR.
Figura 2: Gráfico de series temporales de respuesta espectral para un solo píxel con representación de Landsat TM, ETM+ y OLI. Las imágenes de TM y ETM+ se armonizan con OLI mediante una transformación lineal. Las imágenes son de julio y agosto, y se filtran para garantizar su alta calidad.
Crea un gráfico de serie temporal que muestre la mediana anual
Para simplificar y quitar el ruido de las series temporales, se puede aplicar la reducción de observaciones intraanuales. Aquí se usa la mediana.
El primer paso es agrupar las imágenes por año. Agrega una nueva propiedad "year" a cada imagen asignando la propiedad "year" de la configuración de la colección de cada ee.Image.Date de la imagen.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
La nueva propiedad "año" se puede usar para agrupar imágenes del mismo año con una unión. La unión entre una colección de años de imágenes distintos y la colección completa proporciona listas del mismo año, a partir de las cuales se pueden realizar reducciones de la mediana. Crea un subconjunto de la colección completa con un conjunto de representantes de años distintos.
var distinctYearCol = col.distinct('year');
Define un filtro y una unión que identifiquen los años coincidentes entre la colección de años distintos (distinctYearCol) y la colección completa (col). Las coincidencias se guardarán en una propiedad llamada "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Aplica la unión y convierte el FeatureCollection resultante en un ImageCollection.
var joinCol = ee.ImageCollection(join.apply(distinctYearCol, col, filter));
// Apply median reduction among matching year collections.
var medianComp = joinCol.map(function(img) {
var yearCol = ee.ImageCollection.fromImages(img.get('year_matches'));
return yearCol.reduce(ee.Reducer.median())
.set('system:time_start', img.date().update({month: 8, day: 1}));
});
Cambia el ee.Reducer anterior si prefieres una estadística que no sea la mediana para representar la colección de imágenes de un año determinado.
Por último, crea un gráfico que muestre los valores de la mediana del NBR de conjuntos de observaciones de verano anuales. Cambia la estadística de reducción de la región como quieras.
var chartMedianComp = ui.Chart.image
.series({
imageCollection: medianComp,
region: aoi,
reducer: ee.Reducer.median(),
scale: 30,
xProperty: 'system:time_start',
})
.setSeriesNames(['NBR Median'])
.setOptions({
title: 'Intra-annual Median',
colors: ['619cff'],
hAxis: {title: 'Date'},
vAxis: {title: 'NBR'},
lineWidth: 6
});
print(chartMedianComp);
Después de un tiempo de procesamiento, aparecerá en la consola un gráfico similar al de la figura 3.
Figura 3: Gráfico de series temporales de la respuesta espectral mediana del verano para un solo píxel con representación de Landsat TM, ETM+ y OLI. Las imágenes de TM y ETM+ se armonizan con OLI a través de una transformación lineal. Las imágenes son de julio y agosto, y se filtran para garantizar su alta calidad.
Información adicional
Funciones de transformación alternativas
En la tabla 2 de Roy et al. (2016), se proporcionan los coeficientes de regresión de MCO y RMA para transformar la reflectancia de la superficie del ETM+ en la reflectancia de la superficie del OLI y viceversa. En el instructivo anterior, solo se muestra la transformación de ETM+ a OLI por medio de la regresión de OLS. Las funciones para todas las opciones de traducción se pueden encontrar en la secuencia de comandos del Editor de código o en la secuencia de comandos fuente de GitHub.
El incendio de Dollar Lake
El incendio que se menciona en este instructivo es el de Dollar Lake, que ocurrió a mediados de septiembre de 2011 en las laderas norte del monte Hood, Oregón, EE.UU. Visita esta entrada de blog y consulta Varner et al. (2015) para obtener más información.