AVISO: os procedimentos descritos neste tutorial para harmonizar dados do Landsat ETM+ e do OLI estão desatualizados e não são recomendados nem necessários ao trabalhar com dados de refletância da superfície da Coleção 2 do Landsat. Para mais informações, consulte a entrada de perguntas frequentes sobre a harmonização do Landsat. Os dados da Coleção 1 foram descontinuados, e o código deste tutorial não vai mais funcionar nem ser atualizado.
Este tutorial trata da harmonização da refletância da superfície do Landsat ETM+ com a refletância da superfície do Landsat OLI. Ele oferece:
- Uma função de transformação espectral
- Funções para criar dados prontos para análise
- Exemplo de visualização de série temporal
Ele foi criado para ser um guia completo sobre como harmonizar e visualizar uma série temporal regional de mais de 35 anos de dados do Landsat, que pode ser aplicada imediatamente às suas regiões de interesse.
Os coeficientes para transformar ETM+ em OLI também se aplicam ao TM. Portanto, neste tutorial, a referência a ETM+ é sinônimo de TM.
Sobre o Landsat
O Landsat é um programa de imagens de satélite que coleta imagens da Terra de resolução moderada desde 1972. Como o programa de observação da Terra mais longo com base no espaço, ele fornece um registro temporal valioso para identificar tendências espaço-temporais na mudança da paisagem. Os dados do instrumento Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) e Operational Land Imager (OLI) são usados neste tutorial. Eles são intimamente relacionados e relativamente fáceis de consolidar em uma série temporal consistente que produz um registro contínuo de 1984 até o presente, em uma cadência de 16 dias por sensor com resolução espacial de 30 metros. O instrumento Multispectral Scanner estende o registro do Landsat até 1972, mas não é usado neste tutorial. Os dados são bem diferentes, o que dificulta a integração com sensores posteriores. Consulte Savage et al. (2018) e Vogeler et al. (2018) para exemplos de harmonização em todos os sensores.
Por que a harmonização
Roy et al. (2016) demonstram que há diferenças pequenas, mas potencialmente significativas, entre as características espectrais do Landsat ETM+ e do OLI, dependendo da aplicação. Motivos para harmonizar conjuntos de dados: produzir uma série temporal longa que abrange Landsat TM, ETM+ e OLI, gerar composições intraanuais quase na data para reduzir os efeitos de observações ausentes de lacunas ETM+ SLC-off e mascaramento de nuvem/sombra ou aumentar a frequência de observação em uma série temporal. Consulte o manuscrito vinculado acima para mais informações.
Instruções
Funções
A seguir, há uma série de funções necessárias para harmonizar o ETM+ com o OLI e gerar dados prontos para análise que serão usados na seção de aplicação deste tutorial para visualizar o histórico espectro-temporal de um pixel.
Harmonização
A harmonização é feita por transformação linear do espaço espectral ETM+ para o espaço espectral OLI de acordo com os coeficientes apresentados na Tabela 2 de Roy et al. (2016), coeficientes de regressão OLS.
Os coeficientes de banda são definidos no seguinte dicionário com constantes de imagem de inclinação (slopes) e intercepto (itcps). Os valores de intercepto y são multiplicados por 10.000 para corresponder ao escalonamento dos dados de refletância da superfície do 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])
};
Os nomes de banda do ETM+ e do OLI para a mesma janela de resposta espectral não são iguais e precisam ser padronizados. As funções a seguir selecionam apenas bandas de refletância e a banda pixel_qa de cada conjunto de dados, renomeando-as de acordo com o intervalo de comprimento de onda que representam.
// 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 fim, defina a função de transformação, que aplica o modelo linear aos dados de ETM+, converte o tipo de dados como Int16 para consistência com o OLI e reconecta a banda pixel_qa para uso posterior na mascaramento de nuvem e sombra.
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'));
}
Mascaramento de nuvens e sombras
Os dados prontos para análise precisam ter nuvens e sombras de nuvens mascaradas. A função a seguir usa a CFmask (Zhu et al., 2015)
A banda pixel_qa é incluída em cada imagem de refletância da superfície do Landsat USGS para definir como nulos os pixels identificados como nuvem e sombra de nuvem.
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 do índice espectral
O aplicativo usa o índice espectral de razão de queimada normalizada (NBR) para representar o histórico espectral de um pixel florestado afetado por incêndio florestal. O NBR é usado porque Cohen et al. (2018) descobriram que, entre 13 índices/bandas espectrais, o NBR tinha a maior relação sinal-ruído em relação à perturbação florestal (sinal) em todo os EUA. Ele é calculado pela seguinte função.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
No seu aplicativo, você pode usar um índice espectral diferente. Aqui é onde você altera ou adiciona um índice.
Combinar funções
Defina uma função de wrapper para cada conjunto de dados que consolide todas as funções acima para facilitar a aplicação às respectivas coleções de imagens.
// 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()));
}
No seu aplicativo, você pode incluir ou excluir funções. Altere essas funções conforme necessário.
Exemplo de série temporal
As funções de wrapper prepOli e prepEtm acima podem ser mapeadas em coleções de refletância da superfície do Landsat para criar dados prontos para análise entre sensores e visualizar a cronologia espectral de um pixel ou uma região de pixels. Neste exemplo, você vai criar uma série temporal de mais de 35 anos e mostrar o histórico espectral de um único pixel. Esse pixel específico relaciona o histórico recente de um trecho de floresta de coníferas madura do noroeste do Pacífico (Figura 1) que passou por alguma perturbação nas décadas de 1980 e 1990 e um incêndio de grande magnitude em 2011.
Figura 1. Localização e caractere do site para exemplo de área de interesse. Floresta de coníferas madura do noroeste do Pacífico na encosta norte do Monte Hood, Oregon, EUA. Imagens cortesia de: Google Earth, USDA Forest Service, Landsat e Copernicus.
Definir uma área de interesse (AOI)
O resultado dessa aplicação é uma série temporal de observações do Landsat para um pixel. Um objeto ee.Geometry.Point é usado para definir o local do pixel.
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
No seu aplicativo, você pode selecionar um pixel diferente. Para isso, substitua as coordenadas de longitude e latitude acima. Como alternativa, você pode resumir o histórico espectral de um grupo de pixels usando outras definições de objetos ee.Geometry, como ee.Geometry.Polygon() e ee.Geometry.Rectangle().
Consulte a seção Geometrias do Guia para desenvolvedores para mais informações.
Recuperar coleções de sensores do Landsat
Receba coleções de refletância de superfície do Landsat USGS para OLI, ETM+ e TM.
Acesse os links para saber mais sobre cada conjunto de dados.
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');
Definir um filtro de coleção de imagens
Defina um filtro para restringir as coleções de imagens pela fronteira espacial da área de interesse, pela temporada de pico de fotossíntese e pela qualidade.
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)));
Altere esses critérios de filtragem como quiser. Identifique outras propriedades para filtrar na guia "Propriedades das imagens" das páginas de descrição de dados vinculadas acima.
Preparar as coleções
Mescle as coleções, aplique o filtro e mapeie a função prepImg em todas as imagens. O resultado do snippet a seguir é um único ee.ImageCollection
que contém imagens das coleções de sensores OLI, ETM+ e TM que atendem aos critérios
de filtro e são processadas para NBR prontas para análise.
// 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);
Crie um gráfico de série temporal mostrando todas as observações
A harmonia entre sensores pode ser avaliada qualitativamente representando todas as observações como um diagrama de dispersão em que o sensor é discernido por cor. A função ui.Chart.feature.groups do Earth Engine fornece essa utilidade, mas primeiro o valor de NBR observado precisa ser adicionado a cada imagem como uma propriedade. Mapeie uma função de redução de região na coleção de imagens que calcula a mediana de todos os pixels que se cruzam com a região de interesse e define o resultado como uma propriedade de imagem.
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'));
});
No seu aplicativo, se a AOI for uma área grande, considere aumentar o scale e especificar os argumentos bestEffort, maxPixels e tileScale para a função reduceRegion. Assim, você garante que a operação não exceda os limites máximos de pixels, memória ou tempo limite. Você também pode substituir o redutor de mediana pela estatística de sua preferência. Consulte a seção Estatísticas de uma região de imagem do guia para desenvolvedores e saiba mais.
A coleção agora pode ser aceita pela função ui.Chart.feature.groups.
A função espera uma coleção e os nomes das propriedades do objeto de coleção
para mapear no gráfico. Aqui, a propriedade "system:time_start" serve como a variável do eixo x e "NBR" como a variável do eixo y. A propriedade NBR foi nomeada pela função reduceRegion acima com base no nome da banda na imagem sendo reduzido. Se você usar uma banda diferente (não "NBR"), mude o argumento do nome da propriedade do eixo y de acordo. Por fim, defina a variável de agrupamento (série) como "SATELLITE", a que cores distintas são atribuídas.
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);
Um gráfico semelhante à Figura 2 vai aparecer no console após algum tempo de processamento. Leve isto em conta:
- Há várias observações por ano.
- A série temporal é composta por três sensores diferentes.
- Não há vieses de sensor discerníveis.
- A frequência de observação varia anualmente.
- Há uma variância intra-anual consistente em todos os anos e relativamente pequena. Ao tentar isso com NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) mostra uma variabilidade de resposta intra e interanual muito maior.
- O NBR permanece alto na maioria das séries temporais, com algumas pequenas perturbações.
- Uma redução grande e rápida na resposta de NBR resulta de um incêndio florestal evidente em
- No entanto, o incêndio (Dollar Lake) ocorreu em setembro de 2011. O erro no ano de detecção ocorre porque o período anual composto é de julho a agosto. As mudanças que ocorrem após esse período não são detectadas até a próxima observação não nula disponível, que pode ser no ano seguinte ou depois.
- A recuperação da vegetação (aumento da resposta do NBR) começa dois anos após a detecção da perda significativa do NBR.
Figura 2. Gráfico de série temporal de resposta espectral para um único pixel com representação do Landsat TM, ETM+ e OLI. As imagens TM e ETM+ são harmonizadas com o OLI por transformação linear. As imagens são de julho a agosto e foram filtradas para alta qualidade.
Criar um gráfico de série temporal mostrando a mediana anual
Para simplificar e remover ruídos da série temporal, é possível aplicar a redução de observações intraanuais. Aqui, a mediana é usada.
A primeira etapa é agrupar as imagens por ano. Adicione uma nova propriedade "year" a cada imagem mapeando a configuração de coleção "year" do ee.Image.Date de cada imagem.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
A nova propriedade "year" pode ser usada para agrupar imagens do mesmo ano por uma junção. A junção entre uma coleção de anos de imagens distintas e a coleção completa fornece listas do mesmo ano, em que as reduções de mediana podem ser realizadas. Crie um subconjunto da coleção completa para um conjunto de representantes distintos de ano.
var distinctYearCol = col.distinct('year');
Defina um filtro e uma junção que identifiquem anos correspondentes entre a coleção de anos distintos (distinctYearCol) e a coleção completa (col). As correspondências serão salvas em uma propriedade chamada "year_matches".
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Aplique a junção e converta o FeatureCollection resultante em um 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}));
});
Mude o ee.Reducer acima se preferir uma estatística diferente da mediana para representar a coleção de imagens em um determinado ano.
Por fim, crie um gráfico que mostre os valores medianos de NBR de conjuntos de observações de verão intraanuais. Mude a estatística de redução de região como quiser.
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);
Um gráfico semelhante à Figura 3 vai aparecer no console após algum tempo de processamento.
Figura 3. Gráfico de série temporal do tempo de resposta espectral médio de verão para um único pixel com representação do Landsat TM, ETM+ e OLI. As imagens do TM e do ETM+ são harmonizadas com o OLI por uma transformação linear. As imagens são de julho a agosto e são filtradas para alta qualidade.
Informações adicionais
Funções de transformação alternativas
A Tabela 2 de Roy et al. (2016) fornece coeficientes de regressão de OLS e RMA para transformar a refletância da superfície ETM+ em OLI e vice-versa. O tutorial acima demonstra apenas a transformação de ETM+ para OLI por regressão OLS. As funções para todas as opções de tradução podem ser encontradas no script do Editor de código ou no script de origem do GitHub.
O incêndio em Dollar Lake
O incêndio mencionado neste tutorial é o Dollar Lake, que ocorreu em meados de setembro de 2011 nas encostas do norte do Monte Hood, em Oregon, EUA. Acesse esta postagem do blog e consulte Varner et al. (2015) para mais informações sobre o assunto.