Опыт классификации космоснимка Landsat с помощью Semi-Automatic Classification Plugin в QGIS
по адресу http://gis-lab.info/qa/landsat_qgis_scp.html
Данная статья описывает опыт работы с Semi-Automatic Classification Plugin для QGIS для классификации снимка Landsat с целью выявления лесонарушений на примере НП "Орловское полесье", а также содержит пошаговую инструкцию для лесного дешифрирования снимка с помощью данного плагина.
Карпачев Андрей Петрович, научный сотрудник НП "Орловское полесье"
Цели проекта
1) Выявление участков усохшего и нарушенного леса на территории НП «Орловское полесье».
2) Общее дешифрирование фрагмента территории Орловско-Брянско-Калужского региона.
Для реализации классификации были выбраны тестовые участки на основе материалов лесоустройства. Примеры этих участков показаны на рисунке ниже.
Определение и формирование макроклассов, проведение классификации
1. Подготовка
а) На этом этапе мы должны подготовить область (тестовый участок) на базе космоснимка Landsat 8 за сентябрь 2014 года на территорию Орловско-Брянско-Калужского региона включая границы НП «Орловское полесье». Космоснимок можно "собрать" непосредственно в QGIS, как показано здесь. Для данного проекта нам понадобятся 6 диапазонов: R, G, B, ближний ИК, коротковолновый ИК-1 и ИК-2. В работе необходимо использовать изображение, полученное из шести каналов (2-7), иначе плагин отказывается строить спектрограммы (однако, классификация выполнена будет), а этот момент важен для последующих расчётов и анализа.
2. Начало работы
Полученный растр добавляем в QGIS с помощью соответствующей кнопки. 
Чтобы добавить наш космоснимок непосредственно в плагин, нажимаем на кнопку "обновить" в плагине (рисунок ниже).
Для более точной визуализации объекта наблюдения можно использовать функцию "RGB: представить снимок в синтезированных и натуральных цветах".
3.Формирование директорий
После загрузки космоснимка нам необходимо сформировать директории, куда впоследствии запишутся файлы ROI и файлы сигнатур.
Переходим в раздел интерфейса плагина, отвечающего за регионы интереса (ROI):
Нажимаем на кнопку "New shp" и задаём на рабочем столе папку с названием ROI, затем открываем и сохраняем документ с названием ROI (сюда будут записываться наши регионы интереса). Переходим в раздел интерфейса плагина, отвечающего за формирование сигнатур и классификацию.
Нажимаем на кнопку "SAVE" и на рабочем столе создаем папку SIG, в которой будет находиться документ с названием SIG. Здесь будут храниться наши сигнатуры для дальнейшей классификации.
4. ROI и SIG
Приступим к формированию файлов "регионов интереса" и сигнатур (рисунок 4).
а) Для объектов с большой площадью, таких как зеркала озёр, используем кнопку полигонального захвата пикселей.
После задаём выбранному объекту id и название и нажимаем на "SAVE ROI".
После записи ROI автоматически изменится и лист сигнатур.
б) Теперь перейдём к захвату пикселей, характеризующих застройку населенных пунктов, и используем коэффициент радиуса захвата.
В строчке "Range radius" меняем коэффициент на 2000, нажимаем на кнопку с плюсом и переходим к захвату пикселей на снимке.
Даём id и название элементу ROI (2, Built-up / 2, Buildings) и нажимаем "save roi" (снова произойдёт отображение в листе сигнатур).
в) Теперь выделим участки голого грунта (с/х угодия) таким же образом. Ввиду большого разнообразия цветности открытого грунта (в примере от фиолетового до слабого марганцового) увеличим радиус захвата до 3500. Выбираем оптимальный вариант отображения.
Даём id и название элементу ROI (3, Bare soil / 3, sh) и нажимаем "save roi" (снова произойдёт отображение в листе сигнатур).
г) Теперь мы переходим непосредственно к главной задаче проекта: поиску и идентификации леса и участков нарушенного леса (Vegetation:Veg). В данном варианте нам нужно разбить территорию леса на хвойную составляющую (h), лиственную (l) и нарушенную (d: defect).
Начнём с нарушенного леса. Выберем коэффициент радиуса захвата равным 550.
Даём id и название элементу ROI (4, Veg / 4, d) и нажимаем "save roi" (снова произойдёт отображение в листе сигнатур). Далее таким же образом захватываем тёмно-зеленые пиксели, характерные для хвои (range radius = 1000), и светло-зелёные, характерные для листвы (range radius = 2000) Лиственный и хвойные лес в данном случае были отнесены к разным ROI путём изменения параметров записи регионов интереса и сохранения по порядку, т.е. 5, Veg / 5, h и 6, Veg / 6, l соответственно.
д) Переходим в лист сигнатур.
По мере записи ROI в листе SIG сформировались строчки с соответствующими названиями и классами регионов интереса, а также задался цвет. Теперь мы должны изменить цвет для финальной визуализации классифицированного изображения. Изменяем цвета примерно вот так:
5. Классификация
Теперь выбираем метод классификации (в разделе интерфейса плагина, отвечающего за сигнатуры), выбираем "Spectral Angle Mapping", "Size = 500" и наводим курсор на космоснимок.
После нажимаем на кнопку "perform classification", создаем папку "classification" и называем в ней документ "classification.tif". Классификация выполнена.
Хочется отметить следующее: если для демонстрации результатов классификации более ничего не требуется (т.е. визуализация проекта отвечает реальным полевым объектам), то дальнейшую более углубленную обработку проводить необязательно. Что касается пиксельного захвата – то можно посоветовать пользоваться кнопкой "show", т.е. определить пиксели, отметить их и далее посмотреть, как они ложатся и не втягиваются ли лишние. Поэтому с коэффициентом захвата можно "поиграть" до определения оптимального отображения (не обязательно использовать числа, указанные в руководстве). Также нелишним, если позволяет снимок (т.е. без проведения атмосферной коррекции), будет отметить дороги и травянистые сообщества. Ошибки, связанные с втягиванием пикселей в какой-то другой класс, предлагается не исправлять, а добавить новый (по очереди в соответствии с id) класс, корректно отображающий объект. Появление жёлтой строчки в QGIS (предупреждение) будет свидетельствовать о пересечении выбранных классов; проверку можно осуществить по графику спектральных кривых. Тот класс, который получился ошибочным, предлагается просто забелить (изменить цвет на белый): так как в проекте нас интересует нарушенная и ненарушенная растительность, белым у нас будет: почва, застройка; воду оставляем без изменений.
Конкретизация и корректировка результатов классификации; представление результатов, вывод
Для этого раздела нам нужно будет подготовить shp-файл (квадрат) тренировочной территории, т.е. часть территории, на которой мы будет проводить дальнейшую классификацию. Шейп-файл лучше сразу же подготовить в QGIS; обратите внимание на проекцию, используемую в проекте (в данном проекте проекция EPSG: 32636, WGS 84/ UTM zone 36N).
Загрузка исходного снимка Landsat. Атмосферная коррекция
а) Вызовем панель загрузки снимка иконкой 
или так:
Далее мы оказываемся на вкладке параметров снимка, здесь нам предлагается создать или выбрать директорию загрузки снимка. Сделать это можно, нажав на кнопку "Select Landsat database directory" и выбрав путь к папке (или создать). Hазовём папку "LandsatDB", нажмём "OK". Теперь ставим галочку в графе "Only Landsat 8" и жмём "Update database" для обновления (ОК). Далее мы должны ввести Image ID (смотрим из предыдущей части анализа на номер архива космоснимка, это оно и есть; в данном проекте это "LC81790232014255LGN00"). Продвинутый пользователь может внести координаты интересующей территории. Нажимаем "Find images". В "image list" появится строчка с записью о космоснимке. Нажимаем на "display image preview" для отображения снимка.
При увеличении масштаба мы видим: качество детализации низкое. Теперь нам необходимо задать параметры загрузки космоснимка. Карточка параметров загрузки должна выглядеть вот так (в нашем примере):
Нажимаем на "download images from list" (выбираем рабочий стол) и ждём загрузки. В зависимости от особенностей вашего компьютера и подключения к Интернету загрузка может занять разное (иногда значительное) время. Ход загрузки мы можем наблюдать в интерфейсе QGIS.
Конец загрузки будет сопровождаться звуковым эффектом. Удаляем из QGIS слой космоснимка.
Нажимаем на иконку 
, в открывшейся форме заходим на вкладку "Landsat". Жмём на кнопку "Select directory" и на рабочем столе ищем и открываем папку со снимком (ОК) – в форме откроются строчки с каналами. Далее ставим галочку в "Apply DOS1 atmospheric correction". Форма должна иметь следующий вид:
Жмём на "perform conversion"; в папке со снимком создаём папку с названием "REF" (ok-ok). Ждём окончания процесса (звуковой сигнал).
В данном проекте выполнение атмосферной коррекции имеет характер проверки (демонстрационный характер). Было необходимым проверить работу данной функции плагина.
В QGIS отобразятся 7 каналов с префиксом "RT". Далее нам будет необходимо открыть заранее подготовленный shp-файл тренировочной территории и произвести прикрепление этого квадрата к снимку (clip).
Кнопкой 
добавляем shp-файл (в текущем проекте этот файл назван "test_location").
В форме плагина нажимаем на вкладку "Clip multiple rasters".
Далее нажимаем на "refresh list" – "select all".
Ставим галочку на "use shapefile for clipping", жмём "refresh list" - в строчке отобразится шейп-файл.
Теперь нажимаем "Clip selected rasters" (на рабочем столе создаем папку "clip"; ок-ок) – звуковой сигнал.
Полученный вид в QGIS:
Нажимаем на любой слой с префиксом "clip_RT", правой кнопкой выходим в свойства и нажимаем "увеличить до слоя". Таким образом, мы ограничиваем территорию исследования. Далее мы отключаем и удаляем все слои с префиксом "RT".
Следующим этапом мы создаем "band set", т.е. набор каналов. Нажимаем на иконку 
(открывается форма) – нажимаем "select all", потом "add rasters to set". Справа будут кнопки со стрелочками (control bands) – ими можно передвигать порядок каналов.
В строчке "quick wavelength settings" выбираем Landsat 8 Oli [bands 2, 3, 4, 5, 6, 7].
В плагине у нас теперь отобразился набор каналов, и в дальнейшем процессе в строчке RGB мы сможем составить свою комбинацию каналов для наилучшей визуализации объектов (ну или воспользоваться предложенными).
Теперь, после настройки дополнительных параметров, мы можем продолжить классификацию. Иконкой открываем shp-файл из папки ROI (из предыдущего раздела классификации) и загружаем регионы интереса. В разделе интерфейса плагина, отвечающем за регионы интереса, нажимаем на кнопку 
выбираем файл. В ROI list отобразились записи.
В разделе интерфейса плагина, отвечающем за сигнатуры, нажимаем "open", ищем папку из предыдущего раздела классификации (SIG) и из неё открываем XML-файл. В списке отобразятся сигнатуры.
Теперь настроим цвета композита (в примере каналы 5-4-6). После окраски визуально определим травянистые сообщества и, выбрав оптимально отражающий объект (коэффициент захвата), так же, как в предыдущем разделе классификации, создадим и сохраним новый регион интереса. Важный момент: теперь к макроклассу Veg (vegetation) добавим новый ROI (grassland). В "MC ID" ставим "4" (в предыдущей части так мы назвали класс вегетации); в "MC Info" ставим "Veg"; в "C ID" ставим номер по порядку (в данном случае это номер 14, потому что ранее автор проводил несколько раз классификацию для закрепления навыка); в "C Info" вписываем название "grassland" и жмём "Save ROI".
Нажмём на строчку "grassland", далее на кнопку 
. Отобразится спектральная кривая:
Проведя курсором по снимку, мы также сможем отметить некоторые закономерности: если выбрать все строки в листе сигнатур и нажать на кнопку то мы получим спектральный график, который в дальнейшем можно использовать в математическом объяснении результатов.
Следующим пунктом выделим дороги и добавим их к классу застройки "Built-up" (2), в "C ID" ставим номер по порядку - 15, в "C Info" вписываем название "road" и жмём "Save ROI". В листе ROI при нажатии "C ID" стрелочками можно поправить нумерацию.
Осталось финализировать классификацию (дополнительную). По желанию (и для поиска наилучшего отображения результата классификации) можно выбрать другой классификационный алгоритм и размер окна. Нажимаем +. В конце процесса нажимаем "perform classification", выбираем папку для сохранения. Данный способ (второй раздел) позволяет нам конкретизировать объекты городской застройки. В ходе данного проекта, когда требовалось провести анализ вегетационной составляющей территории, способ имел ведущее значение. Как отмечалось выше, элементы голой почвы, застройки и некоторые погрешности в определении класса были окрашены белым цветом.
Выводы
Semi-Automatic Classification Plugin для QGIS успешно справился с поставленной задачей классификации космического снимка Landsat с целью выявления усыханий леса. Поэтому ее можно использовать в качестве свободной альтернативы платному специализированному ПО. Кроме того, одним из плюсов плагина является его тесная интерграция с QGIS, что облегчает визуализацию и анализ результатов классификации.
Литература
1.http://semiautomaticclassificationmanual-v4.readthedocs.org/en/latest/Tutorials.html
2.Крылов А. М., Соболев А. А., Владимирова Н. А. Выявление очагов короеда-типографа в Московской области с использованием снимков Landsat //Вестник Московского государственного университета леса - Лесной вестник. – 2011. – №. 4. – С. 54-60
3.Крылов А.М., Владимирова Н.А., Малахова Е.Г. Использование свободных ГИС в системе дистанционного лесопатологического мониторинга // Вестник Московского государственного университета леса - Лесной Вестник, №1 2012. С. 148-152
4.Крылов А.М., Владимирова Н.А. Дистанционный мониторинг состояния лесов по данным космической съемки // "ГЕОМАТИКА" №3(12), 2011 г. с. 53-57