Время прочтения: 4 минуты
Версия QGIS: 3.40.1
Уровень владения QGIS: необходимо уметь ориентироваться в основных кнопках калькулятора растров
Уровень владения QGIS: необходимо уметь ориентироваться в основных кнопках калькулятора растров
Определение температуры поверхности земли (Land Surface Temperature) — классическая и при этом трендовая задача с использованием геоинформационных систем. Уже не надо быть ученым, считающим доли градусов, чтобы заметить, что климат стал гораздо теплее. Жители центральной части России начинают узнавать, что такое сплит-система, а картографы и не только — все чаще интересоваться, насколько жарко стало в их среде обитания и где же находятся те самые острова тепла.
В этой статье вы узнаете один из классических способов определения температуры поверхности земли с помощью снимков Landsat, а также о подводных камнях, с которыми могли бы столкнуться, если бы не прочитали статью.
Исходные данные
Чтобы определить температуру поверхности земли, нам потребуются снимки Landsat 8 или 9. Скачать их можно со старого доброго EarthExplorer. Для нашей работы потребуются снимки Landsat Collection 2 Level-1. Если вы скачаете Level-2, то результаты будут просто огонь в прямом смысле этого слова, но это не наш вариант, нам нужно получить корректные вычисления.
Для работы нам потребуются 4, 5 и 10 каналы снимка, а также файл метаданных. Названия файлов будут выглядеть примерно так:
- LC09_L1TP_172026_20240712_20240712_02_T1_B4
- LC09_L1TP_172026_20240712_20240712_02_T1_B5
- LC09_L1TP_172026_20240712_20240712_02_T1_B10
- LC09_L1TP_172026_20240712_20240712_02_T1_MTL
Скачивая снимки, старайтесь выбирать те, на которых как можно меньше облаков. В идеале нужно, чтобы их вообще не было, но тут уж как вам повезет.
Что нас ожидает
Для того, чтобы рассчитать температуру поверхности земли, нам потребуется провести ряд вычислений, чтобы перейти от сырых данных 10 канала снимка к яркостной температуре, а потом скорректировать ее с учетом излучательной способности (emissivity). Мы сделаем это с помощью формул, которые будем вбивать в калькулятор растров.
Рассчитываем температуру
Расчет спектральной энергетической яркости
1. Первым шагом мы рассчитываем спектральную энергетическую яркость (the top atmospheric spectral radiance):
Lλ = ML × Qcal + AL, где
ML — «RADIANCE_MULT_BAND_10» из файла метаданных, того, что заканчивается на MTL;
Qcal — значение пикселя 10 канала снимка;
AL — «RADIANCE_ADD_BAND_10» из файла метаданных.
Яркостная температура
2. Теперь с помощью Lλ вычислим яркостную температуру (brightness temperature, BT) в градусах Цельсия:
BT = (K2 / ln( (K1/Lλ) + 1) ) - 273.15, где
K1 — «K1_CONSTANT_BAND_10» из файла метаданных;
K2 — «K2_CONSTANT_BAND_10» из файла метаданных.
Температуры на этом этапе уже должны выглядеть правдоподобно времени году и климату исследуемой местности. Если вы анализируете территорию в Мурманске осенью, а у вас получаются температуры под 60 градусов, то формула сработала не так, либо в этом месте во время пролета спутника было что-то интересное.
Время, в которое спутник сделал снимок территории, можно узнать тоже в файле метаданных, оно скрывается на строчке «SCENE_CENTER_TIME»
Коррекция на излучательную способность
3. Теперь проведем коррекцию на излучательную способность поверхности. Для этого нам требуется сделать ряд вычислений.
3.1. Сначала рассчитываем нормализованный индекс вегетации (NDVI):
(NIR - RED) / (NIR + RED), где
NIR — значение пикселя 5 канала;
RED — значение пикселя 4 канала
3.2. Теперь рассчитываем долю растительного покрова:
Pv =( (NDVI - NDVIs) / (NDVIv - NDVIs) )^2, где
NDVI — слой с NDVI, который мы рассчитали в шаге 3.1
NDVIs = 0.2
NDVIv = 0.5
3.3. Предпоследний шаг, рассчитываем излучательную способность земной поверхности ε:
ε = 0.004 * Pv + 0.986
3.4. Финальный шаг, как финальный босс с самой длинной формулой, рассчитываем температуру поверхности земли (LST):
LST = BT / (1 + (λ * BT / ρ) × ln(ε)), где помимо слоев, которые мы использовали:
λ — длина волны излученной радиации, среднее значение которой принято 10.895;
ρ — постоянная, равная 14388 мкК
λ — длина волны излученной радиации, среднее значение которой принято 10.895;
ρ — постоянная, равная 14388 мкК
Интерпретация и что еще почитать
Температура поверхности земли — это один из индикаторов городских островов тепла, а также с помощью нее определяют интенсивность поверхностного острова тепла. Благодаря спутниковым снимкам это очень доступный для изучения показатель, но только на основании него нельзя точно делать выводы, что в каком-то месте жарче, чем в другом, ведь приповерхностный слой воздуха может рассеиваться ветром или смешиваться с другими потоками прохладного воздуха, которые могут быть вызваны как локальными явлениями, так и явлениями большего масштаба.
Если вы хотите узнать больше об этой теме, то могу посоветовать вам для начала эти ресурсы:
- Плейлист с лекциями ученых из МГУ с курса «Климатология с основами метеорологии» от Кислова А.В. и Константинова П.И., там есть и про городские острова тепла и про климатообразование, как в городе, так и в целом. Ютуб, но есть и на rutube, и ВК Видео.
- Книга «Urban Climates»” от T.R. Oke, G. Mills, A. Christen and J.A. Voogt
- «Algorithm for Automated Mapping of Land Surface Temperature Using LANDSAT 8 Satellite Data» Ugur Avdan and Gordana Jovanovska — статья с методикой определения LST, представленной в статье.
Материал подготовил Александр Зуев