GEE-10. Библиотека Awesome Spectral Indices

Код примера

В прошлый раз мы вычисляли вегетационные индексы. Есть способ сделать это еще проще — с помощью библиотеки Awesome Spectral Indices (документация).

Создадим коллекцию снимков Sentinel-2 L2A, покрывающих интересующий район (ROI — region of interest):

var col = ee.ImageCollection('COPERNICUS/S2_SR')
            .filterDate('2021-06-01', '2021-09-01')
            .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE',10))
            .filter(ee.Filter.bounds(ROI))
            .select('B.*');

Вычислим медиану коллекции и обрежем ее по границам ROI. Медиана коллекции — снимок, каналы которого являются медианами соответствующих каналов коллекции (вычисления в каждом канале происходят попиксельно). В итоге получим некий типичный летний снимок Небраски.

var image = col.median().clip(ROI);

Загружаем библиотеку Awesome Spectral Indices. Путь к ней указан в параметре функции require.

var spectral = require("users/dmlmont/spectral:spectral");

Awesome Spectral Indices — сторонняя библиотека. Она создана не разработчиками Earth Engine, а пользователем с ником dmlmont. Вы можете создать свою библиотеку и загрузить ее на GEE, но об этом поговорим позже. Сейчас же масштабируем значения пикселей.

var image = spectral.scale(image, "COPERNICUS/S2_SR");

Помните, у нас была функция applyScaleFactors для масштабирования значений пикселей? scale из библиотеки spectral делает то же самое.

Зададим параметры для расчета индексов. Мы хотим рассчитать NDVI, NDWI и EVI. Для этого понадобится несколько каналов и числовых констант.

Список спектральных индексов можно посмотреть здесь. В списке приведена формула, каждой букве которой мы сопоставим канал или число.

var parameters = {
    "N": image.select("B8"),
    "R": image.select("B4"),
    "G": image.select("B3"),
    "B": image.select("B2"),
    "L": 1,
    "g": 2.5,
    "C1": 6,
    "C2": 7.5
};

Вот формула для расчета EVI:

2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))

А вот как она выглядит в списке:

g * ((N - R) / (N + C1 * R - C2 * B + L))

Думаю, принцип понятен. Зеленый канал нужен нам для расчета NDWI.

spectral.computeIndex вычисляет индексы по списку и добавляет к снимку.

var image = spectral.computeIndex(image, ["NDVI","NDWI","EVI"], parameters);

Для GEE Python API есть несколько библиотек, реализующих расчет спектральных индексов.

#GEE #индексы

Vegetation Condition Index

Normalized Difference Vegetation Index (NDVI) является хорошим показателем объема зеленой (фотосинтезирующей) биомассы — чем он больше, тем больше объем зеленой биомассы.* Мы уже обсуждали и использовали этот индекс. Так что здесь просто напомним о нем.

NDVI вычисляется по формуле:

NDVI = (NIR – RED) / (NIR + RED).

Когда солнечный свет попадает на растение, излучение в красной области спектра (0.4-0.7 мкм) поглощается хлорофиллом листьев, тогда как клеточные образования в листьях отражают большую часть излучения в ближней инфракрасной области спектра (NIR) (0.7-1.1 мкм). Таким образом, здоровая растительность поглощает красный свет (RED) и отражает NIR-излучение.

Значения NDVI варьируется в диапазоне от -1 до +1, при этом значения ниже нуля означают отсутствие зеленой растительности. На участках сельскохозяйственных культур NDVI в течение года меняется примерно в диапазоне от 0.3 до 0.9.

С помощью NDVI удобно следить за состоянием растений в течение сезона. Однако, он ничего не скажет о том, нормально ли такое состояние растений в данной местности. Представление о норме будем изменяться в зависимости от климата, почвы, типа растительности и рельефа. Например, значения NDVI, нормальные для высокогорья, будут слишком низкими для тропического леса.

Состояние растений в данной местности можно оценить, сравнивая его с результатами многолетних наблюдений в этой местности. Один их способов такого сравнения дает Vegetation Condition Index (VCI):

VCI_t = (NDVI_t - NDVI_min) / (NDVI_max - NDVI_min ) × 100.

Он сравнивает состояние растительности в данный момент времени (NDVI_t) с максимальным и минимальным значениями NDVI, которые достигались в том же пикселе поверхности и в ту же дату за все предыдущие годы наблюдений.

VCI, как и NDVI, является нормализованным разностным индексом. Обычно он изменяется в диапазоне от 0 (NDVI_t = NDVI_min) до 100 (NDVI_t = NDVI_max), но в рекордные годы может принимать значения за пределами этого диапазона.

Чтобы получить качественные значения VCI нужны многолетние наблюдения. В нашем примере использованы данные MODIS за 18 лет. Феликс Коган, предложивший VCI, использовал данные радиометра AVHRR, временной ряд которых начинается в 1978 году. Не углубляясь в вопросы статистики, примем, что 10 лет достаточно, чтобы оценить норму состояния растительности в данной местности.

Можно считать, что высокие значения VCI соответствуют хорошему состоянию растительности и отсутствию действия на нее стрессовых факторов. Развивая тему стрессовых факторов, можно построить Temperature Condition Index (TCI)

TCI_t = (T_t - TCI_min) / (TCI_max - TCI_min) × 100.

Высокий TCI говорит о том, что на участке необычно жарко. В сочетании с пониженным VCI, это может свидетельствовать о наступлении засухи.

Рассчитывается TCI совершенно аналогично VCI. В качестве исходных данных можно взять коллекцию MOD11A2.

Мониторинг сельскохозяйственной засухи и ее последствий подробнее описаны здесь. Выявлять засуху по индексам VCI и TCI — это тоже идея Когана и его коллег.

Код примера: https://code.earthengine.google.com/8fc468a7ca7ee9e893fa2ec741948312


*NDVI — хороший и простой показатель зеленой биомассы, но не идеальный. Он, как и все нормализованные индексы, входит в насыщение в области пиковых значений. Он зависит от почвы. Например, светлая и сухая почва вполне может дать NDVI > 0.3 в отсутствии растений. Но сейчас нам важно, что NDVI — это показатель состояния растений.

#GEE #индексы

Vegetation Health Index
1. VCI

Сегодня мы научимся рассчитывать индекс здоровья растительности — Vegetation Health Index (VHI) — на основе индексов температурного режима (Temperature Condition Index, TCI) и состояния растительности (Vegetation Condition Index, VCI), используя Google Earth Engine.

Все три индекса — VCI, TCI и VHI — служат косвенными показателями засухи, изменения климата и общего состояния растительности. Они показывают, насколько текущие показатели отличаются от типичных для данной местности в это время года.

Код скрипта

Начнем с указания области интереса и периода наблюдений: Смоленская и Брянская области России, а также сопредельная территория Республики Беларусь, июнь 2023 года. Область интереса можно создать в GEE при помощи Инспектора (Inspector) или построить в geojson.io и скопировать координаты в свой код.

Индекс состояния растительности VCI (мы уже касались его здесь) оценивает текущее значение NDVI в сравнении с диапазоном значений NDVI, наблюдавшихся в тот же период в предыдущие годы. VCI показывает, где находится текущее значение между минимальным и максимальным значениями NDVI прошлых лет:

VCI = (NDVI – NDVI_min) / (NDVI_max – NDVI_min)

Чем выше VCI, тем лучше состояние растительности. Как правило, VCI изменяется в пределах от 0 до 1, но если текущий NDVI окажется меньше минимального или больше максимального, то мы получим значения ниже 0 и выше 1 соответственно.

Для значений NDVI возьмем коллекцию MOD13Q1.061 Terra Vegetation Indices 16-Day Global 250m и выберем из нее канал NDVI

var NDVI = ee.ImageCollection('MODIS/061/MOD13Q1').select('NDVI');

Пройдем по коллекции, рассчитывая VCI с помощью функции:

function calcVCI(image) {
  var date = image.date();
  var history = NDVI.filterDate('2000-01-01', ee.Date.fromYMD(date.get('year').subtract(1), 1, 1)); 
  var min = history.min();
  var max = history.max();
  var vci = image.subtract(min).divide(max.subtract(min));
  return vci.rename('VCI').copyProperties(image, ['system:time_start','system:time_end']);
}

#GEE #индексы

Vegetation Health Index
2. TCI

Индекс температурного режима TCI позволяет сравнить текущее значение температуры поверхности суши (land surface temperature, LST) с диапазоном значений, наблюдавшихся в тот же период в предыдущие годы. Индекс показывает, где находится наблюдаемое значение температуры поверхности относительно крайних значений (минимума и максимума) в предыдущие годы:

TCI = (LST_max – LST) / (LST_max – LST_min)

Низкие значения TCI показывают, что температура поверхности близка к максимальной для данной местности, а высокие — что к минимальной. Диапазон изменений TCI — от 0 до 1, с теми же оговорками, что и для VCI.

Температуру берем из коллекции MOD11A2.061 Terra Land Surface Temperature and Emissivity 8-Day Global 1km

var LST = ee.ImageCollection('MODIS/061/MOD11A2').select('LST_Day_1km');

Периодичность данных в коллекциях отличается: NDVI вычисляется каждые 16 суток, а LST — каждые 8 суток. Нам нужно получать данные LST каждые 16 суток. Для этого мы перебираем снимки MOD11A2, сделанные за период создания одного композита NDVI (16 суток) и вычисляем среднее значение LST с помощью ee.ImageCollection.mean().

function convert_lst_dates_to_modis_dates(ndvimg) {
  var start = ndvimg.get('system:time_start');
  var end = ndvimg.get('system:time_end');
  var composite = LST.filterDate(start, end).mean();
  return composite
    .set('system:time_start', start)
    .set('date', ee.Date(start).format())
    .set('system:time_end', end)
    .set('empty', composite.bandNames().size().eq(0));
}

После этого, рассчитываем TCI аналогично VCI.

Как видно из рисунка, июнь выдался жарким, хотя и не рекордным.

#GEE #индексы

Vegetation Health Index
3. Объединение коллекций и вычисление VHI

Теперь, зная VCI и TCI, мы можем объединить эти коллекции для расчета VHI. Последний вычисляется по формуле:

VHI = α * VCI + (1 - α) * TCI,

то есть является взвешенным средним значений VCI и TCI. Повышенная температура и плохое состояние растительности (низкий NDVI) являются признаками наступившей засухи, что вполне логично.

Для простоты, примем α = 0.5.

Объединяем коллекции изображений VCI и TCI, и вычисляем VHI, как показано здесь

var filter = ee.Filter.equals({leftField:'system:time_start', rightField:'system:time_start'});
var join = ee.Join.saveFirst('match');

var bothCol = ee.ImageCollection(join.apply({primary: VCI,
                                             secondary: TCI,
                                             condition: filter}))
  .map(function(img) { return img.addBands(img.get('match')).set('date', img.date().format('YYYY_MM_dd')) });

var VHI = bothCol.map(function(img){
  return img.addBands(
  img.expression('a/2 + b/2', {
  'a': img.select('VCI'),
  'b': img.select('TCI'),
  }).rename('VHI'));
}
);

Код выполняет левое внешнее объединение коллекций VCI и TCI с помощью функции ee.Join.saveFirst(). В результате создается коллекция, содержащая пары изображений VCI и TCI, у которых совпадают свойства 'system:time_start'.

К полученной коллекции применяется функция map(). Она используется для добавления к каждому изображению нового канала (ee.Image.addBands()), содержащего совпадающее изображение из другой коллекции. Используем ee.Image.get() для получения совпадающего изображения.

После этого выполняем еще один map(), чтобы вычислить VHI с помощью ee.Image.expression().

В итоге, получим коллекцию изображений, каждое из которых содержит каналы VCI, TCI и VHI.

Код скрипта

Индексы VCI и TCI предложены Феликсом Коганом в работе: Kogan, F. N. (1995). Application of vegetation index and brightness temperature for drought detection. Advances in Space Research, 15(11), 91–100. https://doi.org/10.1016/0273-1177(95)00079-t Там же была предложена формула для вычисления VHI, хотя название индекса Коган придумал позже.

Первоначально, индексы использовались для наблюдений, сделанных прибором AVHRR — Advance Very High Resolution Radiometer, летавшим на спутниках семейств GOES, METEOSAT, MTSAT и DMSP. Сейчас эти данные находятся на странице STAR - Global Vegetation Health Products. Они являются глобальными, начинаются с 1981 года, имеют пространственное разрешение 4 км (Very High Resolution того времени!) и используются в виде 7-суточного композита.

#GEE #индексы

SAVI. Часть 2

Вернемся к формуле

NDVI = (NIR – Red) / (NIR + Red + 2ε).

‍Обратите внимание, что ε не влияет на числитель NDVI. Если бы мы использовали в качестве индекса разность (NIR – RED), у нас бы не было проблем с изменением цвета почвы. Однако, возникла бы другая сложность: новый вегетационный индекс будет меняться в зависимости от освещенности. Когда общая интенсивность света составляет 50% от нормы (например, из-за умеренной облачности), (NIR – RED) тоже будет на 50% меньше. Отсюда и необходимость нормализации индекса путем деления на общую интенсивность (NIR + RED).

Чтобы уменьшить чувствительность формулы к ε, добавим в знаменатель константу L:

Стабилизированный NDVI = (NIR – Red) / (NIR + Red + 2ε + L).

‍‍Разберемся с тем, как выбрать L. Если L будет очень велико по сравнению с NIR и RED, то формула для вычисления стабилизированного NDVI превратится в

(NIR – Red) / (NIR + Red + 2ε + L) ≈ (NIR – Red) / L

Таким образом, для больших L чувствительность к ε исчезает. Зато формула сводится к масштабированной версии (NIR – RED), которая, как мы знаем, не работает.

Напротив, при L=0, мы вернемся к формуле NDVI, которая слишком чувствительна к ε.

Существует компромисс: выберем L достаточно большим, чтобы знаменатель был чувствителен к ε, но не слишком большим, чтобы не исчез эффект нормализации от деления на (NIR + RED).

Для конкретного региона можно определить свое значение L, откалибровав его на основе реальных данных. Однако в целом исследователи пришли к выводу, что L=0,5 создает вегетационный индекс, значительно менее чувствительный к цвету почвы, чем NDVI, и при этом достаточно нормализованный, чтобы общие изменения в интенсивности света не привели к существенному изменению индекса.

‍Осталось сделать последний шаг, чтобы сделать новый индекс удобным. Мы хотим, чтобы его значения находилась в интервале между –1 и 1. Увеличив знаменатель на L, мы добились того, что максимальное значение индекса (теоретически встречающееся при NIR=1, RED=0) стало равно 1/1+L. Если мы хотим, чтобы максимальное значение было равно 1, нам нужно умножить числитель на (1+L):

SAVI = (1 + L)(NIR – Red) / (NIR + Red + L), где L = 0,5

Новый индекс называется Soil-Adjusted Vegetation Index или SAVI. Он предложен в работе: Huete, A. R. (1988). A soil-adjusted vegetation index (SAVI). Remote Sensing of Environment, 25(3), 295–309. https://doi.org/10.1016/0034-4257(88)90106-x

SAVI варьируется в пределах от –1 до 1, малочувствителен к цвету почвы, нечувствителен к общей интенсивности света и увеличивается пропорционально состоянию и плотности растительности. Он лучше NDVI подходит для районов с разреженной растительностью, где цвет почвы может меняться.

#индексы #сельхоз #основы

⭐️ СТРАНЫ / КОМПАНИИ / СПУТНИКИ

Страны: #австралия #германия #индия #иран #испания #канада #китай #португалия #россия #США #япония и т. п.
Но:
#корея обозначает Северную и Южную Кореи
#РБ — Республика Беларусь
#UK — Великобритания

Компании: #planet #maxar

Спутники: #landsat #sentinel1 #sentinel2

⭐️ ДЗЗ

Методы и приборы
#альтиметр
#гиперспектр — гиперспектральная оптическая съемка
#лидар
#оптика — мультиспектральная оптическая съемка
#радиометр — микроволновой радиометр
#dnb — ночная съёмка (day / night band)
#SIF — солнечно-индуцированная флуоресценция хлорофилла
#ro — радиозатменный метод
#SAR — радарная съемка
#InSAR — радарная интерферометрия
#LST — съемка в тепловом инфракрасном диапазоне
#GNSSR — ГНСС-рефлектометрия
#sigint — радиоэлектронная разведка

Виды орбит: #ГСО — геостационарная, #VLEO — сверхнизкая

#основы — обучающие материалы по ДЗЗ
#обучение курсы, обучающие сервисы и т. п.
#история — в основном, история ДЗЗ
#индексы — спектральные индексы
#комбинация — комбинации каналов

Данные
#данные — коллекции данных ДЗЗ, наземных данных, карты и т.п.
#датасет — набор данных для машинного обучения
Дополнительные хештеги, описывающие данные:
#LULC — Land Use & Land Cover
#осадки
#SST — Sea Surface Temperature
#nrt — (near real time) изображения, получаемые в режиме, близком к реальном времени
#debris — космический мусор
#границы — административные границы
#DEM — цифровая модель рельефа (ЦМР)
#keyhole — рассекреченные снимки разведспутников

Литература, справочная информация
#справка — спектральные каналы, орбиты спутников, поиск данных и т.п.
#обзор
#книга — текст книги прикреплён к сообщению.
Дополнительные хештеги:
#наблюдение — ресурсы для наблюдения спутников и орбиты спутников
#космодромы

#конференции — анонс конференций/семинаров/школ, посвященных ДЗЗ и анализ их материалов.
#конкурсы — анонс конкурсов/чемпионатов/олимпиад.
#МВК — материалы заседаний Межведомственной комиссии (МВК) по использованию результатов космической деятельности.

#снимки — поучительные (хоть в чем-то интересные) снимки, первые снимки

Программные инструменты / Языки
#нейронки #софт #GEE #R #tool #python
#ИИ #FM — Foundation Model (Remote Sensing Foundation Model)

⭐️ ОТРАСЛИ / ТЕМАТИЧЕСКИЕ ЗАДАЧИ

#археология #атмосфера #вода #война #засуха #климат #лед #лес #нефть #океан #оползни #наводнение #пожары #почва #растительность #севморпуть #сельхоз #снег
#AGB — надземная биомасса
#ЧС — мониторинг стихийных бедствий и катастроф
#GHG — парниковые газы
Отдельные газы: #CO2 #NO2
#энергетика — космическая энергетика
#SSA — Space Situational Awareness

Метод картографирования цветения рапса

(d’Andrimont et al., 2020) предложили метод картографирования цветения масличного рапса на основе временных рядов оптических данных Sentinel-2 (S2) и радарных данных Sentinel-1 (S1). Район исследования включал в себя северные регионы Германии (N) и южную Баварию (S). Метод использует нормализованный разностный индекс жёлтого цвета (Normalized Difference Yellow Index, NDYI) для S2 (помните, как выглядит цветущий рапс на снимках из космоса?) и локальный минимум коэффициентов обратного рассеяния в поляризации VV для S1. Пик цветения определялся с точностью от 1 до 4 суток. При определении цветения по данным S1 наблюдалась систематическая задержка на 1 сутки, по сравнению с результатами по S2.

📸 Пространственно-усредненные временные ряды S2 NDYI для всех участков на севере и юге Германии с медианными датами начала, пика и конца цветения, полученными по наземным наблюдениям. Пунктирная красная линия соответствует медианной дате начала цветения (BBCH61), пунктирная синяя линия — медианной дате пика цветения (BBCH65), а пунктирная розовая линия — медианной дате окончания цветения (BBCH69).

#сельхоз #индексы

rsi — загрузка данных из STAC и расчет спектральных индексов [ссылка]

Пакет rsi (от repeated spatial infelicities) предоставляет пользователю:

- Интерфейс к проекту Awesome Spectral Indices project, который содержит список спектральных индексов в виде таблицы tibble.
- Метод эффективного вычисления этих спектральных индексов.
- Метод загрузки данных с любого сервера STAC, с дополнительными настройками для загрузки популярных данных Landsat, Sentinel-1 и Sentinel-2 с бесплатных и публичных серверов STAC.
- Метод объединения нескольких растров, содержащих различные наборы данных, в единый растровый стек.

Функция spectral_indices() возвращает таблицу спектральных индексов.

Функция get_stac_data() позволяет загружать изображения из любого доступного каталога STAC. Например, можно загрузить композит каналов Landsat с маской облачности:

aoi <- sf::st_point(c(-74.912131, 44.080410))
aoi <- sf::st_set_crs(sf::st_sfc(aoi), 4326)
aoi <- sf::st_buffer(sf::st_transform(aoi, 5070), 1000)

landsat_image <- get_stac_data(
  aoi,
  start_date = "2022-06-01",
  end_date = "2022-06-30",
  pixel_x_size = 30,
  pixel_y_size = 30,
  asset_names = c("red", "blue", "green"),
  stac_source = "https://planetarycomputer.microsoft.com/api/stac/v1/",
  collection = "landsat-c2-l2",
  mask_band = "qa_pixel",
  mask_function = landsat_mask_function,
  output_filename = tempfile(fileext = ".tif"),
  item_filter_function = landsat_platform_filter,
  platforms = c("landsat-9", "landsat-8")
)

Для популярных данных, например для снимков Landsat, есть отдельные функции, где большинство параметров настроено по умолчанию:

landsat_image <- get_landsat_imagery(
  aoi,
  start_date = "2022-06-01",
  end_date = "2022-06-30",
  output_filename = tempfile(fileext = ".tif")
)

По умолчанию, данные загружаются из Microsoft's Planetary Computer API.

Теперь на основе полученных каналов снимков Landsat рассчитаем спектральные индексы при помощи calculate_indices():

indices <- calculate_indices(
  landsat_image,
  available_indices,
  output_filename = tempfile(fileext = ".tif")
)

Наконец, в rsi есть утилита для эффективного объединения растров, содержащих различные данные об одном и том же месте, в VRT, что позволяет программам типа GDAL рассматривать эти отдельные источники данных как единый файл.

Например, мы можем объединить наши снимки Landsat с полученными индексами:

raster_stack <- stack_rasters(
  c(landsat_image, indices),
  tempfile(fileext = ".vrt")
)

#R #индексы