ВНИМАНИЕ: Процедуры, описанные в этом руководстве по гармонизации данных Landsat ETM+ и OLI, устарели и не рекомендуются или не требуются при работе с данными об отражательной способности поверхности Landsat Collection 2. Подробнее см. в разделе часто задаваемых вопросов о гармонизации Landsat Collection 1. Данные Collection 1 устарели, и код в этом руководстве больше не будет работать или обновляться.
В этом руководстве рассматривается согласование отражательной способности поверхности Landsat ETM+ с отражательной способностью поверхности Landsat OLI. В нем представлены:
- Функция спектрального преобразования
- Функции для создания данных, готовых к анализу
- Пример визуализации временного ряда
Он призван стать комплексным руководством по гармонизации и визуализации более чем 35-летнего регионального временного ряда данных Landsat, который можно немедленно применить к интересующему вас региону (регионам).
Обратите внимание, что коэффициенты преобразования ETM+ в OLI применимы и к TM. Поэтому в рамках данного руководства ETM+ является синонимом TM.
О компании Landsat
Landsat — это программа спутниковой съемки, которая собирает изображения Земли среднего разрешения с 1972 года. Будучи самой продолжительной программой космических наблюдений Земли, она предоставляет ценную временную запись для выявления пространственно-временных тенденций в изменении ландшафта. В этом руководстве используются данные инструментов Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+) и Operational Land Imager (OLI). Они тесно связаны и относительно легко объединяются в последовательный временной ряд, который создает непрерывную запись с 1984 года по настоящее время, с частотой 16 дней на датчик с пространственным разрешением 30 метров. Инструмент Multispectral Scanner продлевает запись Landsat до 1972 года, но не используется в этом руководстве. Его данные существенно отличаются, что затрудняет интеграцию с более поздними датчиками. См. Savage et al. (2018) и Vogeler et al. (2018) для примеров гармонизации по всем датчикам.
Зачем нужна гармонизация
Рой и др. (2016) демонстрируют наличие небольших, но потенциально значимых различий между спектральными характеристиками Landsat ETM+ и OLI, в зависимости от области применения. Гармонизация наборов данных может быть целесообразна при: создании длительных временных рядов, охватывающих Landsat TM, ETM+ и OLI; формировании внутригодовых композитных данных с близкими датами для снижения влияния пропусков наблюдений из-за пропусков SLC-off ETM+ и маскирования облачностью/тенью; или увеличении частоты наблюдений в пределах временного ряда. Подробнее см. в рукописи, ссылка на которую приведена выше.
Инструкции
Функции
Ниже приведен ряд функций, необходимых для согласования ETM+ с OLI и создания готовых к анализу данных, которые будут использоваться в разделе «Применение» этого руководства для визуализации спектрально-временной истории пикселя.
Гармонизация
Гармонизация достигается посредством линейного преобразования спектрального пространства ETM+ в спектральное пространство OLI в соответствии с коэффициентами, представленными в работе Роя и др. (2016). Таблица 2. Коэффициенты регрессии OLS. Коэффициенты, соответствующие полосам, определены в следующем словаре с константами изображения: уклон ( slopes ) и пересечение ( itcps ). Обратите внимание, что значения пересечения по оси Y умножаются на 10 000 для соответствия масштабированию данных об отражательной способности поверхности, полученных со спутника 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])
};
Названия полос ETM+ и OLI для одного и того же окна спектрального отклика не совпадают и должны быть стандартизированы. Следующие функции выбирают только полосы отражения и полосу pixel_qa из каждого набора данных и переименовывают их в соответствии с диапазоном длин волн, который они представляют.
// 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']);
}
Наконец, определите функцию преобразования, которая применяет линейную модель к данным ETM+, приводит тип данных к Int16 для согласованности с OLI и повторно присоединяет полосу pixel_qa для последующего использования при маскировании облаков и теней.
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'));
}
Маскировка облаков и теней
В готовых к анализу данных облака и тени от них должны быть замаскированы. Следующая функция использует канал pixel_qa CFmask ( Zhu et al., 2015 ), включенный в каждый снимок отражения поверхности Landsat USGS, чтобы обнулить пикселы, идентифицированные как облака и тени от них.
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);
}
Расчет спектрального индекса
В готовящемся приложении для представления спектральной истории пикселя, затронутого лесным пожаром, используется спектральный индекс нормализованного коэффициента выгорания (NBR). NBR используется, поскольку Коэн и др. (2018) обнаружили, что среди 13 спектральных индексов/диапазонов NBR имеет наибольшее отношение сигнал/шум в отношении лесных нарушений (сигнала) по всей территории США. Он рассчитывается с помощью следующей функции.
function calcNbr(img) {
return img.normalizedDifference(['NIR', 'SWIR2']).rename('NBR');
}
В вашем приложении вы можете использовать другой спектральный индекс. Здесь вы можете изменить или добавить дополнительный индекс.
Объединить функции
Определите функцию-оболочку для каждого набора данных, которая объединяет все вышеперечисленные функции для удобства их применения к соответствующим коллекциям изображений.
// 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()));
}
Вы можете включить или исключить функции из своего приложения. Измените эти функции по мере необходимости.
Пример временного ряда
Описанные выше функции-оболочки prepOli и prepEtm можно применять к наборам данных отражения поверхности Landsat для создания данных, готовых к кросс-сенсорному анализу, и визуализации спектральной хронологии пикселя или области пикселей. В этом примере мы создадим временной ряд продолжительностью более 35 лет и отобразим спектральную историю для одного пикселя. Этот конкретный пиксель описывает недавнюю историю участка зрелого хвойного леса на северо-западе Тихого океана (рис. 1), который испытал некоторые нарушения в 1980-х и 90-х годах, а также сильный пожар в 2011 году.
Рисунок 1. Местоположение и характеристики участка для примера исследуемой территории. Зрелый хвойный лес Тихоокеанского Северо-Запада на северном склоне горы Худ, штат Орегон, США. Изображения предоставлены: Google Earth, Лесной службой Министерства сельского хозяйства США, Landsat и Copernicus.
Определить область интереса (ОИ)
Результатом работы этого приложения является временной ряд наблюдений Landsat для пикселя. Для определения местоположения пикселя используется объект ee.Geometry.Point .
var aoi = ee.Geometry.Point([-121.70938, 45.43185]);
В вашем приложении вы можете выбрать другой пиксель. Для этого замените указанные выше координаты долготы и широты. Кроме того, вы можете суммировать спектральную историю группы пикселей, используя другие определения объектов ee.Geometry , такие как ee.Geometry.Polygon() и ee.Geometry.Rectangle() . Подробнее см. в разделе «Геометрия» Руководства разработчика.
Получить коллекции датчиков Landsat
Получите коллекции данных об отражательной способности поверхности Landsat USGS для OLI , ETM+ и TM .
Чтобы узнать больше о каждом наборе данных, перейдите по ссылкам.
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');
Определить фильтр коллекции изображений
Определите фильтр, ограничивающий коллекции изображений по пространственной границе интересующей области, пиковому сезону фотосинтеза и качеству.
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)));
Измените эти критерии фильтрации по своему усмотрению. Другие свойства для фильтрации можно определить на вкладке «Свойства изображений» на страницах описания данных, ссылки на которые приведены выше.
Подготовьте коллекции
Объедините коллекции, примените фильтр и примените функцию prepImg ко всем изображениям. Результатом следующего фрагмента кода является единый ee.ImageCollection , содержащий изображения из коллекций датчиков OLI, ETM+ и TM, соответствующие критериям фильтра и обработанные для анализа в формате NBR.
// 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);
Создайте диаграмму временного ряда, отображающую все наблюдения.
Межсенсорную гармонию можно оценить качественно, построив диаграмму рассеяния для всех наблюдений, где датчики различаются по цвету. Функция ui.Chart.feature.groups в Earth Engine предоставляет такую возможность, но сначала необходимо добавить наблюдаемое значение NBR к каждому изображению в качестве свойства. Примените функцию сокращения области к коллекции изображений, которая вычисляет медиану всех пикселей, пересекающих область интереса, и задаёт результат как свойство изображения.
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'));
});
Если в вашем приложении область интереса (AOI) представляет собой большую область, вы можете увеличить scale и указать аргументы bestEffort , maxPixels , tileScale для функции reduceRegion , чтобы гарантировать, что операция не превысит ограничения по максимальному количеству пикселей, памяти или времени ожидания. Вы также можете заменить медианный редуктор на предпочитаемую вами статистику. Подробнее см. в разделе «Статистика области изображения» Руководства разработчика.
Теперь коллекция может быть принята функцией ui.Chart.feature.groups . Функция ожидает, что коллекция и имена свойств объектов коллекции будут сопоставлены с диаграммой. Здесь свойство 'system:time_start' служит переменной оси X, а 'NBR' — переменной оси Y. Обратите внимание, что свойство NBR было названо функцией reduceRegion выше на основе имени полосы на уменьшаемом изображении. Если вы используете другую полосу (не 'NBR'), измените аргумент имени свойства оси Y соответствующим образом. Наконец, установите для переменной группировки (серии) значение 'SATELLITE', которому будут назначены отдельные цвета.
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);
После некоторой обработки в консоли появится диаграмма, похожая на рисунок 2. Обратите внимание на несколько моментов:
- Проводится несколько наблюдений в год.
- Временной ряд состоит из трех различных датчиков.
- Не существует никаких заметных смещений датчиков.
- Частота наблюдений меняется ежегодно.
- Существует некоторая внутригодовая дисперсия, которая сохраняется по годам и относительно невелика. Проверка этого на NDVI (img.normalizedDifference(['SWIR1', 'NIR']).rename('NDVI');) показывает гораздо большую внутри- и межгодовую изменчивость отклика.
- NBR остается высоким для большей части временного ряда с некоторыми незначительными отклонениями.
- Значительное, быстрое снижение реакции NBR происходит в результате лесного пожара, очевидного в
- Однако следует отметить, что пожар ( Доллар Лейк ) произошел в сентябре 2011 года. Ошибка в годе обнаружения обусловлена тем, что годовой составной диапазон дат составляет период с июля по август; изменения, происходящие после этого периода, не учитываются до следующего доступного ненулевого наблюдения, которое может быть в следующем году или позже.
- Восстановление растительности (усиление реакции NBR) начинается через два года после обнаружения значительной потери NBR.
Рисунок 2. Временной ряд спектрального отклика для одного пикселя с данными Landsat TM, ETM+ и OLI. Изображения TM и ETM+ гармонизированы с OLI посредством линейного преобразования. Изображения получены с июля по август и отфильтрованы для повышения качества.
Создайте диаграмму временного ряда, отображающую годовую медиану
Для упрощения и устранения шума во временном ряду можно применить сокращение внутригодовых наблюдений. В данном случае используется медиана.
Первый шаг — сгруппировать изображения по году. Добавьте новое свойство «year» к каждому изображению, сопоставив его с параметром коллекции «year» из ee.Image.Date каждого изображения.
var col = col.map(function(img) {
return img.set('year', img.date().get('year'));
});
Новое свойство «year» можно использовать для группировки изображений одного года с помощью объединения. Объединение коллекции изображений одного года с полной коллекцией создаёт списки изображений одного года, на основе которых можно выполнить медианное сокращение. Разделите полную коллекцию на подмножество, состоящее из отдельных репрезентативных изображений разных лет.
var distinctYearCol = col.distinct('year');
Определите фильтр и соединение, которые будут определять совпадающие годы между коллекцией отдельных лет ( distinctYearCol ) и полной коллекцией ( col ). Совпадения будут сохраняться в свойстве year_matches.
var filter = ee.Filter.equals({leftField: 'year', rightField: 'year'});
var join = ee.Join.saveAll('year_matches');
Примените объединение и преобразуйте полученную FeatureCollection в 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}));
});
Измените ee.Reducer выше, если вы предпочитаете, чтобы для представления коллекции изображений за определенный год использовалась статистика, отличная от медианы.
Наконец, создайте диаграмму, отображающую медианные значения NBR на основе наборов внутригодовых летних наблюдений. Изменяйте статистику сокращения региона по своему усмотрению.
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);
Через некоторое время обработки в консоли появится диаграмма, аналогичная рисунку 3.
Рисунок 3. Диаграмма временного ряда медианного летнего спектрального отклика для одного пикселя с данными Landsat TM, ETM+ и OLI. Изображения TM и ETM+ гармонизированы с OLI посредством линейного преобразования. Изображения получены с июля по август и отфильтрованы для повышения качества.
Дополнительная информация
Альтернативные функции преобразования
Рой и др. (2016) В таблице 2 приведены коэффициенты регрессии МНК и РМА для преобразования отражательной способности поверхности ETM+ в отражательную способность поверхности OLI и наоборот. В приведенном выше руководстве показано только преобразование ETM+ в OLI с помощью регрессии МНК. Функции для всех вариантов преобразования можно найти в скрипте редактора кода или исходном скрипте на GitHub .
Пожар в Доллар-Лейк
Пожар, упомянутый в этом руководстве, — это пожар на озере Доллар, произошедший в середине сентября 2011 года на северных склонах горы Худ, штат Орегон, США. Подробнее об этом можно узнать в этой публикации блога Varner et al. (2015) .