AVERTISSEMENT : Les procédures décrites dans ce tutoriel pour harmoniser les données Landsat ETM+ et OLI sont obsolètes. Elles ne sont pas recommandées ni nécessaires lorsque vous travaillez avec les données de réflectance de surface Landsat Collection 2. Pour en savoir plus, consultez la section des questions fréquentes sur l'harmonisation Landsat. Les données de la collection 1 sont obsolètes. Le code de ce tutoriel ne fonctionnera plus et ne sera plus mis à jour.
Ce tutoriel porte sur l'harmonisation de la réflectance de surface Landsat ETM+ avec la réflectance de surface Landsat OLI. Ce guide comporte les éléments suivants :
- Fonction de transformation spectrale
- Fonctions permettant de créer des données prêtes pour l'analyse
- Exemple de visualisation de série temporelle
Il s'agit d'un guide complet pour harmoniser et visualiser une série temporelle régionale de plus de 35 ans de données Landsat, qui peut être immédiatement appliquée à vos propres régions d'intérêt.
Notez que les coefficients de transformation d'ETM+ en OLI s'appliquent également à TM. Dans ce tutoriel, la référence à ETM+ est donc synonyme de TM.
À propos de Landsat
Landsat est un programme d'imagerie satellite qui collecte des images de la Terre à résolution modérée depuis 1972. Il s'agit du programme d'observation de la Terre le plus ancien basé dans l'espace. Il fournit un enregistrement temporel précieux pour identifier les tendances spatiotemporelles dans l'évolution du paysage. Les données des instruments Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) et Operational Land Imager (OLI) sont utilisées dans ce tutoriel. Elles sont étroitement liées et relativement faciles à consolider en une série temporelle cohérente qui produit un enregistrement continu de 1984 à aujourd'hui, à une cadence de 16 jours par capteur avec une résolution spatiale de 30 mètres. L'instrument Multispectral Scanner étend l'enregistrement Landsat jusqu'en 1972, mais n'est pas utilisé dans ce tutoriel. Ses données sont très différentes, ce qui rend l'intégration avec les capteurs ultérieurs difficile. Pour obtenir des exemples d'harmonisation sur l'ensemble des capteurs, consultez Savage et al. (2018) et Vogeler et al. (2018).
Pourquoi harmoniser les données ?
Roy et al. (2016) montrent qu'il existe de petites différences, mais potentiellement importantes, entre les caractéristiques spectrales de Landsat ETM+ et OLI, selon l'application. Voici quelques raisons pour lesquelles vous pouvez harmoniser des ensembles de données : produire une longue série temporelle couvrant Landsat TM, ETM+ et OLI, générer des composites intra-annuels quasi-datés pour réduire les effets des observations manquantes des lacunes ETM+ SLC-off et du masquage des nuages/ombres, ou augmenter la fréquence d'observation dans une série temporelle. Pour en savoir plus, veuillez consulter le manuscrit associé ci-dessus.
Instructions
Fonctions
Vous trouverez ci-dessous une série de fonctions nécessaires pour harmoniser ETM+ avec OLI et générer des données prêtes à l'analyse qui seront utilisées dans la section "Application" de ce tutoriel pour visualiser l'historique spectrotemporel d'un pixel.
Harmonisation
L'harmonisation est obtenue par transformation linéaire de l'espace spectral ETM+ en espace spectral OLI selon les coefficients présentés dans le tableau 2 des coefficients de régression OLS de Roy et al. (2016).
Les coefficients respectifs des bandes sont définis dans le dictionnaire suivant avec les constantes d'image de pente (slopes) et d'ordonnée à l'origine (itcps). Notez que les valeurs d'ordonnée à l'origine sont multipliées par 10 000 pour correspondre à la mise à l'échelle des données de réflectance de surface Landsat de l'USGS.
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])
};
Les noms de bandes ETM+ et OLI pour la même fenêtre de réponse spectrale ne sont pas égaux et doivent être standardisés. Les fonctions suivantes sélectionnent uniquement les bandes de réflectance et la bande pixel_qa de chaque ensemble de données, puis les renomment en fonction de la plage de longueurs d'onde qu'elles représentent.
// 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']);
}
Enfin, définissez la fonction de transformation, qui applique le modèle linéaire aux données ETM+, définit le type de données sur Int16 pour assurer la cohérence avec OLI et rattache la bande pixel_qa pour une utilisation ultérieure dans le masquage des nuages et des ombres.
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'));
}
Masquage des nuages et des ombres
Les données prêtes à l'analyse doivent être masquées pour les nuages et les ombres de nuages. La fonction suivante utilise CFmask (Zhu et al., 2015)
pixel_qa inclus avec chaque image de réflectance de surface Landsat USGS pour définir les pixels identifiés comme nuages et ombres de nuages sur "null".
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);
}
Calcul de l'indice spectral
La prochaine application utilise l'index spectral du ratio de brûlure normalisé (NBR) pour représenter l'historique spectral d'un pixel forestier touché par un incendie. Le NBR est utilisé, car Cohen et al. (2018) ont constaté que, parmi 13 indices/bandes spectraux, le NBR présentait le rapport signal/bruit le plus élevé en ce qui concerne les perturbations forestières (signal) dans l'ensemble des États-Unis. Il est calculé à l'aide de la fonction suivante.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
Dans votre propre application, vous pouvez choisir d'utiliser un autre indice spectral. C'est ici que vous pouvez modifier ou ajouter un index.
Combiner des fonctions
Définissez une fonction wrapper pour chaque ensemble de données qui regroupe toutes les fonctions ci-dessus afin de les appliquer facilement à leurs collections d'images respectives.
// 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()));
}
Dans votre application, vous pouvez choisir d'inclure ou d'exclure des fonctions. Modifiez ces fonctions si nécessaire.
Exemple de série temporelle
Les fonctions d'encapsulation prepOli et prepEtm ci-dessus peuvent être mappées sur des collections de réflectance de surface Landsat pour créer des données prêtes à l'analyse multisensorielle afin de visualiser la chronologie spectrale d'un pixel ou d'une région de pixels. Dans cet exemple, vous allez créer une série temporelle de plus de 35 ans et afficher l'historique spectral pour un seul pixel. Ce pixel particulier concerne l'historique récent d'une parcelle de forêt de conifères mature du nord-ouest du Pacifique (figure 1) qui a subi quelques perturbations dans les années 1980 et 1990, ainsi qu'un incendie de grande ampleur en 2011.
Figure 1. Emplacement et caractère du site pour l'exemple de zone d'intérêt. Forêt de conifères mature du nord-ouest du Pacifique sur le versant nord du mont Hood, Oregon, États-Unis. Images fournies par Google Earth, le service forestier de l'USDA, Landsat et Copernicus.
Définir une zone d'intérêt
Le résultat de cette application est une série temporelle d'observations Landsat pour un pixel. Un objet ee.Geometry.Point est utilisé pour définir l'emplacement du pixel.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
Dans votre application, vous pouvez sélectionner un autre pixel. Pour ce faire, remplacez les coordonnées de longitude et de latitude ci-dessus. Vous pouvez également résumer l'historique spectral d'un groupe de pixels à l'aide d'autres définitions d'objets ee.Geometry, telles que ee.Geometry.Polygon() et ee.Geometry.Rectangle().
Pour en savoir plus, consultez la section Géométries du guide du développeur.
Récupérer les collections de capteurs Landsat
Obtenez les collections de réflectance de surface Landsat USGS pour OLI, ETM+ et TM.
Cliquez sur les liens pour en savoir plus sur chaque ensemble de données.
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');
Définir un filtre de collection d'images
Définissez un filtre pour limiter les collections d'images en fonction de la limite spatiale de la zone d'intérêt, de la saison de photosynthèse maximale et de la qualité.
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)));
Modifiez ces critères de filtrage à votre guise. Identifiez d'autres propriétés à filtrer dans l'onglet "Propriétés des images" des pages de description des données dont les liens sont fournis ci-dessus.
Préparer les collections
Fusionnez les collections, appliquez le filtre et mappez la fonction prepImg sur toutes les images. L'extrait de code suivant génère un seul ee.ImageCollection contenant des images des collections de capteurs OLI, ETM+ et TM qui répondent aux critères de filtrage et qui sont traitées pour obtenir un NBR prêt à l'analyse.
// 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);
Créer un graphique de série temporelle affichant toutes les observations
L'harmonie entre les capteurs peut être évaluée de manière qualitative en traçant toutes les observations sous forme de graphique en nuage de points où chaque capteur est représenté par une couleur. La fonction ui.Chart.feature.groups d'Earth Engine fournit cet utilitaire, mais la valeur NBR observée doit d'abord être ajoutée à chaque image en tant que propriété. Mappez une fonction de réduction de région sur la collection d'images qui calcule la médiane de tous les pixels croisant la ZOI et définit le résultat comme propriété d'image.
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'));
});
Dans votre propre application, si votre zone d'intérêt est étendue, vous pouvez envisager d'augmenter scale et de spécifier les arguments bestEffort, maxPixels et tileScale pour la fonction reduceRegion afin de vous assurer que l'opération ne dépasse pas les limites maximales de pixels, de mémoire ou de délai d'expiration. Vous pouvez également remplacer le réducteur de médiane par la statistique de votre choix. Pour en savoir plus, consultez la section Statistiques d'une région d'image du guide du développeur.
La collection peut désormais être acceptée par la fonction ui.Chart.feature.groups.
La fonction attend une collection et les noms des propriétés des objets de collection à mapper au graphique. Ici, la propriété "system:time_start" sert de variable de l'axe X et "NBR" de variable de l'axe Y. Notez que la propriété NBR a été nommée par la fonction reduceRegion ci-dessus en fonction de la réduction du nom de la bande dans l'image. Si vous utilisez une autre bande (pas "NBR"), modifiez l'argument du nom de la propriété de l'axe Y en conséquence. Enfin, définissez la variable de regroupement (série) sur "SATELLITE", à laquelle des couleurs distinctes sont attribuées.
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);
Un graphique semblable à la figure 2 s'affiche dans la console après un certain temps de traitement. Veuillez noter que :
- Il y a plusieurs observations par an.
- La série temporelle est composée de trois capteurs différents.
- Aucun biais de capteur discernable.
- La fréquence des observations varie chaque année.
- La variance intra-annuelle est relativement faible et constante au fil des ans. Si vous essayez cette méthode avec NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');), vous constaterez une variabilité de la réponse intra-annuelle et interannuelle beaucoup plus importante.
- Le NBR reste élevé pour la majorité des séries temporelles, avec quelques perturbations mineures.
- Une forte et rapide diminution de la réponse NBR résulte d'un incendie de forêt visible dans
- Notez toutefois que l'incendie (Dollar Lake) s'est produit en septembre 2011. L'erreur dans l'année de détection est due au fait que la période de composite annuel s'étend de juillet à août. Les changements qui se produisent après cette période ne sont détectés qu'à la prochaine observation non nulle disponible, qui peut avoir lieu l'année suivante ou plus tard.
- La végétation commence à se rétablir (augmentation de la réponse NBR) deux ans après la détection de la perte majeure de NBR.
Figure 2 Graphique de série temporelle de réponse spectrale pour un seul pixel avec une représentation de Landsat TM, ETM+ et OLI. Les images TM et ETM+ sont harmonisées à OLI par transformation linéaire. Les images datent de juillet et août et sont filtrées pour une qualité élevée.
Créer un graphique de série temporelle affichant la médiane annuelle
Pour simplifier les séries temporelles et éliminer le bruit, vous pouvez réduire les observations intra-annuelles. Ici, la médiane est utilisée.
La première étape consiste à regrouper les images par année. Ajoutez une propriété "year" à chaque image en mappant le paramètre de collection "year" à partir de l'ee.Image.Date de chaque image.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
La nouvelle propriété "year" (année) peut être utilisée pour regrouper les images de la même année par une jointure. La jointure entre une collection d'années d'images distinctes et la collection complète fournit des listes de la même année, à partir desquelles des réductions médianes peuvent être effectuées. Sous-ensemble de la collection complète à un ensemble de représentants distincts pour chaque année.
var distinctYearCol = col.distinct('year');
Définissez un filtre et une jointure qui identifieront les années correspondantes entre la collection d'années distinctes (distinctYearCol) et la collection complète (col). Les correspondances seront enregistrées dans une propriété appelée "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Appliquez la jointure et convertissez la FeatureCollection résultante en 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}));
});
Modifiez le ee.Reducer ci-dessus si vous préférez une statistique autre que la médiane pour représenter la collection d'images d'une année donnée.
Enfin, créez un graphique qui affiche les valeurs médianes de NBR à partir d'ensembles d'observations estivales intra-annuelles. Modifiez la statistique de réduction de la région à votre convenance.
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);
Un graphique semblable à la figure 3 s'affiche dans la console après un certain temps de traitement.
Figure 3. Graphique de série temporelle de la réponse spectrale médiane estivale pour un seul pixel avec une représentation de Landsat TM, ETM+ et OLI. Les images TM et ETM+ sont harmonisées à OLI par une transformation linéaire. Les images datent de juillet et août, et sont filtrées pour une qualité élevée.
Informations supplémentaires
Autres fonctions de transformation
Le tableau 2 de Roy et al. (2016) fournit les coefficients de régression OLS et RMA pour transformer la réflectance de surface ETM+ en réflectance de surface OLI et inversement. Le tutoriel ci-dessus ne montre que la transformation ETM+ vers OLI par régression OLS. Les fonctions de toutes les options de traduction sont disponibles dans le script de l'éditeur de code ou le script source GitHub.
Incendie de Dollar Lake
L'incendie mentionné dans ce tutoriel est celui de Dollar Lake, qui s'est produit à la mi-septembre 2011 sur les pentes nord du mont Hood, dans l'Oregon (États-Unis). Pour en savoir plus, consultez cet article de blog et Varner et al. (2015).