Расчёт вегетационного индекса NDVI в геоинформационной системе QGIS

Произведём расчёт популярного вегетационного индекса NDVI в настольной геоинформационной системе QGIS на примере каналов спутниковых данных Landsat 8.

1. Для расчёта индекса NDVI нам понадобится:

2. Пошаговая инструкция

3. Картинки для перчинки

NDVI (Normalized Difference Vegetation Index) — нормализованный относительный индекс растительности — простой количественный показатель количества фотосинтетически активной биомассы (обычно называемый вегетационным индексом). Один из самых распространенных и используемых индексов для решения задач, использующих количественные оценки растительного покрова. Данный индекс используется для количественной оценки биомассы на сельскохозяйственных полях.

1. Для расчёта индекса NDVI нам понадобится:

  1. человек;
  2. компьютер средней мощности;
  3. QGIS посвежее;
  4. прочесть базовые вещи по NDVI
    1. NDVI — {теория} и практика на gis-lab.info;
    2. Функция NDVI — из справки ArcGIS;
  5. формула:
  • NDVI = (NIR — RED) / (NIR + RED)
  • NIR — инфракрасный канал, RED — красный канал.
  1. снимок Landsat 8 в TIF;
  2. цветовая шкала NDVI из вышеупомянутой статьи для QGIS .
NDVI_index_calculation_program_QGIS_ndvi_gis-lab_info
Цветовая шкала NDVI

2. Пошаговая инструкция

  1. скачиваем снимок Landsat 8 в формате TIFF;
  2. распаковываем снимок в одноимённую папку;
  3. распаковываем шкалу ndvi.zip;
  4. запускаем QGIS;
  5. нам нужны каналы NIR и RED, у Landsat 8 это 5-й и 4-й каналы соответственно. Загружаем слои:
  • LC8xxxxxx20xxxxxLGN00B5
  • LC8xxxxxx20xxxxxLGN00B4
  1. делаем вычисления. Меню Растр –> Калькулятор растров,
    1. в поле Выражение пишем формулу NDVI: float («LC8xxxxxx20xxxxxLGN00B5@1» — «LC8xxxxxx20xxxxxLGN00B4@1») / float («LC8xxxxxx20xxxxxLGN00B5@1» + «LC8xxxxxx20xxxxxLGN00B4@1»)
    2. Формат вывода = GeoTIFF;
    3. Слой результатов = имя_получаемого_снимка_NDVI. Можно что-то типа LC8xxxxxx20xxxxxLGN00_NDVI. Полезный совет, если не хотите гадать почему вдруг не срабатывает сохранение файла, тогда избегайте в путях и именах файлов русскоязычных символов и пробелов, и тогда проблем не будет.
    4. [V] Добавить результат в проект;
  2. видим, что получился наш снимок NDVI, отображается он в Grayscale, а нам нужно в цветной удобной шкале увидеть растительность, для этого: заходим в Свойства слоя –> Загрузить стиль… –> [путь_до_папки_с_шаблоном_ndvi]/ndvilinear.qml –> OK;
  3. радуемся результату;
  4. если мы захотим это дело распечатать, то:
    1. меню Проект –> Создать макет –> имя_макета;
    2. на панели инструментов выбираем Добавить карту, помещаем её как надо;
    3. и самое главное — не забудьте поместить шкалу NDVI. Выбираем Добавить легенду;
    4. неплохо бы дать расшифровку значениям, 0.67 — густая растительность, 0.50 — разреженная растительность и т.п.
    5. остальное на свой страх и вкус.

3. Картинки для перчинки

28 Replies to “Расчёт вегетационного индекса NDVI в геоинформационной системе QGIS”

  1. Спасибо за подробную инструкцию. Все получилось. Использовал QGIS 2.6.
    А можете подсказать: я добавил на получившийся NDVI слой точки с навигатора. Как можно узнать конкретные значения NDVI для этих точек? Я уже, вроде, все в программе излазил, но ничего подходящего не нашел.

    1. Поставьте плагин Point Samling Tool (Plugins —> Manage and Install Plugins… —> All —> Point Sampling Tool —> Install Plugin).
      Открываете слой NDVI и слой точек. Запускаете плагин. Выбираете поля которые будут включены в создаваемый слой точек, и здесь же выбираете растр (NDVI). Указываете выходной формат, Ouput point vector layer, —> OK. На выходе получаете shape-файл c точками, у которого атрибутивная таблица состоит из тех полей, что вы выбрали + поле с именем слоя растра, в котором и хранятся нужные числовые значения.

  2. Добрый день! Уже второй раз считаю по вашей инструкции, в первый раз все было отлично, а сейчас, после добавления стиля-легенды, все осталось серым…В чем может быть причина? И еще, у вас написано выражение для калькулятора/, что значит ‘float’?

    1. Добрый день. Раз первый раз всё получилось, значит где-то ошиблись, попробуйте всё заново, может быть шаблон перепутали, нужно выбирать ndvilinear.qml.
      Float значит вещественное число, здесь в выражении мы явно задаём тип числового значения float, для того чтобы быть уверенными что NDVI считается с максимальной точностью.

  3. да, перепутала шаблоны. Спасибо. float-забивается вручную? или это в калькуляторе функция какая-то? Я в калькуляторе ставлю (В5-В4)/(В5+В4) и всё,это неверно?

    1. В калькуляторе мы всё вручную набираем. (B5-B4)/(B5+B4) это формула, а калькулятору мы указываем конкретные растры с которыми работаем.
      Выражение должно выглядеть так:
      float («LC8xxxxxx20xxxxxLGN00B5@1″ — «LC8xxxxxx20xxxxxLGN00B4@1″) / float («LC8xxxxxx20xxxxxLGN00B5@1″ + «LC8xxxxxx20xxxxxLGN00B4@1″).
      Там где LC8xxxxxx20xxxxxLGN00B должно быть ваше название конкретного космоснимка Landsat 8.

      1. это все ясно, меня интересует, «float» вбивается «вручную» или это какая-то функция калькулятора?

  4. Здравствуйте! У меня есть исходные снимки с БПЛА с RGB камеры. Могу ли я их сшить в QGIS 2.18.3 (x64) и рассчитать NDVI и другие снимки?
    Здесь же методика будет другая. Не спутниковый снимок. Как узнать где какой канал?

    1. Методика такая же. Чтобы расчитать индекс NDVI нужен инфракрасный спектр и красный, обычная камера этого не обеспечивает, соответственно нужна такая NIR-камера, которая снимает в инфракрасном и красном спектрах. Второй аспект, космоснимки Landsat чётко ортотрансформированы, устранены геометрические искажения и точно геопозиционированы, а снимки с БПЛА имеют самое различное качество в зависимости от камеры, большие геометрические искажения объектива. Для создания аналога космоснимка — аэрофотоснимка, у которого к слову сказать пространственное разрешение гораздо выше, применяется специальный софт, вроде Pix4D, DroneDeploy. Такой софт готовит полётную миссию — заранее обозначенную область (прямоугольник) съёмки на карте, задаётся метод движения коптера — серпантин и др. Потом эта серия снимков склеивается на серверах этих сервисов, там же производится ортотрансформирование, устранение искажений и т.п. На выходе получается подобие космоснимка — аэроснимок, только с гораздо более высоким пространственным разрешением (до 5 см и выше).
      У меня самого, опыт только с Phantom 4, с помощью указанного софта удавалось создавать геопривязанный аэроснимок высокого разрешения. Соответственно, чтобы это был NDVI-снимок, нужно просто переоборудовать камеру, готовые решения стоят дорого — в цену нового квадрокоптера. Самому интересно, как на фантоме, по простому и за дешёво помимо съёмки в обычном диапазоне, снимать и в инфракрасном.

    1. Географическая система координат является стандартом де-факто для векторных данных и самая популярная это WGS 84 (EPSG:4326).

      Соответственно нужно создать векторный слой с системой координат WGS 84 EPSG:4326 (используется по-умолчанию для векторных данных).

      Потом в свойствах проекта установить формат координат —
      ‘Settings’ > ‘Project Properties’ > ‘General’ > ‘Coordinate display’ — ‘Display coordinates using’ = ‘Degrees, minutes, seconds’ > OK.
      После этого в статусной строке в графе ‘Coordinate’ отобразятся координаты в формате градусы минуты секунды.

      Что касается спутниковых данных, они как правило используют проекционную систему координат UTM.

    2. Спасибо Вам!!!! Я это 3 месяца искал!!… Кто ищет тот найдет..! Спасибо!! Сразу вопрос. А можно снимки сшитые из drone deploy так оценить?

  5. Пожалуйста.
    Да можно, у них на сайте примеры есть (DroneDeploy).

    Просто нужно обеспечить аэросъёмку с квадрокоптера или другого БПЛА в ближнем инфракрасном спектре и в видимом красном.

    Формат данных — стандартный GeoTIFF, когда будут такие данные, просто открываем в QGIS, калькулятор растра, вычисляем NDVI.

    Продаются наборы для переоборудования камер, они очень дорого стоят.

    Есть вариант где камера заменяется, а в другом — цепляется параллельно основной камере.

  6. Здравствуйте!
    Мне нужно определить динамику орошаемых площадей за 12 лет на один бассейн.

    Возможно у Вас была похожая задача, могли бы подсказать.

    Так, как мне нужно просто оконтурить орошаемую площадь я взяла снимки за все 12 лет, месяц где пик растений приходится на (июнь). Я хочу в ArcGISе посчитать NDVI. До этого где то я читала, что нужно сделать предварительную обработку радиометрическую калибровку и атмосферную коррекцию нужно сделать для снимка, потом посчитать NDVI. Или NDVI в таких случаях как у меня считаются на сырой снимок? Использовать я буду Landsat8.
    Как правильно будет сделать? могли бы по шаговую инструкцию написать?

    1. Динамику можно сделать и просто как вы сказали на «сырую», но по фэншую положена атмосферная коррекция. Если её не делать, то просто полученные значения NDVI нельзя будет ни с чем сравнить, данные будут актуальны только в рамках вашей работы.
      Спутниковые данные с атмосферной коррекцией можно получить с сайта earthexplorer.usgs.gov выбрав категорию Landsat Collection 1 Level-2 (On-Demand) > Landsat 8 OLI/TIRS C1 Level-2.
      Извлечь числовые значения из растра, из конкретных точек можно с помощью плагина для QGIS — point sampling tool.
      Про ArcGIS ничего сказать не могу, в нём не работаю.

  7. Здравствуйте! Не могли бы Вы подсказать, как можно получить числовые значения индекса NDVI в точках с помощью плагина для QGIS — point sampling tool. На данный момент построены только карты NDVI по вашему алгоритму.

    1. 1) Открыть QGIS
      2) Установить плагин — point sampling tool.
      меню ‘Plugins’ > ‘Manage and Install Plugins…’ > ‘Search’=’point sampling tool’ > ‘Install plugin’
      3) Открыть снимок NDVI, зайти в свойства, запомнить CRS
      4) Создать точечный слой с той же системой координат (CRS) что и снимок.
      меню ‘Layer’ > ‘Create Layer’ > ‘New Shapefile Layer…’, ‘Geometry type’=»Point», выбрать систему коорд.
      5) На панеле инструментов запустить значок ‘Point sampling tool’.
      В 1й графе выбрать точечный слой, во 2-й из какого слоя точки извлечь (NDVI), указать выходной файл.
      6) В результате получается точечный слой с извлечёнными из растра значениями.

      1. Спасибо большое за оперативный ответ, все очень доступно объяснили! Только еще вопросик, помогите, пожалуйста (так как я с QGIS на вы): если у меня слой NDVI в прямоугольной СК, а точки в географической СК, как их можно привести к одной СК?

      2. 1) В панеле слоёв ‘Layers’ дважды кликнуть по NDVI снимку — попадём в свойства, смотрим и записываем графу CRS.
        2) В панеле слоёв нажимаем правой кнопкой мыши (ПКМ) на точечный слой > ‘Save as…’,
        в графе CRS выбираем запомненную систему координат (СК). Даём имя файла ‘File name’ > ‘…’ —> ‘OK’

      3. А какая причина может быть, что этот новый слой выходит с пустой атрибутивной таблицей, все точки на месте, а в поле значений NDVI пусто?

  8. Атрибутивная таблица получилась пустой. Что можно сделать в данном случае? Cпасибо за ответ.

    1. Если вы про плагин point sampling tool, то пользоваться им очень просто. Нужен слой растра и точек. Слой с точками должен иметь систему координат (CRS) такую же как и растра, об этом кстати сам плагин предупреждает.

Оставьте комментарий