Спутник ДЗЗ
3.11K subscribers
2.43K photos
139 videos
187 files
2.18K links
Человеческим языком о дистанционном зондировании Земли.

Обратная связь: @sputnikDZZ_bot
加入频道
ПРОСТЕЙШИЙ ИНДЕКС ОТКРЫТОЙ ПОЧВЫ
В видимой области спектра отражательная способность почвы выше, чем у растительности. Напротив, в NIR отражательная способность выше у растительности. На этом различии основан новый индекс, предназначенный для разделения этих классов поверхности. Индекс назван NDVISI (Normalized Difference VISible Index) и записывается как

NDVISI = (VIS − NIR)/(VIS + NIR),

где VIS = (Blue + Green + Red), а Blue, Green, Red и NIR – каналы Landsat или Sentinel-2.

Открытая почва имеет отрицательные или очень небольшие положительные значения NDVISI. Классы Urban и Water также имеют более низкие значения NDVISI, чем растительность, но их можно замаскировать при помощью других индексов (MNDWI – для воды, границы городов из OpenStreetMap – для Urban).

Ershov et al. (2022). Natural Afforestation on Abandoned Agricultural Lands during Post-Soviet Period: A Comparative Landsat Data Analysis of Bordering Regions in Russia and Belarus. Remote Sensing, 14(2), 322. https://doi.org/10.3390/rs14020322

#индексы
Дистанционное зондирование Земли (ДЗЗ)
Термины и сокращения, #термины
Организации: NASA, NOAA, DARPA и другие
Спектральные каналы Landsat 8/9 и Sentinel-2, MODIS
Спектральные сигнатуры

📚Основы дистанционного зондирования Земли, #основы
#индексы (спектральные, вегетационные, ...)
#комбинация каналов
#история ДЗЗ
Научно-популярные лекции по ДЗЗ
Лекции школы молодых учёных (ИКИ РАН): 2015-2017, 2018-2019, 2020-2021, 2022-2023
Рекомендованные практики мониторинга ЧС (UN-SPIDER)
Космическое образование в России: раз, два.

Поиск / Справочная информация
Общий каталог искусственных космических объектов (GCAT)
Спутники и съемочная аппаратура
Российские спутники ДЗЗ, #МВК
Информация о запусках
Орбиты спутников
#наблюдение за спутниками
Где взять научную литературу #книга #журнал
ИИ-поиск, патентный поиск, поиск наборов данных
#справка

Google Earth Engine
📚Учебник по Google Earth Engine
🌍 Список всех данных Google Earth Engine
Проекты и примеры кода
Учебные ресурсы
Полезные ссылки
#GEE

📚🖥 Работа с пространственными данными в R

Спутниковые и другие данные#данные
Бесплатные спутниковые снимки, в т.ч. высокого разрешения
🛰 Sentinel-1, Радары на GEE
🛰 Sentinel-2
🛰 Landsat Collection 2, снимки Landsat
🛰 CBERS
#LULC — Land Use & Land Cover
#DEM
#границы
#nrt — Земля из космоса в реальном времени
Международная хартия по космосу и крупным катастрофам: список активаций
Погода: фактическая, реанализ, прогнозы
#ЧС

Тематические задачи
#лес, #AGB (надземная биомасса)
#пожары
#вода — водные объекты, наводнения, качество воды
#лед
#погода, #климат
#атмосфера
#археология
#сельхоз
#LST — температура земной поверхности

Типы данных
#гиперспектр
#SAR #InSAR
#лидар
#LST
#GNSSR
#ro
#SIF

Конференции, школы, семинары
#конференции

Конкурсы и чемпионаты
#конкурс

Новости военного ДЗЗ
#война #sigint #SSA

⭐️Все хештеги
ИНДЕКС ПОЖАРНОЙ ОПАСНОСТИ FWI
Индекс пожароопасной погоды Fire Weather Index (FWI) основан на погодных переменных и учитывает влияние влажности топлива и скорости ветра на возникновение и распространение лесного пожара. Чем выше FWI, тем более благоприятны метеорологические условия для возникновения пожара.

Всего в мире существует около полутора десятков погодных индексов пожарной опасности (в России принят индекс Нестерова). FWI разработан в Канаде и получил широкое распространение во всем мире. Так, его используют европейские системы информации о лесных пожарах EFFIS и GWIS .

Для оперативного мониторинга пожарной опасности можно применить:

* GEOS-5 (температура, влажность, скорость ветра) + GPM IMERG Late (осадки)
* GEOS-5 (температура, влажность, скорость ветра, осадки)

Первый продукт имеет самое высокое пространственное разрешение среди ему подобных (0.1°) и подходит для оценки текущей ситуации. Во втором продукте реализован прогноз FWI на 7 суток (при разрешении 0.25°). Продукты поставляются в файлах netCDF, FWI является одним из датасетов.

#данные #пожары #индексы
Вегетационные индексы

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

Расчет многих вегетационных индексов базируется на двух наиболее стабильных участках кривой спектральной отражательной способности растений: на красную область спектра (0.62–0.75 мкм) приходится максимум поглощения солнечной радиации хлорофиллом, а на ближнюю инфракрасную область (0.75–1.3 мкм) — максимальное отражение энергии клеточной структурой листа. Таким образом, высокая фотосинтетическая активность, связанная, как правило, с большой фитомассой растительности, ведет к более низким значениям коэффициентов отражения в красной области и большим высоким — в ближней инфракрасной.

Оценить фотосинтетическую активность мы можем, сравнивая отражение в красном (Red) и в ближнем инфракрасном (NIR) каналах. Например, можно построить простое отношение каналов

SR = Red/NIR,

то есть значение каждого пикселя канала Red делится на значение соответствующего пикселя канала NIR. При больших объемах фотосинтезирующей биомассы SR будет иметь низкие значения, а при малых — высокие.

Вегетационный индекс создается так, чтобы он был связан с выбранным биофизическим параметром растительности (например, биомассой) и обеспечивал максимальную чувствительность к изменениям этого параметра. Например, можно оценивать объем фотосинтезирующей (то есть зеленой) биомассы, сравнивая отражение в красном и в зеленом каналах. Однако разница между коэффициентами отражения в этих каналах гораздо меньше, чем между красным и инфракрасным каналами, а значит чувствительность такого индекса будет ниже. Связь между индексом и биофизическим параметром в идеале должна быть линейной.

Вместо отношения каналов, для оценки фотосинтетической активности можно использовать разность: NIR - Red. Она будет тем больше, чем выше фотосинтетическая активность. Однако значения такого индекса от участка к участку будут сильно отличаться. Удобно, поэтому нормализовать эти значения, то есть разделить разность на сумму NIR + Red:

NDVI = (NIR - Red) / (NIR + Red).

Значения NDVI всегда находятся в диапазоне [-1; 1], что очень удобно для сравнения разных участков. Оказалось, что использование нормализованной разности между минимумом и максимумом отражений позволяет также уменьшить влияние различий в освещенности снимка, облачности, дымки, поглощение радиации атмосферой и пр.

NDVI расшифровывается как Normalized Difference Vegetation Index — нормализованный разностный вегетационный индекс. На сегодня, это самый популярный вегетационный индекс.

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

* Index DataBase: A database for remote sensing indices
* Alphabetical List of Spectral Indices

В первом содержатся все-все-все индексы, во втором — только самые популярные.

Примеры применения NDVI: http://www.geol.vsu.ru/ecology/ForStudents/4Graduate/RemoteSensing/Lection06.pdf

Индексы применяются не только для оценки параметров растительности, но и при анализе других явлений. В этом случае они называются не вегетационными, а спектральными.

#основы #индексы
GEE-9. Индексы. Арифметика каналов

“Работа, работа — иди в штат Дакота. Оттуда — в Небраску, с Небраски — на Аляску“. Если что, мы находимся в Небраске. Загрузим два снимка этого штата:

var image  = ee.Image("COPERNICUS/S2_SR/20220711T171859_20220711T172712_T14TNL");
var image2 = ee.Image("COPERNICUS/S2_SR/20220721T171859_20220721T172912_T14TNL");

Вычислим нормализованный разностный вегетационный индекс NDVI. Для этого есть специальная функция normalizedDifference()

var ndvi = image.normalizedDifference(['B8', 'B4']);

Чем выше NDVI, тем больше в пикселе здоровой (фотосинтезирующей) растительности. Водные объекты имеют отрицательный NDVI.

А теперь рассчитаем один из “водных” индексов — Normalized Difference Water Index (NDWI). Вычисляют его по формуле:

NDWI = (Green - NIR) / (Green + NIR).

По сути, это перевернутый NDVI. Вода в нем имеет высокие значения, а растительность — низкие (отрицательные). Вычислить его можно аналогично NDVI:

var ndwi = image.normalizedDifference(['B3', 'B8']);

Но мы не ищем легких путей, а рассчитаем NDWI так, будто не знаем о существовании normalizedDifference():

var ndwi = (image.select('B3').subtract(image.select('B8'))
.divide(image.select('B3').add(image.select('B8'))));

Математические операции с каналами задаются функциями: add (+), subtract (-), divide (/). Нетрудно догадаться, что есть еще multiply (*). Полный список операций здесь.

Добавим индексы к одному из снимков:

function addInd(image) {
var ndvi = image.normalizedDifference(['B8', 'B4']).rename('NDVI');
var ndwi = image.normalizedDifference(['B8', 'B4']).rename('NDWI');
return image.addBands([ndvi, ndwi], null, true);
}

image = addInd(image);

Обратите внимание на переименование каналов с помощью rename(). Без этого, каналы будут иметь ничего не говорящие имена nd, nd_1, …

А теперь рассчитаем индекс EVI. Он использует сразу три канала, так что normalizedDifference() не поможет. Запись через add/subtract и т. п. будет слишком длинной. К счастью, в Earth Engine можно записать операции над каналами в виде математических выражений

var evi = image.expression(
'2.5 * ((NIR - RED) / (NIR + 6 * RED - 7.5 * BLUE + 1))', {
'NIR': image.select('B5'),
'RED': image.select('B4'),
'BLUE': image.select('B2')
});

То есть, сначала записали выражение, а потом описали обозначения входящих в него каналов.

Можно выполнять математические операции над разными снимками. Например, вычислим разность значений каналов B8 (NIR):

var difB8 = (image2.select('B8')).subtract(image.select('B8'));

Создадим новый снимок из подготовленных нами индексов.

var newBands = [ndvi, ndwi, evi, difB8];
var image3 = ee.Image(newBands).rename(['NDVI','NDWI','EVI','DIF']);

Теперь у нас будет снимок, состоящий из каналов NDVI, NDWI, EVI и DIF.

Код примера

#GEE #индексы
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 #индексы
Soil-Adjusted Vegetation Index (SAVI)

Нормализованный разностный вегетационный индекс NDVI — это удобная и универсальная мера оценки состояния растительности, но — не идеальная. Когда поле покрыто редкой растительностью, NDVI может колебаться, даже если состояние растительности не меняется. Это происходит потому что почва на поле меняет яркость в зависимости от того, насколько она влажная или сухая.

Напомним формулу NDVI

NDVI = (NIR – Red) / (NIR + Red),

где NIR и Red обозначают отражательную способность в ближнем инфракрасном и в красном диапазонах соответственно.

Предположим, что 20% поля покрыто растительностью, а остальные 80% представляют собой открытую почву 1️⃣. После дождя эта почва станет влажной и, как следствие, более темной. При этом отражение в ближнем инфракрасном и красном диапазонах снизится примерно на одну и ту же величину 2️⃣. В результате NDVI поля увеличится. Напротив, сухая почва становится светлее, и это приводит к уменьшению NDVI 3️⃣.

Итак, исходный NDVI поля равнялся 0,21. После дождя он увеличился до 0,25, а для высохшей почвы упал до 0,17. И все это — без изменения состояния растительности!

Посмотрим, как “исправить” NDVI.

Изменение отражательной способности поля при изменении цвета почвы приводит к тому, что отражение в инфракрасном и в красном диапазонах увеличиваются или уменьшаются примерно на одинаковые величины. Предположим, что эти изменения действительно одинаковы, и обозначим величину изменения через ε. Тогда, с изменением отражательной способности почвы, NDVI будет равен

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

или

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

ε зависит от доли открытой почвы на поле и от цвета почвы. Если поле в основном покрыто растительностью и почвы не видно, то ε будет мало по сравнению с (NIR+RED), так что им можно пренебречь и мы получим обычную формулу NDVI. То есть, когда почвы не видно, NDVI не чувствителен к изменениям ее цвета.

Нас же интересует ситуация, когда почву видно хорошо...

#индексы #сельхоз #основы
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 подходит для районов с разреженной растительностью, где цвет почвы может меняться.

#индексы #сельхоз #основы
SAVI. Часть 3

‍Сравним NDVI и SAVI.

На рисунке 1️⃣ показана разница между NDVI и SAVI при различной плотности растительности. При редкой и умеренной растительности NDVI имеет значительный разброс значений из-за чувствительности к цвету почвы. SAVI имеет гораздо меньший разброс. Разброс NDVI уменьшается, когда плотность растительности высока и почва закрыта растительным покровом.

Кроме малой чувствительности к цвету почвы, существует еще одно преимущество, которое SAVI приобретает благодаря дополнительному члену L: порог насыщения SAVI выше, чем NDVI.

Вспомним формулу NDVI

NDVI = (NIR – Red) / (NIR + Red),

Когда Red приближается к нулю, NDVI будет расти вместе с NIR, асимптотически приближаясь к единице. Все большее увеличение NIR оказывает все меньшее влияние на NDVI, то есть происходит насыщение.

Рассмотрим эффект от добавления L к знаменателю

SAVI = (1 + L)(NIR – Red) / (NIR + Red + L)

SAVI тоже достигнет насыщения, но это произойдет при большем значении NIR. Чтобы понять как это работает, зададим RED = 0 и сделаем L гораздо больше, чем NIR. В таком случае, в знаменателе формулы будет преобладать L, а вся формула сводится к (1+L)NIR/L. Она линейно возрастает с увеличением NIR, выпуклости графика нет, а значит, нет и насыщения 2️⃣. Таким образом, точка насыщения увеличивается вместе с увеличением L.

‍Таким образом, дополнительный член в знаменателе позволяет SAVI регистрировать увеличение плотности растительности даже после того, как NDVI вошел в насыщение.

Подведем итоги:

1. SAVI, как правило, используется вместо NDVI в ситуации, когда видна значительная часть почвы и возможны изменения ее яркости.
2. SAVI не имеет однозначного превосходства перед NDVI: смягчая эффекты яркости почвы, в нем нарушается нечувствительность индекса к общей интенсивности света (за счет добавления L).
3. Меньшая чувствительность SAVI к насыщению позволяет использовать этот индекс для оценки состояния густой растительности.

#индексы #сельхоз #основы
⭐️ СТРАНЫ / КОМПАНИИ / СПУТНИКИ

Страны: #австралия #германия #индия #иран #испания #канада #китай #португалия #россия #США #япония и т. п.
Но:
#корея обозначает Северную и Южную Кореи
#РБ — Республика Беларусь
#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).

#сельхоз #индексы
GEE-44. qualityMosaic или самый зелёный пиксель Беларуси

Самым зелёным будем считать пиксель с максимальным значением вегетационного индекса NDVI, который используется как показатель количества зелёной растительности. Необходимо найти максимальное значение NDVI в заданной области в течение календарного года, установить день года, в который оно достигается, и, наконец, построить карту.

Областью исследований стала Республика Беларусь. Выбор спутниковых данных в задаче не важен. Мы используем снимки Landsat 8.

Основная работа заключается в: 1) создании коллекции изображений l8, 2) добавлении к снимкам коллекции слоя NDVI, а также слоя, содержащего номер дня года (doy), в который получен снимок. Самую важную работу выполняет функция qualityMosaic:

var l8 = ee.ImageCollection('LANDSAT/LC08/C02/T1_L2')
.filterBounds(AOI)
.filterDate(start_date, end_date);

function addNDVI(image) {
var ndvi = image.normalizedDifference(['SR_B5', 'SR_B4']).rename('NDVI');
return image.addBands(ndvi);
}

function addDOY(image) {
var img_date = ee.Date(image.date());
var img_doy = ee.Number.parse(img_date.format('D'));
return image.addBands(ee.Image(img_doy).rename('doy').toInt());
}

var greenest = l8.map(addNDVI).map(addDOY).qualityMosaic('NDVI');


Функция ImageCollection.qualityMosaic(qualityBand) формирует из коллекции изображений мозаику (итоговое изображение), состоящую из пикселей с наивысшей оценкой качества, то есть с максимальными значениями слоя qualityBand. У нас в качестве такого слоя выступает 'NDVI', а значит мозаика будет состоять из максимальных (за год) значений NDVI в данном пикселе.

При этом qualityMosaic включает в мозаику не только значения канала qualityBand, но и значения всех остальных каналов снимка с наивысшей оценкой качества. Таким образом, в итоговое изображение попадет и канал 'doy' — день года, в который достигается максимум NDVI.

1️⃣ RGB-мозаика в естественных цветах, 2️⃣ Максимумы NDVI, 3️⃣ Дни года, когда достигается максимум NDVI.

🌍Код скрипта GEE

Мы уже рассказывали о создании мозаик здесь и здесь. Вот ещё немного полезной информации:

🔗 Создание мозаики из коллекции изображений Sentinel-2, где качество пикселя основано на оценке вероятности облаков
🔗 Пример использования qualityMosaic для PythonAPI от Q. Wu

#GEE #индексы
This media is not supported in your browser
VIEW IN TELEGRAM
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 #индексы