Оценка температуры поверхности из снимка Landsat-8 при помощи Land Surface Temperature QGIS Plugin

Эта страница является черновиком статьи.

В статье описан алгоритм оценки температуры земной поверхности (LST) по снимку Landsat-8 полуавтоматическим способом при помощи плагина QGIS Land Surface Temperature. А также в целом рассмотрен процесс оценки температуры земной поверхности по данным Дистанционного зондирования Земли (ДЗЗ). По важным вопросам приводится спектр мнений специалистов, из которых делаются выводы.

1.Какие данные ДЗЗ о температуре земной поверхности существуют в открытом доступе?

Температуру земной поверхности можно оценить по термальным спектральным каналам сенсора, то есть таким каналам, которые снимают земную поверхность в диапазоне Thermal Infrared Radiation (10-15 микрометров). Соответственно, эти каналы называются по-английски Thermal bands, TIRS Bands или Thermal Infrared bands.
Подробнее можно узнать в статье ДЗЗ для экологических задач. Часть 1: Введение в теорию ДЗЗ (рис. 3)

Термальные каналы имеются лишь у некоторых сенсоров. Поэтому немного и данных ДЗЗ, по которым можно оценить температуру земной поверхности за наше время и за прошлый период. В целом, существующие данные можно разделить на два класса масштаба (в скобках указано число метров поверхности в 1 пикселе снимка термального канала и годы получения данных):

Данные ДЗЗ температуры земной поверхности мелкого масштаба

• MODIS (500 и 1000 м, 2000 – н.в.), сенсоры спутников Terra, Aqua

Температуру земной поверхности можно получить в готовом виде через продукты Thermal Anomalies and Fire. Подробнее: MODIS Products Table

• Suomi NPP (650 м, 2011 – н.в.)

Подробнее: Suomi NPP: краткая характеристика

Данные ДЗЗ температуры земной поверхности среднего масштаба

• ASTER (90 м, 2000 – н.в.), сенсор спутника Terra

Температуру земной поверхности можно получить в готовом виде через продукт AST_08 (ASTER L2 Surface Temperature V003). Подробнее: AST_08: ASTER Surface Kinetic Temperature V003

• Данные спутников Landsat:

Landsat-8 (100 м, 2013 – н.в.)
Landsat-7 (60 м, 1999- н.в.)
Landsat-5 (120 м, 1984-2013)
Landsat-4 (120 м, 1982-1993)

Для научных и практических исследований регионального масштаба наибольшую ценность представляют данные температуры земной поверхности среднего масштаба. На время написания статьи, Landsat-8 дает наилучшие такие данные. В целом данные о температуре земной поверхности можно получить с 1982 года (год запуска Landsat-4).

В последние годы сенсор ASTER по неясным причинам значительно (вероятно, во много раз) сократил объем снимков. На многие или на все участки земной поверхности уже нельзя получить данные за любой интересующий период, наблюдаются провалы в несколько лет без снимков вообще. До 2013 года такого за данными ASTER не замечалось.

Снимки Landsat-7 с 31 мая 2003 г идут с дефектом, который не отражается на точности данных, но значительно сокращает их объем и снижает удобство их использования (усложняет обработку). Однако, разрешение термального канала Landsat-7 (60 м) при этом остается самым высоким из когда-либо существовавших в открытом доступе.

В статье мы рассмотрим расчет данных температуры земной поверхности по снимку Landsat-8 при использовании Land Surface Temperature QGIS plugin. Этот же плагин и изложенные аспекты оценки температуры земной поверхности можно применять и при расчете температуры по Landsat-4, Landsat-5 и Landsat-7 (в последнем не забывая о техническом дефекте, решение которого данная статья не рассматривает).

2.Факторы, которые необходимо учитывать при оценке температуры земной поверхности по снимкам Landsat

При оценке температуры земной поверхности по данным Landsat нужно учитывать два фактора, которые определяют разнообразие методов расчета температуры и, соответственно, различающийся итоговый результат и точность оценки. Это:

• Коррекция снимков Landsat
• Учет разного характера излучательной способности земной поверхности (emissivity)

Эти факторы мы рассмотрим в ходе соответствующих шагов по построению карты температуры.

Из-за наличия этих факторов не существует единственного алгоритма расчета температуры земной поверхности по снимкам Landsat, на который просто можно дать ссылку.
В принципе, температуру земной поверхности по Landsat-8 (и Landsat-4,5,7) можно рассчитать в любом растровом калькуляторе. Именно это предлагают делать многочисленные руководства (видео-мануалы, блоги) легко находимые в Google по поисковым словам. Однако, если вы хотите воспользоваться таким руководством, разберитесь в сути формул, которые они используют, поскольку разные расчеты дают разные результаты и точность оценки температуры. Иначе говоря, разные пособия предлагают разные методы расчета температуры по снимкам Landsat. Пользователям при этом не сообщается о заложенном в данном пособии принципе и о том, что есть другие варианты расчетов.

3. Алгоритм вычисления температуры земной поверхности из снимка Landsat-8 при помощи Land Surface Temperature QGIS Plugin

Шаг 1. Скачиваем снимок

Скачиваем снимок Landsat-8 на интересующую территорию и дату. Желательно выбирать снимки с отсутствием облаков или с небольшой (до 3-5%) облачностью, или же следить (через функции предпросмотра), чтобы облака не покрывали территорию интереса. Ссылки на некоторые популярные ресурсы выбора и скачивания снимков Landsat-8 представлены ниже:

• EarthExplorer https://earthexplorer.usgs.gov
• GloVis https://glovis.usgs.gov/app?fullscreen=1
• Landsat on AWS https://landsatonaws.com
• Данные Landsat-8, Sentinel-2, CBERS-4 (удобный предпросмотр) https://search.remotepixel.ca/#3/40/-70.5
• Libra: https://libra.developmentseed.org

Мы будем строить карту по данным снимка LC08_L1TP_220076_20180501_20180501_01_RT, который можно скачать или по данной ссылке, или, например, в EarthExplorer, пройдя по Data Sets > Landsat > Landsat Collection 1 Level-1

Информация о дате и точном времени съемки можно найти в текстовом файле MLT, который приложен к общему набору скачиваемых данных.
Открываем его в Блокноте и видим, что снимок сделан 1 мая 2018 года в 13:09:42.

DATE_ACQUIRED = 2018-05-01, SCENE_CENTER_TIME = "13:09:42.0926400Z"

Время съемки указано в GMT (Greenwich Mean Time), о чем также говорит буква Z ("Zulu time", аналог GMT). Подробнее – здесь.

Шаг 2. Узнаем про коррекции снимков Landsat

Рассмотрим первый "усложняющий жизнь" фактор – коррекции снимков Landsat.

По ранним данным Landsat мы привыкли, что снимки Landsat поступают к пользователю в сыром виде и их нужно долго и страшно как-то корректировать. Однако с 2016 года всё стало значительно лучше. И теперь все данные Landsat, то есть Landsat-1,4,5,7 и 8 поставляются уже с геометрической и радиометрической коррекцией. Такими к нам поступают снимки, называющиеся в наборах Landsat Level-1 Data Processing Levels или Landsat Level-1 data product или L1TP.

"L1TP" - входит в название скачиваемого продукта, например, как на нашем снимке по которому мы будем строить температуру:

LC08_L1TP_220076_20180501_20180501_01_RT

Если же снимки Landsat имеют в названии L1GT или L1GT, то это означает, что они не выдерживают критерии по пройденным коррекциям, которые имеет основная часть снимков Landsat (L1TP), то есть надо дополнительно изучать, как и для чего можно использовать такие данные. Для неискушенного пользователя можно сделать вывод – их не надо использовать. :)

Подробнее:
Landsat Processing Details
Landsat Collections

Таким образом, пользователь может провести сам еще только одну - атмосферную коррекцию снимков Landsat.

О том, что продукты L1TP не имеют атмосферной коррекции можно заключить из официального ответа на такой вопрос по ссылке:

Does Landsat Level-1 data processing include atmospheric correction?

По этой же ссылке упоминается продукт Landsat Level-2 или Surface Reflectance Product для которого атмосферная коррекция проведена. Данный продукт содержит каналы с 3 по 7, и он не содержит термальных каналов. Поэтому его нельзя использовать для построения температуры земной поверхности.

Вывод: Для построения температуры земной поверхности нам нужен только продукт L1TP. Для более точного построения температуры, этот продукт можно провести через атмосферную коррекцию (как это сделать при помощи плагина QGIS "Land Surface Temperature" рассмотрено ниже).

ПРИМЕЧАНИЕ-1
Дискуссия о необходимости атмосферной коррекции в снимках Landsat
Некоторые специалисты считают, что атмосферную коррекцию проводить бессмысленно, а то и вредно для точности оценки температуры земной поверхности и другого использования Landsat, поскольку еще никто не научился определять точно такую вещь, как, например, количество жидкости в атмосфере (влажность), а этот фактор влияет на данные. Специалисты говорят, что давно ведутся попытки разработки метода оценки показателей атмосферы для коррекции, но пока ученые не достигли в этом особых результатов. Также есть мнение, что той базовой коррекции, которую делает Landsat (L1TP) для всех своих данных достаточно, и что дополнительная коррекция необходима сценам с большим количеством облачности. Соответственно, можно сделать вывод, что, если вы считаете по сценам без облачности, атмосферная коррекция не нужна. Данные мнения собраны при обсуждении вопроса атмосферной коррекции в сообществе GIS-Lab. Пытливый читатель может найти их подтверждения (или опровержения) в литературе и внести свой вклад в точность материала, передав ссылки автору.

Шаг 3. Устанавливаем QGIS плагин "Land Surface Temperature" и разбираемся с формулами

Открываем проект QGIS и загружаем плагин: В верхней панели Модули > Управление модулями > в окне Поиск вбиваем "Land Surface Temperature" >Установить модуль. После установки модуль открывается через верхнюю панель Растр > Land Surface Temperature

Принцип работы плагина заключается в следующем. Двигаясь слева направо по вкладкам: Radiance > Brightness Temperature > NDVI > Land Surface Emissivity > Land Surface Temperature Algorithm мы последовательно выполняем в каждой из них расчет требуемого растра на основе подгружения требуемых каналов снимка Landsat (скачанных и сохраненных на компьютере в отдельной папке "LC08_L1TP_220076_20180501_20180501_01_RT"). На последней вкладке производится расчет самой температуры земной поверхности (Land Surface Temperature, LST) на базе подгружаемых растров, созданных на предыдущих шагах.

Плагин выполняет расчет LST по алгоритму и формулам, с которыми можно ознакомиться, например, здесь:
Estimation of Land Surface Temperature
или здесь:
Tutorial: Estimation of Land Surface Temperature with Landsat and ASTER

Рассмотрим упрощенный алгоритм и формулы расчета температуры земной поверхности по снимку Landsat-8, которые можно выполнить в растровом калькуляторе QGIS или ArcGIS, разобранный на форуме Александром Черепановым:

источник
Здесь:
DN - Исходные значения растра в каналах TIRS.
TOA brightness temperature - Top Of Atmosphere brightness temperature (температура черного тела на поверхности атмосферы).
LST - Land Surface Temperature (температуры земной поверхности).
emissivity - Отражательная способность земной поверхности.
NDVI - Normalized Difference Vegetation Index (вегетационный индекс). Подробнее об NDVI можно узнать здесь: NDVI - теория и практика

• Выберем один канал для расчета температуры (band 10 или band 11 в Landsat-8). Ряд источников рекомендует считать температуру по 10-му каналу. Например: Why using only Landsat 8 band 10 in the estimation of surface temperature?¶

• Мы выбрали 10 канал и для него находим в метаданных параметр RADIANCE_MULT_BAND_10 и RADIANCE_ADD_BAND_10, обычно эти параметры равны 3.3420E-04 и 0.10000 и считаем по формуле:

Lλ=RADIANCE_MULT_BAND_10∗BAND10+RADIANCE_ADD_BAND_10

where:
Lλ-TOA spectral radiance (Watts/(m2*srad*µm))

• Находим в метаданных параметр K1_CONSTANT_BAND_10 и K2_CONSTANT_BAND_10, обычно эти параметры равны 774.8853 и 1321.0789. Теперь мы можем посчитать TOA brightness temperature по формуле:

TB=K2/ln[(K1/Lλ)+1]

where:
TB-TOA brightness temperature (K)
Lλ-TOA spectral radiance (Watts/(m2*srad*µm))
K1-Band-specific thermal conversion constant from the metadata (K1_CONSTANT_BAND_x, where x is the thermal band number)
K2-Band-specific thermal conversion constant from the metadata (K2_CONSTANT_BAND_x, where x is the thermal band number)

• Переход от DN к TOA brightness temperature делается только так( перейти) никаких других вариантов быть не может. Для уменьшения числа операций можно объединить формулы:

TB=K2/ln[(K1/(RADIANCE_MULT_BAND∗BAND10+RADIANCE_ADD_BAND))+1]

• Готовимся к переходу к Land surface temperature и считаем по формуле:

LST=TB/[1+(λ∗TB/c2)∗ln(emissivity)]

where:
λ-wavelength of emitted radiance, для 10 канала - значение 10.8
c2=h∗c/s=1.4388∗10−2 m K = 14388 µm K
h = Planck’s constant = 6.626∗10−34 J s
s = Boltzmann constant = 1.38∗10−23 J/K
c = velocity of light = 2.998∗108 m/s

На этом этапе мы подошли ко второму "усложняющему жизнь фактору" в оценке LST, а именно – к оценке emissivity (отражающей способности земной поверхности).
Для того, чтобы рассчитать LST нам надо определить, что будет использоваться в качестве значения emissivity – классификация типов земной поверхности, значение на основе вегетационного индекса NDVI или константа.
Так, если emissivity ввести в виде константы равной 0.98 и добавить перевод в градусы Цельсия, то формула расчета LST получится такой:

LST=(TB/[1+(10.8∗TB/14388)∗ln(0.98))-273.15
where:
TB-TOA brightness temperature (K)
LST -Land surface temperature (C)

ПРИМЕЧАНИЕ-2
Причина разнообразия методов оценки температуры земной поверхности по данным Landsat
Разнообразие методов оценки температуры земной поверхности обусловлено тем, что значения, которые показывает тепловой канал (после пересчета TOA brightness temperature) и температура земной поверхности (LST) связаны весьма косвенно. Представьте - у нас нет никакой информации о состояние поверхности и атмосферы, есть только излучение зафиксированное съемочной аппаратурой и мы должны определить как это излучение исказилось около земной поверхности и пройдя через атмосферу, прежде чем было зафиксировано аппратурой. Отсюда и множество подходов к оценке LST и подключение дополнительных источников информации о состоянии поверхности и атмосферы. Одни более простые и требуют меньше сторонних данных, другие более сложные (используют NDVI, или\и классификацию типов земной поверхности). Качеством получаемых результатов LST в абсолютных значениях, вероятно, далеки от заявленных 0.1-1 градуса чувствительности аппаратуры.
Именно по этой причине, нельзя использовать для расчета LST "первый попавшийся мануал" не разбираясь в том, какой алгоритм расчета он вам предлагает. Одни пособия считают значение emissivity постоянной величиной, другие используют NDVI, третьи предлагают вам предварительно заняться классификацией земной поверхности.

Шаг 4. Забываем про формулы и просто считаем LST в плагине QGIS "Land Surface Temperature"

Разбор формул был необходим, чтобы понять, что именно делает плагин QGIS "Land Surface Temperature", что предлагают рассчитывать различные пособия и модули других ГИС софтов. Чтобы показать, что разными пособиями, тулбоксами, скриптами дается не один единственный способ вычисления LST из снимков Landsat, а даются разные способы такого вычисления. Разнообразие методов расчета обусловлено двумя факторами – способ учета emissivity и способ учета влияния атмосферы (атмосферная коррекция).

При выборе пособия расчета LST или модуля для полуавтоматического расчета LST из снимков Landsat, следует узнать, какой именно алгоритм расчета LST используется и, если этот метод устраивает для выбранной задачи, то его использовать.
Плагин QGIS "Land Surface Temperature" был позитивно оценен специалистами дистанционного зондирования Земли (в частности, Александром Черепановым).
Видно, что плагин использует при расчете LST индекс NDVI (через него определяя значение emissivity), не использует для значения emissivity классификацию земной поверхности (для оценки растительного покрова NDVI это лучшее, для других сред, возможно, также его достаточно), и также плагин дает возможность внести атмосферные параметры (то есть провести и атмосферную коррекцию), но и имеет способ расчета LST без атмосферной коррекции.

Теперь, когда мы все это поняли и выбрали плагин QGIS "Land Surface Temperature" как подходящий для наших задач оценки температуры земной поверхности из снимков Landsat-8, мы можем забыть о формулах и просто его использовать. Выполняя последовательно всё, что он от нас хочет. Там все очень просто. Приведу принт-скрины моего расчета с краткими пояснениями.

Итак, считаем температуру земной поверхности (LST) по снимку LC08_L1TP_220076_20180501_20180501_01_RT, все данные по ссылке скачаны и записаны в папку на компьютере с таким же названием. С учетом атмосферной коррекции нам нужно сделать 6 последовательных расчетов:

Р1. Расчет Radiance

Р2. Расчет Brightness Temperature

Р3. Расчет NDVI

Для OLI (Landsat-8) каналы для NDVI это: NIR band 5 (0.85 - 0.88 µm, Red band 4 (0.64-0.67 µm).
Источник: каналы Landsat

Р4. Расчет Land Surface Emissivity

Выбрала способ NDVI Threshhold Algorithm

Р5. Расчет атмосферного профиля в Калькуляторе атмосферных параметров

Для того, чтобы получить атмосферные параметры, которые дадут возможность посчитать LST по способу с атмосферной коррекцией "Radiative transfer equation", воспользуемся
Atmospheric Correction Parameter Calculator
Координаты указала той территории, для которой мне важно получиться температуру (а не для центра сцены). Это -22.8194, -47.1105 По карте рельефа и сводке погоды на 1 мая 2018 г в ближайшем городе я внесла в калькулятор эти параметры: Surface altitude 0.7 km, Surface pressure 1019 mb, Surface temperature 29 (C), Surface relative humidity 41%
Получилось так:

Мы получили атмосферный профиль, из которого возьмем нужные параметры для атмосферной коррекции:

Р6. Расчет температуры земной поверхности (LST) с использованием атмосферной коррекции методом Radiative transfer equation

Для того, чтобы построить LST с учетом атмосферной коррекции на последней стадии расчетов выберем метод "Radiative transfer equation". Я выбрала этот метод, поскольку только для него можно получить все необходимые параметры при помощи Калькулятора атмосферных параметров. Другие способы (Mono-Window Algorithm, Single Channel Algorithm)- содержат в себе параметры, которые нужно искать в других источиках. Для того, чтобы построить LST без атмосферной коррекции, на последней стадии выберете метод "Planck Equation".
Вношу нужные параметры из полученного атмосферного профиля:

Upwelling Radiance: 1.7, Downwelling Radiance: 2.81, Atmospheric Transmission: 0.80.
В графу Top of Atmosphere Radiance (TOA) подгружаю растр radiance, рассчитанный в первом расчете плагина.
Выбираю единицы температуры Celsius.


Итоговый растр LST_1may2018atm можно обрабатывать дальше в QGIS, а можно из QGIS пересохранить в Geotiff и открыть в ArcMap. Полученная карта температур может выглядеть таким образом:

Разрешение этой карты – 30 метров земной поверхности в 1 пикселе изображения. Термальный канал Landsat-8 имеет разрешение 100 метров, но карта температур строится при использовании параметров и других каналов, имеющих разрешение 30 метров (для NDVI использовались каналы 4 и 5). Поэтому выходное разрешение карты становится 30 метров.