AVVISO: le procedure descritte in questo tutorial per armonizzare i dati Landsat ETM+ e OLI sono obsolete e non consigliate o necessarie quando si lavora con i dati di riflettanza di superficie di Landsat Collection 2. Per ulteriori informazioni, consulta la voce delle domande frequenti sull'armonizzazione di Landsat. I dati della raccolta 1 sono ritirati e il codice in questo tutorial non funzionerà più né verrà aggiornato.
Questo tutorial riguarda l'armonizzazione della riflettanza di superficie di Landsat ETM+ con la riflettanza di superficie di Landsat OLI. Nella guida troverai:
- Una funzione di trasformazione spettrale
- Funzioni per creare dati pronti per l'analisi
- Esempio di visualizzazione di una serie temporale
È pensata per essere una guida end-to-end per armonizzare e visualizzare una serie temporale regionale di dati Landsat di oltre 35 anni, che può essere applicata immediatamente alle tue regioni di interesse.
Tieni presente che i coefficienti per la trasformazione di ETM+ in OLI si applicano anche a TM. Pertanto, in questo tutorial, il riferimento a ETM+ è sinonimo di TM.
Informazioni su Landsat
Landsat è un programma di imaging satellitare che raccoglie immagini della Terra a risoluzione moderata dal 1972. In quanto programma di osservazione della Terra dallo spazio più lungo, fornisce un prezioso record temporale per identificare le tendenze spaziotemporali nel cambiamento del paesaggio. In questo tutorial vengono utilizzati i dati degli strumenti Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) e Operational Land Imager (OLI). Sono strettamente correlati e relativamente facili da consolidare in una serie temporale coerente che produce una registrazione continua dal 1984 a oggi, con una cadenza di 16 giorni per sensore con una risoluzione spaziale di 30 metri. Lo strumento Multispectral Scanner estende il record Landsat fino al 1972, ma non viene utilizzato in questo tutorial. I suoi dati sono molto diversi, il che rende difficile l'integrazione con i sensori successivi. Consulta Savage et al. (2018) e Vogeler et al. (2018) per esempi di armonizzazione tra tutti i sensori.
Perché l'armonizzazione
Roy et al. (2016) dimostrano che esistono piccole, ma potenzialmente significative differenze tra le caratteristiche spettrali di Landsat ETM+ e OLI, a seconda dell'applicazione. I motivi per cui potresti voler armonizzare i set di dati includono: produrre una serie temporale lunga che comprenda Landsat TM, ETM+ e OLI, generare compositi intra-annuali quasi contemporanei per ridurre gli effetti delle osservazioni mancanti dovuti alle lacune SLC-off di ETM+ e alla maschera di nuvole/ombre o aumentare la frequenza delle osservazioni all'interno di una serie temporale. Per ulteriori informazioni, consulta il manoscritto collegato sopra.
Istruzioni
Funzioni
Di seguito è riportata una serie di funzioni necessarie per armonizzare ETM+ con OLI e generare dati pronti per l'analisi che verranno utilizzati nella sezione dell'applicazione di questo tutorial per visualizzare la cronologia spettrale-temporale di un pixel.
Armonizzazione
L'armonizzazione viene ottenuta tramite la trasformazione lineare dello spazio spettrale ETM+ in spazio spettrale OLI in base ai coefficienti presentati nella
tabella 2 dei coefficienti di regressione OLS di Roy et al. (2016).
I coefficienti specifici per banda sono definiti nel seguente dizionario con le costanti immagine pendenza (slopes) e intercetta (itcps). Tieni presente che i valori dell'intercetta
vengono moltiplicati per 10.000 per corrispondere alla scalabilità dei dati di riflettanza
di superficie 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])
};
I nomi delle bande ETM+ e OLI per la stessa finestra di risposta spettrale non sono uguali e
devono essere standardizzati. Le seguenti funzioni selezionano solo le bande di riflettanza
e la banda pixel_qa di ogni set di dati e le rinominano in base all'intervallo di
lunghezza d'onda che rappresentano.
// 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']);
}
Infine, definisci la funzione di trasformazione, che applica il modello lineare ai dati ETM+, esegue il cast del tipo di dati come Int16 per coerenza con OLI e ricollega la banda pixel_qa per un utilizzo successivo nella maschera di nuvole e ombre.
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'));
}
Mascheramento di nuvole e ombre
I dati pronti per l'analisi devono avere le nuvole e le ombre delle nuvole mascherate. La
funzione seguente utilizza CFmask (Zhu et al., 2015)
pixel_qa inclusa in ogni immagine di riflettanza di superficie USGS di Landsat per impostare
i pixel identificati come nuvole e ombre delle nuvole su 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);
}
Calcolo dell'indice spettrale
L'applicazione in arrivo utilizza l'indice spettrale del rapporto di combustione normalizzato (NBR) per rappresentare la storia spettrale di un pixel forestale interessato da un incendio boschivo. L'NBR viene utilizzato perché Cohen et al. (2018) hanno scoperto che, tra 13 indici/bande spettrali, l'NBR aveva il rapporto segnale/rumore più elevato per quanto riguarda il disturbo forestale (segnale) in tutti gli Stati Uniti. Viene calcolato mediante la seguente funzione.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
Nella tua applicazione, puoi scegliere di utilizzare un indice spettrale diverso. Qui puoi modificare o aggiungere un indice aggiuntivo.
Funzioni combinate
Definisci una funzione wrapper per ogni set di dati che consolidi tutte le funzioni precedenti per facilitarne l'applicazione alle rispettive raccolte di immagini.
// 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()));
}
Nella tua applicazione, puoi decidere di includere o escludere le funzioni. Modifica queste funzioni in base alle esigenze.
Esempio di serie temporale
Le funzioni wrapper prepOli e prepEtm riportate sopra possono essere mappate sulle raccolte di riflettanza di superficie Landsat per creare dati pronti per l'analisi tra sensori per visualizzare la cronologia spettrale di un pixel o di una regione di pixel. In questo
esempio, creerai una serie temporale di oltre 35 anni e visualizzerai la cronologia spettrale
per un singolo pixel. Questo pixel specifico si riferisce alla storia recente di una zona di foresta di conifere matura del Pacifico nord-occidentale (Figura 1) che ha subito alcune perturbazioni negli anni'80 e'90 e un incendio di grande intensità nel 2011.
Figura 1. Posizione e carattere del sito per l'area di interesse di esempio. Foresta di conifere matura del Pacifico nord-occidentale sul versante nord del Monte Hood, Oregon, Stati Uniti. Immagini per gentile concessione di: Google Earth, USDA Forest Service, Landsat e Copernicus.
Definire un'area di interesse
Il risultato di questa applicazione è una serie temporale di osservazioni Landsat per un
pixel. Un oggetto ee.Geometry.Point viene utilizzato per definire la posizione del pixel.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
Nella tua applicazione, potresti selezionare un pixel diverso. Per farlo, sostituisci le coordinate di longitudine e latitudine riportate sopra. In alternativa, puoi riepilogare
la cronologia spettrale di un gruppo di pixel utilizzando altre definizioni
di oggetti ee.Geometry, come ee.Geometry.Polygon() e ee.Geometry.Rectangle().
Per saperne di più, consulta la sezione
Geometrie
della Guida per gli sviluppatori.
Recupera le raccolte di sensori Landsat
Ottieni le raccolte di riflettanza di superficie USGS di Landsat per OLI, ETM+ e TM.
Visita i link per scoprire di più su ogni set di dati.
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');
Definisci un filtro per la raccolta di immagini
Definisci un filtro per vincolare le raccolte di immagini in base al confine spaziale dell'area di interesse, alla stagione di massima fotosintesi e alla 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)));
Modifica questi criteri di filtro come preferisci. Identifica altre proprietà da filtrare nella scheda "Proprietà immagini" delle pagine di descrizione dei dati collegate sopra.
Preparare le raccolte
Unisci le raccolte, applica il filtro e mappa la funzione prepImg su tutte le
immagini. Il risultato del seguente snippet è un singolo ee.ImageCollection
che contiene immagini delle raccolte di sensori OLI, ETM+ e TM che soddisfano i
criteri di filtro e vengono elaborate in NBR pronte per l'analisi.
// 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 grafico delle serie temporali che mostri tutte le osservazioni
L'armonia tra i sensori può essere valutata qualitativamente tracciando tutte le osservazioni
come un grafico a dispersione in cui il sensore è distinto per colore. La funzione ui.Chart.feature.groups di Earth Engine fornisce questa utilità, ma prima il valore NBR osservato deve essere aggiunto a ogni immagine come proprietà. Mappa una funzione di riduzione
della regione sulla raccolta di immagini che calcola la mediana di tutti i pixel
che intersecano l'AOI e imposta il risultato come proprietà dell'immagine.
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'));
});
Nella tua applicazione, se l'AOI è un'area di grandi dimensioni, potresti prendere in considerazione
l'aumento di scale e la specifica degli argomenti bestEffort, maxPixels, tileScale
per la funzione reduceRegion per assicurarti che l'operazione non
superi i limiti massimi di pixel, memoria o timeout. Puoi anche sostituire il
riduttore della mediana con la statistica che preferisci. Per saperne di più, consulta la sezione
Statistiche di una regione di un'immagine
della Guida per gli sviluppatori.
La raccolta ora può essere accettata dalla funzione ui.Chart.feature.groups.
La funzione prevede una raccolta e i nomi delle proprietà degli oggetti della raccolta
da mappare al grafico. In questo caso, la proprietà "system:time_start" funge da variabile dell'asse x e "NBR" da variabile dell'asse y. Tieni presente che la proprietà NBR è stata
denominata dalla funzione reduceRegion sopra in base al nome della banda nell'immagine
che viene ridotta. Se utilizzi una banda diversa (non "NBR"), modifica l'argomento del nome della proprietà dell'asse y di conseguenza. Infine, imposta la variabile di raggruppamento (serie) su "SATELLITE", a cui sono assegnati colori distinti.
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);
Dopo un po' di tempo di elaborazione, nella console verrà visualizzato un grafico simile a quello della Figura 2. Tieni presente diverse cose:
- Ci sono più osservazioni all'anno.
- La serie temporale è composta da tre sensori diversi.
- Non sono presenti distorsioni dei sensori rilevabili.
- La frequenza delle osservazioni varia ogni anno.
- Esiste una certa varianza intra-annuale, che è coerente nel corso degli anni e relativamente ristretta. Se provi a farlo con l'NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) la variabilità della risposta intra- e interannuale è molto maggiore.
- L'NBR rimane elevato per la maggior parte della serie temporale con alcune piccole perturbazioni.
- Una grande e rapida riduzione della risposta NBR è il risultato di un incendio boschivo evidente in
- Tuttavia, tieni presente che l'incendio (Dollar Lake) si è verificato a settembre 2011. L'errore nell'anno di rilevamento è dovuto al fatto che l'intervallo di date composito annuale va da luglio ad agosto; le modifiche che si verificano dopo questo periodo non vengono rilevate fino alla successiva osservazione non nulla disponibile, che potrebbe essere l'anno successivo o successivo.
- Il recupero della vegetazione (aumento della risposta NBR) inizia due anni dopo il rilevamento della perdita NBR principale.
Figura 2. Grafico delle serie temporali del tempo di risposta spettrale per un singolo pixel con rappresentazione di Landsat TM, ETM+ e OLI. Le immagini TM ed ETM+ vengono armonizzate a OLI mediante trasformazione lineare. Le immagini sono di luglio e agosto e sono filtrate per l'alta qualità.
Crea un grafico delle serie temporali che mostri la mediana annuale
Per semplificare e rimuovere il rumore dalla serie temporale, è possibile applicare la riduzione delle osservazioni intra-annuali. In questo caso viene utilizzata la mediana.
Il primo passaggio consiste nel raggruppare le immagini per anno. Aggiungi una nuova proprietà "anno" a ogni
immagine mappando l 'impostazione della raccolta "anno" da ogni immagine
ee.Image.Date.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
La nuova proprietà "anno" può essere utilizzata per raggruppare le immagini dello stesso anno tramite un'unione. L'unione tra una raccolta di anni di immagini distinte e la raccolta completa fornisce elenchi dello stesso anno, dai quali è possibile eseguire le riduzioni mediane. Sottoinsieme della raccolta completa in un insieme di rappresentanti distinti per anno.
var distinctYearCol = col.distinct('year');
Definisci un filtro e un'unione che identifichino gli anni corrispondenti tra la raccolta di anni distinti (distinctYearCol) e la raccolta completa (col). Le corrispondenze verranno salvate in una proprietà denominata "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Applica l'unione e converti la FeatureCollection risultante in 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}));
});
Modifica il valore ee.Reducer sopra se preferisci che una statistica diversa dalla mediana
rappresenti la raccolta di immagini in un determinato anno.
Infine, crea un grafico che mostri i valori NBR mediani di insiemi di osservazioni estive intra-annuali. Modifica la statistica di riduzione della regione come preferisci.
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);
Dopo un po' di tempo di elaborazione, nella console verrà visualizzato un grafico simile alla Figura 3.
Figura 3. Grafico delle serie temporali del tempo di risposta spettrale estivo mediano per un singolo pixel con rappresentazione di Landsat TM, ETM+ e OLI. Le immagini TM ed ETM+ vengono armonizzate a OLI mediante una trasformazione lineare. Le immagini sono state scattate tra luglio e agosto e sono filtrate per garantire l'alta qualità.
Informazioni aggiuntive
Funzioni di trasformazione alternative
La tabella 2 di Roy et al. (2016) fornisce i coefficienti di regressione OLS e RMA per trasformare la riflettanza di superficie ETM+ in riflettanza di superficie OLI e viceversa. Il tutorial precedente mostra solo la trasformazione da ETM+ a OLI mediante la regressione OLS. Le funzioni per tutte le opzioni di traduzione sono disponibili nello script dell'editor di codice o nello script di origine di GitHub.
L'incendio di Dollar Lake
L'incendio menzionato in questo tutorial è l'incendio di Dollar Lake, che si è verificato a metà settembre 2011 sui versanti settentrionali del Monte Hood, in Oregon, Stati Uniti. Visita questo post del blog e consulta Varner et al. (2015) per saperne di più.