ACHTUNG: Die in dieser Anleitung beschriebenen Verfahren zum Abstimmen von Landsat ETM+- und OLI-Daten sind veraltet und werden bei der Arbeit mit Landsat Collection 2-Daten zur Oberflächenreflektanz nicht empfohlen oder sind nicht erforderlich. Weitere Informationen finden Sie im FAQ-Eintrag zur Landsat-Harmonisierung. Die Daten aus Collection 1 sind veraltet und der Code in dieser Anleitung funktioniert nicht mehr und wird auch nicht mehr aktualisiert.
In dieser Anleitung geht es darum, die Oberflächenreflexion von Landsat ETM+ an die Oberflächenreflexion von Landsat OLI anzupassen. Er enthält Folgendes:
- Eine spektrale Transformationsfunktion
- Funktionen zum Erstellen von bereitstellbaren Daten
- Beispiel für die Visualisierung einer Zeitreihe
Dieser Leitfaden soll Ihnen helfen, eine über 35 Jahre umfassende regionale Zeitreihe von Landsat-Daten zu harmonisieren und zu visualisieren, die Sie sofort auf Ihre eigenen Regionen anwenden können.
Die Koeffizienten für die Transformation von ETM+ zu OLI gelten auch für TM. Daher ist in dieser Anleitung der Verweis auf ETM+ gleichbedeutend mit TM.
Landsat
Landsat ist ein Satellitenbildprogramm, das seit 1972 Bilder der Erde mit mittlerer Auflösung aufnimmt. Als das am längsten laufende raumbasierte Erdbeobachtungsprogramm bietet es einen wertvollen zeitlichen Datensatz zur Identifizierung raumzeitlicher Trends bei Landschaftsveränderungen. In dieser Anleitung werden Daten der Instrumente Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) und Operational Land Imager (OLI) verwendet. Sie sind eng miteinander verbunden und lassen sich relativ einfach zu einer konsistenten Zeitreihe zusammenfassen, die einen kontinuierlichen Datensatz von 1984 bis heute mit einer Häufigkeit von 16 Tagen pro Sensor und einer räumlichen Auflösung von 30 Metern ergibt. Das Multispectral Scanner-Instrument (MSS) erweitert die Landsat-Aufzeichnungen bis ins Jahr 1972, wird in diesem Tutorial jedoch nicht verwendet. Die Daten sind sehr unterschiedlich, was die Integration mit späteren Sensoren erschwert. Beispiele für die Harmonisierung über alle Sensoren hinweg finden Sie bei Savage et al. (2018) und Vogeler et al. (2018).
Warum Harmonisierung?
Roy et al. (2016) zeigen, dass es je nach Anwendung kleine, aber potenziell signifikante Unterschiede zwischen den spektralen Eigenschaften von Landsat ETM+ und OLI gibt. Gründe für die Harmonisierung von Datensätzen sind beispielsweise die Erstellung einer langen Zeitreihe, die Landsat TM, ETM+ und OLI umfasst, die Generierung von intra-annuellen Composites aus dem gleichen Zeitraum, um die Auswirkungen fehlender Beobachtungen aus ETM+ SLC-off-Lücken und Cloud-/Schattenmaskierung zu reduzieren, oder die Erhöhung der Beobachtungsfrequenz innerhalb einer Zeitreihe. Weitere Informationen finden Sie im oben verlinkten Manuskript.
Anleitung
Funktionen
Im Folgenden finden Sie eine Reihe von Funktionen, die erforderlich sind, um ETM+ an OLI anzupassen und analysierbare Daten zu generieren, die im Anwendungsabschnitt dieses Tutorials verwendet werden, um die spektral-zeitliche Entwicklung eines Pixels zu visualisieren.
Harmonisierung
Die Harmonisierung erfolgt durch lineare Transformation des ETM+-Spektralraums in den OLI-Spektralraum gemäß den in Roy et al. (2016), Tabelle 2, angegebenen OLS-Regressionskoeffizienten.
Die bandbezogenen Koeffizienten werden im folgenden Dictionary mit den Bildkonstanten für die Steigung (slopes) und den Achsenabschnitt (itcps) definiert. Die Werte für den y-Achsenabschnitt werden mit 10.000 multipliziert,um der Skalierung der USGS-Landsat-Daten zur Oberflächenreflexion zu entsprechen.
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])
};
Die ETM+- und OLI-Bandnamen für dasselbe spektrale Reaktionsfenster sind nicht gleich und müssen standardisiert werden. Mit den folgenden Funktionen werden nur Reflexionsbänder und das pixel_qa-Band aus jedem Dataset ausgewählt und entsprechend dem Wellenlängenbereich, den sie darstellen, umbenannt.
// 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']);
}
Definieren Sie schließlich die Transformationsfunktion, die das lineare Modell auf ETM+-Daten anwendet, den Datentyp zur Konsistenz mit OLI als Int16 festlegt und das pixel_qa-Band zur späteren Verwendung bei der Maskierung von Wolken und Schatten wieder anhängt.
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'));
}
Wolken- und Schattenmaskierung
Bei analysereifen Daten sollten Wolken und Wolkenschatten maskiert sein. Die folgende Funktion verwendet die CFmask (Zhu et al., 2015)pixel_qa-Band, das in jedem Landsat-Bild mit USGS-Oberflächenreflexion enthalten ist, um Pixel, die als Wolken und Wolkenschatten identifiziert wurden, auf „null“ zu setzen.
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);
}
Berechnung des Spektralindex
In der kommenden Anwendung wird der spektrale Index „Normalized Burn Ratio“ (NBR) verwendet, um die spektrale Entwicklung eines bewaldeten Pixels darzustellen, das von einem Waldbrand betroffen ist. NBR wird verwendet, weil Cohen et al. (2018) herausgefunden haben, dass NBR unter 13 Spektralindizes/-bändern das größte Signal-Rausch-Verhältnis in Bezug auf Waldstörungen (Signal) in den gesamten USA aufwies. Sie wird mit der folgenden Funktion berechnet.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
In Ihrer eigenen Anwendung können Sie einen anderen Spektralindex verwenden. Hier können Sie einen zusätzlichen Index ändern oder hinzufügen.
Funktionen kombinieren
Definieren Sie für jedes Dataset eine Wrapper-Funktion, die alle oben genannten Funktionen zusammenfasst, damit sie bequem auf die jeweiligen Bildsammlungen angewendet werden können.
// 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()));
}
In Ihrer Anwendung können Sie Funktionen ein- oder ausschließen. Ändern Sie diese Funktionen nach Bedarf.
Beispiel für eine Zeitreihe
Die Wrapper-Funktionen prepOli und prepEtm oben können auf Landsat-Sammlungen für die Oberflächenreflektanz angewendet werden, um sensorübergreifende analysierbare Daten zu erstellen, mit denen die spektrale Chronologie eines Pixels oder einer Region von Pixeln visualisiert werden kann. In diesem Beispiel erstellen Sie eine Zeitreihe von mehr als 35 Jahren und zeigen den spektralen Verlauf für einen einzelnen Pixel an. Dieses spezielle Pixel bezieht sich auf die jüngste Geschichte eines reifen Nadelwaldes im pazifischen Nordwesten (Abbildung 1), der in den 1980er- und 1990er-Jahren einigen Störungen ausgesetzt war und 2011 einen Brand mit hoher Intensität erlitt.
Abbildung 1. Standort- und Standortmerkmale, z. B. Interessenbereich. Reifer Nadelwald im pazifischen Nordwesten am Nordhang des Mount Hood, Oregon, USA. Bilder mit freundlicher Genehmigung von: Google Earth, USDA Forest Service, Landsat und Copernicus.
Interessengebiet definieren
Das Ergebnis dieser Anwendung ist eine Zeitreihe von Landsat-Beobachtungen für ein Pixel. Mit einem ee.Geometry.Point-Objekt wird der Standort des Pixels definiert.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
In Ihrer Anwendung wählen Sie möglicherweise ein anderes Pixel aus. Ersetzen Sie dazu die oben angegebenen Längen- und Breitengradkoordinaten. Alternativ können Sie die spektrale Historie einer Gruppe von Pixeln mit anderen ee.Geometry-Objektdefinitionen wie ee.Geometry.Polygon() und ee.Geometry.Rectangle() zusammenfassen.
Weitere Informationen finden Sie im Abschnitt „Geometries“ des Entwicklerleitfadens.
Landsat-Sensorsammlungen abrufen
Landsat-Sammlungen mit Oberflächenreflexionsvermögen der USGS für OLI, ETM+ und TM abrufen.
Über die Links erhalten Sie weitere Informationen zu den einzelnen Datasets.
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');
Filter für eine Bildsammlung definieren
Definieren Sie einen Filter, um die Bildsammlungen nach der räumlichen Grenze des betreffenden Gebiets, der Hauptsaison für die Fotosynthese und der Qualität einzuschränken.
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)));
Ändern Sie diese Filterkriterien nach Bedarf. Auf den oben verlinkten Seiten mit Datenbeschreibungen finden Sie auf dem Tab „Bildeigenschaften“ weitere Eigenschaften, nach denen Sie filtern können.
Sammlungen vorbereiten
Führen Sie die Sammlungen zusammen, wenden Sie den Filter an und ordnen Sie die prepImg-Funktion allen Bildern zu. Das Ergebnis des folgenden Snippets ist ein einzelnes ee.ImageCollection, das Bilder aus OLI-, ETM+- und TM-Sensorsammlungen enthält, die die Filterkriterien erfüllen und zu analyseready NBR verarbeitet werden.
// 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);
Zeitreihendiagramm mit allen Beobachtungen erstellen
Die Harmonie zwischen den Sensoren kann qualitativ bewertet werden, indem alle Beobachtungen als Streudiagramm dargestellt werden, in dem der Sensor durch die Farbe unterschieden wird. Die ui.Chart.feature.groups-Funktion von Earth Engine bietet dieses Dienstprogramm. Zuerst muss jedoch jedem Bild der beobachtete NBR-Wert als Eigenschaft hinzugefügt werden. Wenden Sie eine Funktion zum Reduzieren von Regionen auf die Bildsammlung an, um den Median aller Pixel zu berechnen, die sich mit dem AOI überschneiden, und legen Sie das Ergebnis als Bildeigenschaft fest.
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'));
});
Wenn Ihr AOI in Ihrer eigenen Anwendung ein großes Gebiet ist, sollten Sie die scale erhöhen und die Argumente bestEffort, maxPixels, tileScale für die Funktion reduceRegion angeben, um sicherzustellen, dass der Vorgang die Beschränkungen für maximale Pixel, Speicher oder Zeitüberschreitung nicht überschreitet. Sie können den Median-Reducer auch durch eine andere Statistik ersetzen. Weitere Informationen finden Sie im Abschnitt Statistiken einer Bildregion des Entwicklerleitfadens.
Die Sammlung kann jetzt von der Funktion ui.Chart.feature.groups akzeptiert werden.
Die Funktion erwartet eine Sammlung und die Namen der Eigenschaften von Sammlungsobjekten, die dem Diagramm zugeordnet werden sollen. Hier dient die Eigenschaft „system:time_start“ als Variable für die X-Achse und „NBR“ als Variable für die Y-Achse. Die NBR-Eigenschaft wurde von der reduceRegion-Funktion oben basierend auf dem Bandnamen im Bild benannt. Wenn Sie ein anderes Band als „NBR“ verwenden, ändern Sie das Argument für den Namen der Y-Achsen-Eigenschaft entsprechend. Legen Sie schließlich die Gruppierungsvariable (Reihe) auf „SATELLITE“ fest, der unterschiedliche Farben zugewiesen sind.
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);
Nach einer gewissen Verarbeitungszeit wird in der Konsole ein Diagramm ähnlich Abbildung 2 angezeigt. Beachten Sie Folgendes:
- Es gibt mehrere Beobachtungen pro Jahr.
- Die Zeitreihe besteht aus drei verschiedenen Sensoren.
- Es gibt keine erkennbaren Sensorabweichungen.
- Die Häufigkeit der Beobachtungen variiert jährlich.
- Es gibt einige jährliche Abweichungen, die über die Jahre hinweg konsistent und relativ gering sind. Wenn Sie dies mit NDVI versuchen (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');), sehen Sie eine viel größere intra- und interannuelle Reaktionsvariabilität.
- Der NBR bleibt für die meisten Zeitreihen hoch, mit einigen geringfügigen Störungen.
- Eine deutliche, schnelle Reduzierung der NBR-Reaktion ist auf einen Waldbrand zurückzuführen, der in
- Das Feuer (Dollar Lake) ereignete sich jedoch im September 2011. Der Fehler im Jahr der Erkennung ist darauf zurückzuführen, dass der jährliche zusammengesetzte Zeitraum von Juli bis August ist. Änderungen, die nach diesem Zeitraum auftreten, werden erst bei der nächsten verfügbaren Beobachtung erfasst, die nicht null ist. Das kann im nächsten Jahr oder später sein.
- Die Erholung der Vegetation (zunehmende NBR-Reaktion) beginnt zwei Jahre nach dem Erkennen des großen NBR-Verlusts.
Abbildung 2. Zeitachsendiagramm für die spektrale Reaktion für ein einzelnes Pixel mit Daten von Landsat TM, ETM+ und OLI. TM- und ETM+-Bilder werden durch lineare Transformation an OLI angepasst. Die Bilder stammen aus dem Zeitraum von Juli bis August und wurden nach hoher Qualität gefiltert.
Zeitreihendiagramm mit jährlichem Median erstellen
Um die Zeitreihe zu vereinfachen und Rauschen zu entfernen, kann die Anzahl der Beobachtungen innerhalb eines Jahres reduziert werden. Hier wird der Median verwendet.
Im ersten Schritt werden die Bilder nach Jahr gruppiert. Fügen Sie jedem Bild ein neues Attribut „year“ hinzu, indem Sie das Attribut „year“ aus den ee.Image.Date jedes Bildes zuordnen.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
Mit dem neuen Attribut „year“ (Jahr) können Bilder aus demselben Jahr über einen Join gruppiert werden. Durch die Verknüpfung einer Sammlung mit unterschiedlichen Bildjahren und der vollständigen Sammlung werden Listen für dasselbe Jahr erstellt, aus denen Medianreduktionen vorgenommen werden können. Teilen Sie die vollständige Sammlung in eine Reihe von repräsentativen Jahren auf.
var distinctYearCol = col.distinct('year');
Definieren Sie einen Filter und einen Join, mit denen übereinstimmende Jahre zwischen der Sammlung mit eindeutigen Jahren (distinctYearCol) und der vollständigen Sammlung (col) ermittelt werden. Übereinstimmungen werden in einer Eigenschaft namens „year_matches“ gespeichert.
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Führen Sie den Join aus und konvertieren Sie die resultierende FeatureCollection in eine 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}));
});
Ändern Sie die ee.Reducer oben, wenn Sie für die Darstellung der Bildersammlung in einem bestimmten Jahr eine andere Statistik als den Median verwenden möchten.
Erstellen Sie schließlich ein Diagramm, in dem die medianen NBR-Werte aus Gruppen von intra-annuellen Sommerbeobachtungen dargestellt werden. Ändern Sie die Statistik zur Regionsreduzierung nach Bedarf.
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);
Nach einiger Zeit wird in der Konsole ein Diagramm ähnlich Abbildung 3 angezeigt.
Abbildung 3. Diagramm mit der mittleren spektralen Reaktionszeit im Sommer für einen einzelnen Pixel mit Daten von Landsat TM, ETM+ und OLI. TM- und ETM+-Bilder werden durch eine lineare Transformation an OLI angepasst. Die Bilder stammen aus dem Zeitraum von Juli bis August und wurden nach hoher Qualität gefiltert.
Weitere Informationen
Alternative Transformationsfunktionen
In Roy et al. (2016), Tabelle 2, finden Sie OLS- und RMA-Regressionskoeffizienten zum Transformieren der Oberflächenreflexion von ETM+ in die Oberflächenreflexion von OLI und umgekehrt. Im obigen Tutorial wird nur die Transformation von ETM+ zu OLI durch OLS-Regression demonstriert. Funktionen für alle Übersetzungsoptionen finden Sie im Code Editor-Script oder im GitHub-Quellscript.
Der Dollar Lake-Brand
Das in dieser Anleitung erwähnte Feuer ist das Dollar Lake-Feuer, das Mitte September 2011 an den Nordhängen des Mount Hood in Oregon, USA, ausbrach. Weitere Informationen finden Sie in diesem Blogpost und bei Varner et al. (2015).