OSTRZEŻENIE: procedury opisane w tym samouczku w celu harmonizacji danych Landsat ETM+ i OLI są przestarzałe i nie są zalecane ani konieczne podczas pracy z danymi odbicia powierzchniowego Landsat Collection 2. Więcej informacji znajdziesz w odpowiedzi na najczęstsze pytania dotyczące harmonizacji danych Landsat. Dane z kolekcji 1 są przestarzałe, a kod w tym samouczku nie będzie już działać ani nie będzie aktualizowany.
Ten samouczek dotyczy harmonizacji współczynnika odbicia światła od powierzchni w przypadku Landsat ETM+ ze współczynnikiem odbicia światła od powierzchni w przypadku Landsat OLI. Oto co w nim znajdziesz:
- funkcja przekształcenia widmowego,
- Funkcje tworzenia danych gotowych do analizy
- Przykład wizualizacji ciągu czasowego
Jest to kompleksowy przewodnik po harmonizacji i wizualizacji regionalnych ciągów czasowych danych Landsat z ponad 35-letniego okresu, które można od razu zastosować w interesujących Cię regionach.
Pamiętaj, że współczynniki przekształcania danych z ETM+ na OLI mają zastosowanie również do TM. Dlatego w tym samouczku odniesienia do ETM+ są synonimem TM.
Informacje o satelicie Landsat
Landsat to program obrazowania satelitarnego, który od 1972 roku gromadzi zdjęcia Ziemi w średniej rozdzielczości. Jest to najdłuższy program obserwacji Ziemi z kosmosu, który zapewnia cenne dane czasowe do identyfikowania trendów czasowo-przestrzennych w zakresie zmian krajobrazu. W tym samouczku używane są dane z instrumentów Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) i Operational Land Imager (OLI). Są one ze sobą ściśle powiązane i stosunkowo łatwo je połączyć w spójny szereg czasowy, który zapewnia ciągłe dane od 1984 r. do chwili obecnej w odstępach 16-dniowych na czujnik i rozdzielczości przestrzennej 30 m. Skaner wielospektralny rozszerza dane Landsat do 1972 roku, ale nie jest używany w tym samouczku. Jego dane są dość odmienne, co utrudnia integrację z późniejszymi czujnikami. Przykłady harmonizacji danych ze wszystkich czujników znajdziesz w artykułach Savage i in. (2018) oraz Vogeler i in. (2018).
Dlaczego harmonizacja
Roy i współautorzy (2016) wykazali, że w zależności od zastosowania występują niewielkie, ale potencjalnie istotne różnice między charakterystyką spektralną Landsat ETM+ i OLI. Powody, dla których warto harmonizować zbiory danych, to m.in.: tworzenie długich ciągów czasowych obejmujących Landsat TM, ETM+ i OLI, generowanie kompozytów wewnątrzrocznych z datami bliskimi bieżącej, aby zmniejszyć wpływ brakujących obserwacji z luk ETM+ SLC-off i maskowania chmur/cieni, lub zwiększenie częstotliwości obserwacji w ciągu czasowym. Więcej informacji znajdziesz w rękopisie, do którego link znajduje się powyżej.
Instrukcje
Funkcje
Poniżej znajdziesz serię funkcji, które są potrzebne do ujednolicenia danych ETM+ i OLI oraz wygenerowania danych gotowych do analizy, które zostaną użyte w sekcji dotyczącej aplikacji w tym samouczku do wizualizacji historii spektralno-czasowej piksela.
Harmonizacja
Harmonizację uzyskuje się poprzez liniową transformację przestrzeni spektralnej ETM+ do przestrzeni spektralnej OLI zgodnie ze współczynnikami przedstawionymi w tabeli 2 współczynników regresji OLS w artykule Roy i wsp. (2016).
Współczynniki dla poszczególnych pasm są zdefiniowane w tym słowniku ze stałymi obrazu nachylenia (slopes) i punktu przecięcia (itcps). Pamiętaj, że wartości punktu przecięcia z osią Y są mnożone przez 10 000, aby dopasować je do skali danych o odbiciu powierzchniowym 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])
};
Nazwy pasm ETM+ i OLI dla tego samego okna odpowiedzi spektralnej nie są równe i wymagają ujednolicenia. Poniższe funkcje wybierają z każdego zbioru danych tylko pasma odbicia i pasmo pixel_qa, a następnie zmieniają ich nazwy zgodnie z zakresem długości fali, który reprezentują.
// 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']);
}
Na koniec zdefiniuj funkcję przekształcenia, która stosuje model liniowy do danych ETM+, przekształca typ danych na Int16, aby zachować spójność z OLI, i ponownie dołącza pasmo pixel_qa do późniejszego użycia w maskowaniu chmur i cieni.
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'));
}
Maskowanie chmur i cieni
Dane gotowe do analizy powinny mieć zamaskowane chmury i ich cienie. Poniższa funkcja korzysta z maski CF (Zhu i in., 2015)
pixel_qa dołączony do każdego obrazu odbicia powierzchniowego Landsat USGS, aby ustawić piksele zidentyfikowane jako chmury i cienie chmur na wartość 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);
}
Obliczanie indeksu spektralnego
Nadchodząca aplikacja wykorzystuje znormalizowany wskaźnik wypalenia (NBR) do przedstawiania historii spektralnej piksela lasu dotkniętego pożarem. Wskaźnik NBR jest używany, ponieważ Cohen i współautorzy (2018) stwierdzili, że spośród 13 indeksów i pasm spektralnych wskaźnik NBR miał największy stosunek sygnału do szumu w odniesieniu do zakłóceń w lasach (sygnał) w całych Stanach Zjednoczonych. Oblicza się go za pomocą tej funkcji.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
W swojej aplikacji możesz użyć innego indeksu spektralnego. W tym miejscu możesz zmienić lub dodać kolejny indeks.
Łączenie funkcji
Zdefiniuj funkcję opakowującą dla każdego zbioru danych, która będzie zawierać wszystkie powyższe funkcje, aby ułatwić stosowanie ich do odpowiednich kolekcji obrazów.
// 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()));
}
W aplikacji możesz uwzględniać lub wykluczać funkcje. W razie potrzeby zmień te funkcje.
Przykład serii czasowej
Funkcje opakowujące prepOli i prepEtm powyżej można mapować na kolekcje odbicia powierzchni Landsat, aby tworzyć dane gotowe do analizy między czujnikami i wizualizować chronologię spektralną piksela lub regionu pikseli. W tym przykładzie utworzysz szereg czasowy obejmujący ponad 35 lat i wyświetlisz historię spektralną dla pojedynczego piksela. Ten konkretny piksel odnosi się do niedawnej historii dojrzałego fragmentu lasu iglastego w północno-zachodniej części Pacyfiku (rysunek 1), który w latach 80. i 90. XX wieku uległ pewnym zaburzeniom, a w 2011 r. został dotknięty pożarem o dużej intensywności.
Rysunek 1. Lokalizacja i charakter witryny w przypadku przykładowego obszaru zainteresowań. Dojrzały las iglasty na północnym zboczu góry Hood w Oregonie, USA. Zdjęcia pochodzą z Google Earth, USDA Forest Service, Landsat i Copernicus.
Określanie obszaru zainteresowań
Wynikiem działania tej aplikacji jest szereg czasowy obserwacji Landsat dla piksela. Obiekt ee.Geometry.Point służy do określania lokalizacji piksela.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
W aplikacji możesz wybrać inny piksel. W tym celu zastąp powyższe współrzędne szerokości i długości geograficznej. Możesz też podsumować historię spektralną grupy pikseli, korzystając z innych definicji ee.Geometryobiektówee.Geometry, takich jak ee.Geometry.Polygon() i ee.Geometry.Rectangle().
Więcej informacji znajdziesz w sekcji Geometrie w przewodniku dla programistów.
Pobieranie kolekcji danych z czujników Landsat
Uzyskaj kolekcje odbicia powierzchni USGS Landsat dla OLI, ETM+ i TM.
Kliknij linki, aby dowiedzieć się więcej o poszczególnych zbiorach danych.
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');
Określanie filtra kolekcji obrazów
Określ filtr, aby ograniczyć kolekcje obrazów według granicy przestrzennej obszaru zainteresowania, szczytu sezonu fotosyntezy i jakości.
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)));
Zmień te kryteria filtrowania w dowolny sposób. Na karcie „Właściwości obrazów” na stronach opisu danych, do których linki znajdziesz powyżej, możesz określić inne właściwości, według których chcesz filtrować dane.
Przygotowywanie kolekcji
Połącz kolekcje, zastosuj filtr i przypisz funkcję prepImg do wszystkich obrazów. Wynikiem tego fragmentu kodu jest pojedynczy element ee.ImageCollection
zawierający obrazy z kolekcji czujników OLI, ETM+ i TM, które spełniają kryteria filtra i są przetwarzane do postaci NBR gotowej do analizy.
// 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);
Utwórz wykres z osią czasu, na którym będą widoczne wszystkie obserwacje.
Spójność między czujnikami można ocenić jakościowo, przedstawiając wszystkie obserwacje na wykresie punktowym, na którym czujnik jest oznaczony kolorem. Funkcja ui.Chart.feature.groups w Earth Engine zapewnia tę funkcjonalność, ale najpierw do każdego obrazu należy dodać zaobserwowaną wartość NBR jako właściwość. Zastosuj do kolekcji obrazów funkcję mapRegion reduce, która oblicza medianę wszystkich pikseli przecinających obszar zainteresowania i ustawia wynik jako właściwość obrazu.
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'));
});
Jeśli w Twojej aplikacji obszar zainteresowania jest duży, możesz zwiększyć wartość parametru scale i określić argumenty bestEffort, maxPixels, tileScale
funkcji reduceRegion, aby mieć pewność, że operacja nie przekroczy maksymalnej liczby pikseli, pamięci ani limitów czasu. Możesz też zastąpić funkcję mediany preferowaną statystyką. Więcej informacji znajdziesz w sekcji Statystyki regionu obrazu w Przewodniku dla programistów.
Kolekcja może teraz zostać zaakceptowana przez funkcję ui.Chart.feature.groups.
Funkcja oczekuje kolekcji i nazw właściwości obiektu kolekcji, które mają być mapowane na wykres. W tym przypadku właściwość „system:time_start” służy jako zmienna osi X, a „NBR” jako zmienna osi Y. Zwróć uwagę, że właściwość NBR została nazwana przez funkcję reduceRegion powyżej na podstawie nazwy pasma na obrazie po jej skróceniu. Jeśli używasz innego pasma (nie „NBR”), odpowiednio zmień argument nazwy właściwości osi Y. Na koniec ustaw zmienną grupowania (serii) na „SATELLITE”, do której przypisane są różne kolory.
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);
Po pewnym czasie przetwarzania w konsoli pojawi się wykres podobny do tego na rysunku 2. Zwróć uwagę na kilka kwestii:
- W ciągu roku jest wiele obserwacji.
- Ciąg czasowy składa się z 3 różnych czujników.
- Nie ma zauważalnych odchyleń czujnika.
- Częstotliwość obserwacji zmienia się co roku.
- Występują pewne wahania w ciągu roku, które są spójne w poszczególnych latach i stosunkowo niewielkie. Próba z użyciem NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) wykazuje znacznie większą zmienność wewnątrz- i międzyroczną.
- W przypadku większości szeregów czasowych wskaźnik NBR pozostaje na wysokim poziomie, z niewielkimi odchyleniami.
- Znaczny i szybki spadek wartości NBR w wyniku pożaru lasu widoczny na
- Pamiętaj jednak, że pożar (Dollar Lake) wybuchł we wrześniu 2011 r. Błąd w roku wykrycia wynika z tego, że roczny złożony zakres dat obejmuje okres od lipca do sierpnia. Zmiany zachodzące po tym okresie nie są wykrywane aż do następnej dostępnej obserwacji niezerowej, która może nastąpić w następnym roku lub później.
- Regeneracja roślinności (wzrost wartości NBR) rozpoczyna się 2 lata po wykryciu znacznej utraty wartości NBR.
Rysunek 2. Wykres szeregów czasowych odpowiedzi spektralnej dla pojedynczego piksela z danymi z Landsat TM, ETM+ i OLI. Obrazy TM i ETM+ są harmonizowane z obrazami OLI za pomocą transformacji liniowej. Obrazy pochodzą z okresu od lipca do sierpnia i są filtrowane pod kątem wysokiej jakości.
Utwórz wykres z osią czasu przedstawiający roczną medianę.
Aby uprościć szereg czasowy i usunąć z niego szum, można zastosować redukcję obserwacji w ciągu roku. W tym przypadku używana jest mediana.
Najpierw pogrupuj zdjęcia według roku. Dodaj do każdego obrazu nową właściwość „year”, mapując ustawienie kolekcji „year” z ee.Image.Date każdego obrazu.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
Nowa właściwość „rok” może służyć do grupowania obrazów z tego samego roku za pomocą złączenia. Połączenie kolekcji obrazów z różnych lat z pełną kolekcją zapewnia listy z tego samego roku, z których można obliczyć medianę redukcji. Wyodrębnij z pełnej kolekcji podzbiór zawierający reprezentatywne dane z poszczególnych lat.
var distinctYearCol = col.distinct('year');
Zdefiniuj filtr i złączenie, które będą identyfikować pasujące lata między kolekcją unikalnych lat (distinctYearCol) a pełną kolekcją (col). Pasujące wartości zostaną zapisane we właściwości o nazwie „year_matches”.
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Zastosuj złączenie i przekonwertuj wynikową kolekcję obiektów na kolekcję obrazów.
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}));
});
Jeśli wolisz, aby zbiór obrazów z danego roku był reprezentowany przez inną statystykę niż mediana, zmień wartość powyżej ee.Reducer.
Na koniec utwórz wykres przedstawiający mediany wartości NBR z zestawów obserwacji letnich w ciągu roku. Zmień statystykę zmniejszenia regionu w dowolny sposób.
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);
Po pewnym czasie przetwarzania w konsoli pojawi się wykres podobny do tego na rysunku 3.
Rysunek 3. Wykres szeregów czasowych przedstawiający medianę czasu reakcji spektralnej w okresie letnim dla pojedynczego piksela z danych z satelitów Landsat TM, ETM+ i OLI. Obrazy TM i ETM+ są harmonizowane z obrazami OLI za pomocą transformacji liniowej. Zdjęcia pochodzą z okresu od lipca do sierpnia i są filtrowane pod kątem wysokiej jakości.
Dodatkowe informacje
Alternatywne funkcje przekształcenia
W tabeli 2 w artykule Roy i wsp. (2016) podano współczynniki regresji OLS i RMA, które umożliwiają przekształcenie współczynnika odbicia powierzchniowego ETM+ na współczynnik odbicia powierzchniowego OLI i odwrotnie. Powyższy samouczek pokazuje tylko przekształcenie ETM+ na OLI za pomocą regresji OLS. Funkcje wszystkich opcji tłumaczenia znajdziesz w skrypcie edytora kodu lub skrypcie źródłowym GitHub.
Pożar Dollar Lake
Pożar, o którym mowa w tym samouczku, to pożar Dollar Lake, który wybuchł w połowie września 2011 r. na północnych zboczach góry Mt Hood w stanie Oregon w USA. Więcej informacji znajdziesz w tym poście na blogu oraz w artykule Varner i in. (2015).