Обработка и интерпретация данных Landsat 8 (OLI) средствами GRASS GIS 7

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Эта страница является черновиком статьи.


Использование открытой геоинформационной системы GRASS GIS 7 для реализации стандартных сценариев обработки данных Landsat 8 (сенсор OLI)

В сообществе с завидной регулярностью появляются вопросы, связанные с обработкой данных Landsat 8 в GRASS GIS. В силу специфики этого программного пакета далеко не каждый способен быстро разобраться в нюансах его работы, и часто пользователь либо бросает затею использования GRASS, либо решает задачи интерпретации данных некорректно.

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

Некоторые предварительные замечания:

  • При работе над статьей использовалась последняя на дату публикации стабильная сборка GRASS GIS 7.0.0. Тестирование осуществлялось в операционных системах Windows 8 (x64), Windows 7 (x32), Mac OS Yosemite, Linux Debian Wheezy.
  • Все (или почти все) действия в статье выполняются через графический интерфейс GRASS, несмотря на то, что часто проще и быстрее использовать командную строку. Это сделано для увеличения наглядности и уменьшения порога вхождения. Основные операции сопровождаются текстом соответствующей команды.
  • Предполагается, что вы загрузили и распаковали набор данных Landsat 8 уровня обработки "Level 1 Product" (Например с помощью EarthExplorer, краткие инструкции можно получить здесь или здесь) с одиннадцатью одноканальными изображениями и файлом метаданных.
  • Не гарантируется корректная работа модулей GRASS, если в путях до ваших файлов и их именах используются кириллические символы, а также если имя пользователя в операционной системе кириллическое.


Создание проекта

Итак, запускаем GRASS. Для начала необходимо создать папку, в которой будет храниться проект. Там будет создана директория с именем набора, которое мы зададим позже, и там же будут храниться все файлы, которые будут использоваться в процессе работы.


Стартовое меню GRASS. Выбор директории для хранения проекта

По умолчанию доступна демонстрационная область "demolocation", данные которой хранятся в домашней директории пользователя, папке grassdata. Рекомендуется сразу выбрать новую папку. После выбора новой директории становятся неактивными все кнопки, кроме "New" у блока с областями. В появившемся по её нажатию окне называем location произвольным подходящим именем (Рекомендуется использовать латинские названия без пробелов). Title может быть произвольным. В примере будет использован набор данных LC80110462013131 (Загрузить), охватывающий юго-восточное побережье Кубы, поэтому используется имя Cuba.


Окно создания новой локации


После нажатия на кнопку "Далее" GRASS предложит выбрать систему координат для создаваемой локации. Правильнее всего будет зачитать её напрямую из наших наборов данных (опция Read projection and datum terms from a georeferenced data file). Можно также указать код EPSG или выбрать проекцию из списка.


Выбор системы координат


GRASS выведет параметры импортированной проекции. Проверяем их на адекватность и нажимаем "Finish".


Описание проекции


После этого появляется предложение сразу импортировать растр в созданный набор. Пока лучше отказаться, хотя утвердительный ответ ничего не сломает. Также будет предложено создать набор геоданных (mapset) и задать регион (region). Набор геоданных создаем, а регион лучше зададим позже, после того как импортируем все необходимые растры.


Импорт данных

Теперь, когда среда подготовлена, можно приступать к работе. Для начала следует импортировать все исходные изображения в рабочий каталог. Используем модуль r.in.gdal (File - import raster data - common formats import).


Функция импорта растровых данных


В открывшемся окне указываем путь и выбираем все нужные нам файлы. Активируем флаг Add imported layers into layer tree (Добавить импортированные слои в дерево слоёв).


Окно импорта растровых данных
r.in.gdal input=G:\data\image_proc\Cuba\input\LC80110462013131LGN01_B1.TIF output=LC80110462013131LGN01_B1


В основном окне программы переходим на вкладку Map layers и видим там все добавленные слои. Отключаем их для быстродействия.


Дерево слоёв и окно просмотра

Для визуализации из командной строки можно воспользоваться командами:

d.mon start=wx0
d.rast map=LC80110462013131LGN01_B1


Далее необходимо задать регион с помощью модуля g.region. Для этого обращаемся в Settings - Region - Set Region.


Функция "задать регион"


Поскольку все изображения одного набора данных Landsat точно соответствуют друг другу геометрически, в качестве региона можно выбрать любой из стандартных каналов (один) в поле [multiple] Set Region to match raster map(s) (Установить регион в соответствии с растровой картой). Оговоримся, что на этом этапе следует выбрать канал с разрешением 30 метров (т.е. не панхроматический).


Меню определения региона
g.region raster=LC80110462013131LGN01_B2


Шаг с выбором региона отдельной функцией можно расценивать как избыточный, так как на этапе импорта данных можно установить флаг "Extend region extents based on new dataset" (-e).

Радиометрическая калибровка

Первым важным шагом при обработке данных Landsat является радиометрическая калибровка. Не нужно путать её с более сложной и не реализуемой пользовательскими силами процедуры радиометрической коррекции. В загруженном наборе данных Landsat каждый пиксель хранит безразмерное нормализованное значение (Digital Number / DN), полученное после преобразований над сырыми значениями, зарегистрированными сенсором спутника. В каждом независимом наборе данных (и даже в разных каналах одного и того же набора) эти числа могут оказываться совершенно несопоставляемыми, они не несут непосредственного физического смысла, поэтому обрабатывать такие данные не принято (только в целях визуального дешифрирования). Однако существуют способы, используя сведения, записанные в метаданных набора, пересчитать DN в один из двух физических параметров - отражательную способность (reflectance) или спектральную энергетическую яркость (radiance). Немного мат. части из официальной документации. Осуществить такие преобразования можно двумя способами: с помощью модуля i.landsat.toar или с помощью калькулятора растров (r.mapcalc). Для того, чтобы убедиться, что i.landsat.toar не совершает ничего магического, Radiance рассчитаем обоими способами, сравнив результаты.

Итак, i.landsat.toar. Модуль работает последовательно с каждым загруженным каналом. Запускаем его через Imagery - Satellite images tools - Landsat DN to radiance/reflectance.


Функция пересчета DN в radiance/reflectance


Во вкладке required (обязательные) задаем basename для исходных и выходных снимков. Для исходных это строка, которая должна полностью совпадать с именами загруженных вами снимков, только без цифры канала, которую модуль будет подставлять автоматически. В рассматриваемом примере: LC80110462013131LGN01_B. Для выходных - произвольная, понятная вам строка. Разумеется, желательно, чтобы по ней был понятен характер хранимой информации. Например, в нашем примере toar_refl_B (результирующие растры будут иметь имена toar_refl_B1, toar_refl_B2 и т.д.).


Вкладка reqiured модуля i.landsat.toar


Во вкладке metadata достаточно указать путь до файла метаданных (В примере: LC80110462013131LGN01_MTL.txt). Параметры сенсора и всё прочее модуль автоматически прочитает из него.


Вкладка metadata модуля i.landsat.toar


Во вкладке optional ключевой флаг отвечает за тип результатов. По умолчанию модуль рассчитывает отражательную способность (reflectance), для расчета же энергетической яркости (radiance) нужно отметить флаг "Output at-sensor radiance isntead of reflectance for all bands".

Примечение: если у вас возникает необходимость перезаписать какие-либо из уже существующих файлов, практически в каждом модуле (во вкладке Optional или где-нибудь еще) доступен флаг "Allow output files overwrite existing files" (в командной строке за это отвечает опция --overwrite). Ещё примечание: флаг "verbose module output", также доступный во вкладке Optional, позволяет получать детализированный отчет о ходе выполнения операции. В командной строке этой цели служит опция --verbose.


Вкладка options модуля i.landsat.toar


Обратите внимание на нижнюю часть окна, там приведена полная команда, которую исполнит GRASS: "i.landsat.toar --overwrite --verbose input=LC81740292013121LGN01_B output=toar_refl_B metfile=/.../LC81740292013121LGN01_MTL.txt". Её запись, если изучить её внимательно, полностью соответствует нашим ожиданиям. После запуска программа будет выполнять последовательное преобразование всех растров до последнего существующего в рабочей области (в нашем случае 8-го панхроматического включительно) канала. Если вы не импортировали какие-то каналы на предыдущих этапах, то модуль сообщит об ошибке: Unable to open header file for raster map... Это не должно вас напугать, поверьте, всё что мы добавили - благополучно обработано.

Обратим также внимание на параметр "esun". В документации к i.landsat.toar указано, что результаты без параметра esun дают приблизительно такие же значения, как и с ним, и сообщение вида "WARNING: ESUN evaluated from REFLECTANCE_MAXIMUM_BAND", как и вышеупомянутая ошибка, не стоит ваших переживаний.


Command output модуля i.landsat.toar


Сейчас мы произвели расчет Reflectance. Для расчета Radiance достаточно переименовать выходные слои в поле "Prefix for output raster maps" (toar_radiance_B) и во вкладке optional поставить выше упоминавшийся флаг "Output at-sensor radiance isntead of reflectance for all bands". Запускаем алгоритм еще раз, получаем искомое.

Теперь обещанные развлечения с калькулятором растров. Согласно документации, Radiance рассчитывается по следующей формуле:


Формула DN to Radiance (Landsat 8)


где, Lλ — количество приходящего излучения; Lmin — количество приходящего излучения, которое после масштабирования становится Qmin; Lmax — количество приходящего излучения, которое после масштабирования становится Qmax; Qcalmin — минимальное калиброванное значение DN (0 или 1); Qcalmax — максимальное калиброванное значение DN (255); Qcal — калиброванное значение рассчитываемого пикселя (собственно DN). Все константы можно найти в файле метаданных, например Qcalmin (Qcalmax) записаны как QUANTIZE_CAL_MIN(MAX)_BAND_channel number, а Lmin (Lmax) как RADIANCE_MIN(MAX)_BAND_channel number.

...
    QUANTIZE_CAL_MAX_BAND_1 = 65535
    QUANTIZE_CAL_MIN_BAND_1 = 1
    QUANTIZE_CAL_MAX_BAND_2 = 65535
    QUANTIZE_CAL_MIN_BAND_2 = 1
    QUANTIZE_CAL_MAX_BAND_3 = 65535
    QUANTIZE_CAL_MIN_BAND_3 = 1
...
    RADIANCE_MAXIMUM_BAND_1 = 744.89569
    RADIANCE_MINIMUM_BAND_1 = -61.51373
    RADIANCE_MAXIMUM_BAND_2 = 762.78229
    RADIANCE_MINIMUM_BAND_2 = -62.99081
    RADIANCE_MAXIMUM_BAND_3 = 702.89734
    RADIANCE_MINIMUM_BAND_3 = -58.04549
...

Теперь откроем калькулятор растров (модуль r.mapcalc), Raster — Raster map calculator.


Запуск калькулятора растров


И, для примера, запишем формулу для второго канала, взяв соответствующие значения из файла метаданных. Сохраним результат в файл mapc_radian_B2.


Ввод формулы в калькуляторе растров


Больше ничего настраивать не нужно. Запускаем и получаем новый растр. Добавим radiance-версии второго канала в дерево слоёв через File - Map Display - Add multiple rasters or vectors.


Функция добавления слоёв на карту


Теперь самое интересное - сравнить результаты! Выделяем интересующие нас слои в дереве (toar_radian_B2 и mapc_radian_B2).


Выделенные слои, которые будут сравниваться


В окне Map Display активируем инструмент "query raster/vector maps".


Инструмент получения информации


Кликаем в произвольном месте рабочей области, в появившемся окне видим, что значения в любой точке у обоих растров идентичные. Что и требовалось доказать!


Окно с выводом значений пикселя


Атмосферная коррекция

Получаемые съемочной системой данные о подстилающей поверхности и объектах на ней, очевидно, являются искаженными из-за влияния множества факторов, среди которых основной - атмосфера, сложная и разнообразная. Ясно, что для корректной обработки необходимо максимально подробно учесть эти искажения, восстановив "истинное лицо" объектов на земной поверхности, именно для этого применяется так называемая Атмосферная коррекция. Чтобы учесть атмосферу, её нужно как-то смоделировать, или, выражаясь корректнее, смоделировать поведение отраженной солнечной радиации на пути от объекта до съемочной аппаратуры. В GRASS (модуль i.atcorr) реализована довольно сложная и подробная модель 6S (Second Simulation of a Satellite Signal in the Solar Spectrum). Документация модели. Ей мы и попробуем воспользоваться.

Существует множество других методов атмосферной коррекции, из которых большим авторитетом пользуется FLAASH, реализованный в ENVI. Для хотя бы примерной оценки полученных в GRASS результатов мы сравним значения NDVI, рассчитанные по данным после коррекций методами 6S и FLAASH.

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


В соответствии с документацией к модулю i.atcorr для начала необходимо создать конфигурационный файл для метода 6S. Для улучшения результатов предварительно добудем данные по рельефу и атмосферным условиям на исследуемую территорию (модель 6S во многом опирается на эти данные, хотя они и являются опциональными).

Данные о рельефе можно брать отсюда (SRTM). Импортируем их как обычный растровый набор данных.

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


В качестве одного из параметров необходимо будет указать среднюю высоту (помимо самой ЦМР). Рассчитаем её с помощью модуля r.univar. Запускаем его через вкладку Command Console:


Запуск модуля r.univar


В появившемся окне выбираем слой с ЦМР


Модуль r.univar. Выбор растрового слоя с ЦМР


Исполняем модуль и отслеживаем результат работы во вкладке command output. Запоминаем (записываем) значение средней высоты (mean).


Результат работы модуля r.univar


Также нам понадобятся координаты центра снимка. Если набор исходных данных был загружен с сайта EarthExplorer, то при выборе снимка указывались географические координаты центра сцены. Можете быстро найти свою сцену по имени.


EarthExplorer. Координаты центра сцены


Создаем файл с параметрами для 6S (например atcorr_pars.txt). Это самый обыкновенный текстовый документ, где одна строка - один параметр (в жестком порядке). Разные каналы корректируются отдельно (то есть процедуру придется повторить несколько раз). Рассмотрим файл конфигурации для нашей Кубинской сцены, синий канал (B2).

18					- читаем следующую строку как положение спутника Landsat8
05 11 15.48 -75.722 20.187		- месяц| день| время GMT| долгота | широта 
2					- атмосферная модель
1					- аэрозольная модель
0					- не указываем видимость
0.120					- параметр оптической толщины аэрозолей для диапазона 550 нм
-0.244					- средняя высота над уровнем моря в км(с обратным знаком)
-1000					- означает, что сенсор находится на борту спутника
116					- индекс соответствующего канала landsat 8 (здесь - синий)


Первая строка это "Geometrical conditions". В сущности, здесь - разные сенсоры. Landsat 7 (ETM+) имеет код 8, а используемый нами Landsat 8 (OLI) - код 18. В зависимости от указанного здесь кода формат второй строки будет отличаться (см. таблицу в документации). Во второй строке (Для Landsat 8) указывается время (GMT) и координаты центра сцены в десятичных долях с соответствующим знаком для широты и долготы. Время извлекается из файла метаданных (строка SCENE_CENTER_TIME):

SCENE_CENTER_TIME = 15:28:49.4315057Z

Координаты - указанным выше способом.


Далее следуют важные параметры, требующие от вас некоторой смекалки. Это атмосферная модель, выбираемая из семи стандартных вариантов:

0	Поглощение отсутствует
1	Тропики
2	Средние широты (лето)
3	Средние широты (зима)
4	Субарктические широты (лето)
5	Субарктические широты (зима)
6	Стандарт US 62

Существует возможность задать собственную модель на основе, к примеру, измерений радиозондом.

Следующая строка - аэрозольная модель. Её также можно выбрать из семи стандартных вариантов:

0	Нет аэрозолей
1	Континентальная модель	 
2	Приморская модель
3	Модель для урбанизированных территорий
4	Модель для пустыни (Shettle)	 
5	Горение биомассы (biomass burning)
6	Стратосферная модель	 

Здесь тоже возможно задать собственную модель на основе различных параметров.

Понятно, что в подавляющем большинстве случаев придётся выбирать из стандартных и атмосферную, и аэрозольную модели.

Следующий параметр - Aerosol concentration model (visibility), модель концентрации аэрозолей. Если у вас есть измеренный метеорологическими методами показатель оптической видимости, достаточно указать его (одно число). Если его нет (что вероятно), следует указать в качестве параметра видимости 0, а в следующей строке записать параметр оптической толщины (optical depth), добытый ранее на порталах NASA.

Далее записывается средняя высота исследуемой местности над уровнем моря в километрах (с обратным знаком), извлеченная из ЦМР ранее.

Следующий параметр указывает положение съемочной аппаратуры. В случае спутниковых съемок здесь всегда будет -1000. 0 соответствует наземным наблюдениям.

Последний параметр - индекс используемого канала. Для основных каналов landsat 8 (OLI) используются следующие индексы:

B2 - 116
B3 - 117
B4 - 118
B5 - 120
B6 - 122
B7 - 123
B8 - 119 (Панхроматический)

Подробно обо всех этих параметрах можно прочесть в документации i.atcorr.


Теперь можно запускать, собственно, модуль i.atcorr. Для этого снова воспользуемся вкладкой Command Console:


Запуск модуля i.atcorr


В поле указания входного растрового слоя будем поочередно записывать откалиброванные ранее снимки (рекомендуется использовать reflectance). Указываем путь до файла atcorr_pars.txt, нажимаем кнопку load, указываем соответствующий индекс канала (загруженный файл конфигурации можно редактировать прямо в окне i.atcorr). Это необходимо проделать со всеми каналами, которые в дальнейшем мы собираемся использовать для построения индексов и проведения классификаций (в рассматриваемом примере: со 2 по 7).


Модуль i.atcorr


Во вкладке input необходимо обратить внимание на параметр input range. По умолчанию он установлен в интервале от 0 до 255. Это минимальное и максимальное значения оптической плотности по всему растру (в нашем случае значение reflectance). Дело в том, что алгоритм 6S принимает в качестве входного параметра ограниченный диапазон значений, который в дальнейшем будет приведен к значениям [0..1]. Укажем его, используя метаданные соответствующего канала. Для этого в дереве слоев в контекстном меню слоя выбираем Metadata, вызывая модуль r.info.


Вызов метаданных через контекстное меню слоя


Находим строку Range of data


Метаданные


Обычно этот диапазон (для reflectance) близок к единице. Обратим внимание, что в случае, когда максимальное значение меньше 1, GRASS автоматически приведет входной интервал i.atcorr к виду "от 0 до 255". В нашем случае подобная ситуация происходит сразу с несколькими каналами. Опытным путем было установлено, что модуль осуществляет автоматические отбрасывание знаков после запятой, поэтому мы имеем полное право ввести сразу округленные целочисленные границы. Также не забываем отметить флаг атрибута (r), Input raster map converted to reflectance (без флага, если корректируем radiance).


Настройки модуля i.atcorr


Стоит отметить, что на выходе получаются значения, приведенные к диапазону, указанному в atcorr output.


Настройки модуля i.atcorr - вкладка output


Поэтому, если сравнить результаты до и после коррекции, они будут значительно отличаться. Это необходимо учитывать и каждый канал приводить к одинаковым диапазонам. Опытным путём было установлено, что результаты расчетов индексов и других производных не зависят от выбранного диапазона.

Аналогично следует поступить со всеми остальными каналами, проявив внимательность при указании диапазонов и создании конфигурационного файла 6S. Атмосферную коррекцию можно считать законченной.


Расчет индекса NDVI

О том, что такое и зачем нужен NDVI замечательно написано в этой статье.

Существуют разные взгляды на то, как правильнее строить вегетационные индексы: по radiance или по reflectance (главное не в DN!). Вы можете поэкспериментировать. В примере мы воспользуемся полученными на прошлом шаге скорректированными за атмосферу данными reflectance.

Для расчета индекса снова откроем растровый калькулятор r.mapcalc, в котором, для удобства пользуясь разделами 'insert existing raster map' и 'operators' вводим выражение, не забыв указать имя выходного слоя (в нашем случае ndvi_refl).


Расчет NDVI в r.mapcalc


В сущности, результат получен. Поскольку индексное изображение одноканально, по умолчанию GRASS отрисует его в оттенках серого. Для того, чтобы раскрасить ndvi, воспользуемся модулем r.colors (запуск через Command Console, либо через контекстное меню слоя, пункт Set Color Table). Убедитесь, что в основной вкладке указаны нужные наборы данных, т.е. рассчитанные NDVI.


Окно модуля r.color


Собственно цветовая таблица задается на вкладке Define. Поскольку вегетационный индекс NDVI является классическим продуктом, в GRASS есть заранее заготовленная таблица. Так и называется - ndvi.


Выбор цветовой таблицы NDVI


В результате наше посчитанное по reflectance индексное изображение выглядит следующим образом:


Визуализированный с цветовой таблицей NDVI


Сравнение результатов с расчетами в ENVI после FLAASH

Для проверки результатов мы посчитали индекс NDVI также в среде ENVI 5.1, предварительно применив к набору исходных данных алгоритм атмосферной коррекции FLAASH. В GRASS существует несколько возможностей сравнить растровые изображения: начиная от простого визуального или попиксельного сравнения и заканчивая расчетом корреляции через модуль d.correlate. Мы воспользуемся модулем r.what, который позволяет простым способом осуществить попиксельное сравнение. Запустить его, помимо консольной команды, можно через меню Raster - Query Raster Maps - Query values by coordinates.


Запуск модуля r.what


Выбираем слои для сравнения (соответственно NDVI, рассчитанный в GRASS по reflectance, и в ENVI).


Выбор слоёв в модуле r.what


Результат модуль выдает в консольном виде. Мы обойдемся визуальной оценкой, хотя понятно, что по двум рядам значений вполне можно считать статистику. Видно, что в целом результаты очень близкие (с разницей в 1-2 сотых), хотя присутствуют явные "вылеты", когда разница достигает 3 десятых. Общее количество таких пикселей очень мало.


Попиксельное сравнение индексов после 6S и FLAASH


Pan Sharpening

Pan Sharpening (сокращение от Panchromatic sharpening) - процедура увеличения разрешения стандартных каналов (в случае Landsat 8: 30 метров) по данным панхроматического канала (15 метров). Обычно используется для получения более детального изображения, синтезированного в естественных цветах. Следует отметить, что при создании "улучшенного" изображения отдельно взятого канала происходит, по большому счету, изобретение информации, которой в исходном наборе данных нет. Производить классификацию или какой-либо спектральный анализ по таким изображениям бессмысленно. Используются они для улучшения качества визуального дешифрирования, более точного создания обучающих выборок для классификаций, а также для иллюстрирования различных материалов.

Для работы с панхроматическим снимком необходимо, в первую очередь, изменить region нашего рабочего пространства, так как в нём устанавливается не только пространственный охват, но и размеры изображения (для 8го канала оно вдвое выше, чем у остальных стандартных). Для этого в уже знакомом меню Region - Set Region нужно выбрать новый источник данных. Либо просто воспользоваться командой:

g.region raster=LC80110462013131LGN01_B8

Будем работать с исходными файлами (без калибровки и коррекции) по вышеозвученным причинам. После того, как мы сменили регион, нужно изменить диапазон значений у всех предполагаемых к использованию файлов, т.е. у стандартных каналов (например, для синтеза в естественных цветах, B2,B3,B4) и у B8 на единый диапазон [0;255]. Это связано с особенностями работы модуля i.pansharpen, с помощью которого мы и будем осуществлять увеличение детализации. Согласно документации, в качестве исходных данных для него предполагается использование только 8-битных растров. Мы поэкспериментировали и убедились в том, что без предварительной процедуры изменения диапазона действительно ничего не выйдет. Поэтому запускаем модуль r.rescale через меню Raster - Change category values and labels - rescale и на вкладке Reqiured указываем преобразуемые слои и диапазон.


Модуль r.rescale


Далее запускаем модуль i.pansharpen (меню Imagery - Pan sharpening). Во вкладке Reqiured указываем растры, которые будут использованы как красный, синий и зелёный каналы для результирующего изображения (в нашем случае B4, B3, B2), имя для выходных растров (их будет три), а также pan sharpening метод. Во вкладке Optional можно установить флаг Rebalance, в нашем случае это позволило получить более естественные цвета растительности. После нескольких экспериментов более приемлемым нам показался результат, полученный с помощью метода ihs.


Модуль i.pansharpen


После выполнения команды будут созданы три одноканальных изображения, имеющих имена с ключами _red, _blue, _green. Теперь нужно их объединить для синтезирования детального изображения в естественных цветах с помощью модуля r.composite (Raster - Manage colors - Create RGB), в котором достаточно указать растры соответствующих каналов и имя для выходного слоя.


Модуль r.composite


Создание детализированного изображение окончено. Также, для визуальной оценки эффективности этой и последующих процедур, создадим синтезированное изображение в естественных цветах для исходных, 30-метровых каналов, применив модуль r.composite для соответствующих растровых слоёв (B2, B3, B4).


Исходное изображение. Синтез в естественных цветах


Сравним фрагмент исходного RGB с результатом процедуры Pan Sharpening.


Pan Sharpening. Исходное (слева) и результирующее (справа) изображения


Создание масок воды и облаков

Процедура выделения воды и облаков может показаться ненужной или избыточной, однако подобная задача возникает часто. Такие маски помогают проводить процедуру классификации, оценивать площади недоступных из-за облаков территорий и территорий водоёмов. Существуют различные подходы к решению этой задачи, мы познакомим вас с одним из них. Для выделения воды рассчитаем индекс Water Ratio Index (WRI) с помощью уже знакомого нам растрового калькулятора. Индекс WRI рассчитывается по формуле:

(GREEN + RED) / (NIR + SWIR)

Воспользуемся растровым калькулятором r.mapcalc, а в качестве исходных данных используем откорректированные ранее методом 6S reflectance-изображения.


Расчет индекса WRI


Оценив полученный растр, выбираем диапазон значений, в который попадает вода. Индекс WRI очень явно выделяет её, в нашем случае это были значения больше либо равные 1.0. Для выделения пикселей только со значением большим или равным 1.0 снова воспользуемся калькулятором растров.


Выделение пикселей со значениями >= 1.0


Также выделим в воду значения No Data с любого изображения, полученного после атмосферной коррекции (затем, при создании карты, мы почистим полученные векторные маски). В No Data после атмосферной коррекции обратились все пиксели, которые не попали на сушу и, соответственно, не имели значения рельефа из SRTM. Для снимков, целиком охваченных ЦМР, такой проблемы не возникнет. Осуществляем задуманное снова с помощью калькулятора растров.


Выделение пикселей со значением No Data


Для создания маски облаков воспользуемся изображением с ключом _BQA (в нашем случае LC80110462013131LGN01_BQA). Найдем в таблице стандартных значений пикселей для Landsat 8 нужные значения для Cloud и Cirrus, которые встречаются на нашем снимке, и создадим маску облаков через растровый калькулятор. Выражение выглядит следующим образом:


Выделение пикселей со стандартными значениями из документации Landsat 8 для _BQA


Мы получили маски воды и облаков отдельно друг от друга и в растровом виде. Вода нам пригодится векторная, поэтому следует конвертировать её. Затем объединим полученные слои и зададим результат в качестве векторной маски.

Для преобразования растра в вектор воспользуемся модулем r.to.vect (Raster - Map type conversions - Raster to vector), задав соответствующие параметры. Проделывая эту процедуру дважды, получаем слои water_mask и cloud_mask.


Модуль r.to.vect


Теперь объединяем эти векторные слои с помощью модуля v.patch (Vector - Overlay Vector maps - Patch Vector maps) в слой mask.


Модуль v.patch


Полученный слой:


Векторная маска воды и облаков

Если предполагается сразу перейти к процедуре классификации, следует установить полученный слой как маску через модуль r.mask (Raster - Mask). Нужно выбрать вкладку Vector, указать слой mask в качестве "Name of vector map to use as mask", а также во вкладке Create установить флаг Create inverse mask, иначе нам останется доступным только то, что попадет под маску, хотя мы добиваемся обратного эффекта.


Модуль r.mask


Классификация без обучения

Переходим к задаче классификации, начиная с неконтролируемого метода как наиболее простого. Во-первых, необходимо создать группу и подгруппу снимков, по которым будет проводиться классификация. Запускаем модуль i.group (Imagery - Develop images and groups - create/edit group). Создаем группу atcorr_group, включая туда все растры, которые мы получили после атмосферной коррекции, а также любые другие, какие вам могут пригодиться в дальнейшей работе с группами. Нажимаем Apply. Теперь создаем подгруппу, где выбираем (в нашем случае) все корректированные каналы со 2 по 7: они будут использоваться как факторы классификации.


Модуль i.group


Сама неконтролируемая классификация выполняется в два этапа. Сначала необходимо сгенерировать сигнатуры с помощью модуля i. cluster (Imagery - classify image - clustering input for unsupervised classification). Задаем необходимые параметры на вкладках Required и Settings. Здесь все зависит от характера ваших данных, многие параметры подбираются только через многочисленные эксперименты. Даже сами создатели модуля говорят о существовании “magic” numbers, которые подбираются чисто эмпирически (хотя, безусловно, они имеют конкретный алгоритмический смысл). Основной параметр, конечно, это количество классов (Initial number of classes).


Модуль i.cluster. Вкладка Required
Модуль i.cluster. Вкладка Settings


Далее запускаем модуль i.maxlik (Imagery - Classify image - Maximum likelihood classification (MLC)). Настраивать здесь ничего не нужно, просто указываем имя файла сигнатур, сгенерированного на предыдущем шаге.


Модуль i.maxlik


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


Результат классификации без обучения


Классификация с обучением

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


Итак, для проведения классификации с обучением нам понадобятся обучающие полигоны. Их можно как создать в GRASS, так и импортировать уже существующие. Главное, чтобы тип данных поля, отвечающего за разделение полигонов на классы, был числовым. Импортируются векторные данные аналогично растровым. Рассмотрим создание собственных полигонов в среде GRASS. Создадим векторный файл с помощью меню Vector - Develop Vector maps - Create new vector map. В открывшемся окне введем название слоя, атрибут менять не будем, для нашей задачи подойдут категории.


Создание векторного слоя


Созданный слой добавился в дерево слоев. Для удобства создания обучающих полигонов в подгрузим полученный ранее растр с результатами Pan Sharpening. Откроем режим Digitize в окне отображения и слева в маленьком окне выберем наш слой.


Режим Digitize
Инструмент создания объектов и кнопка настроек


До начала векторизации зайдем в настройки. Нам нужно, чтобы полигоны замыкались по окончанию (флаг close boundary во вкладке General), также нужно настроить присваивание категории созданному объекту (Category mode - manual entry во вкладке Attributes).


Настройки векторизации. Вкладка General
настройки векторизации. Вкладка Attributes


При таком подходе удобно векторизовать по классам: устанавливаем в поле category number (вкладка Attributes) значение категории '1' и создадим столько, сколько нужно, полигонов, относящихся к данному классу. Затем изменяем category number на 2, векторизуем полигоны второго класса - и так далее. С количеством классов необходимо определиться заранее, руководствуясь объемами имеющейся вспомогательной информации и общим характером местности. Для примера мы выбрали всего пять общих категорий: Урбанизированные территории, растительность, сельскохозяйственные угодья, земли без растительного покрова, не классифицируемые без вспомогательных данных территории. По окончанию нажмем кнопку выхода из режима векторизации, согласившись сохранить изменения, и получаем векторный обучающий слой.

Модуль классификации работает с обучающими полигонами в растровом формате. Поэтому нам нужно совершить преобразование с помощью модуля v.to.rast (Vector - Map type conversions - Vector to raster). Никаких особенных настроек он не предполагает:


Модуль v.to.rast


Далее запускаем модуль i.gensigset (Imagery - Classify Image - Input for supervised SMAP), который позволит сгенерировать сигнатуры для классификации с обучением. Особое внимание обратите на параметр 'maximum number of sub-signtures in any class' во вкладке Optional, изменение которого существенно изменяет результат. Первоначальную классификацию можно провести с настройками по умолчанию, задав только опции во вкладке Required: выбрать обучающий слой, группу и подгруппу наборов классифицируемых данных, имя выходного файла. После получения первичного результата можно будет начать эксперименты с другими настройками.


Модуль i.gensigset


Теперь запускаем модуль i.smap (Imagery - Classify image - Sequential maximum a posteriori classification (SMAP)). Задаем знакомые параметры во вкладке General. Во вкладке Optional можно установить флаг 'Use maximum likelihood estimation', возможно, для ваших данных, это будет решающий момент (в любом случае поэкспериментируйте с флагом и без него).


Модуль i.smap


Как только мы получаем, путем всевозможных комбинаций доступных параметров (возможно, вы поймете что следует модифицировать обучающие выборки), удовлетворительный результат - переходим к постклассификации.


Классифицированное изображение


Постклассификация

Постклассификация - это субъективный, но важный этап. Мы воспользуемся способом скользящего окна для того, чтобы избавиться от всевозможных шумов (отдельных пикселей) и генерализовать полученный результат. Запускаем модуль r.neighbors (Raster - Neighborhood analysis - Moving window). Заполняем обязательные поля во вкладке Required (растр для фильтрации, имя для результирующего растра) и переходим к выбору фильтра (вкладка Neighborhood). В нашем случае, когда каждый пиксель имеет всего 5-7 (может быть меньше или больше, в зависимости от количества классов) возможных значений, определенных однозначно, и не может быть никаких промежуточных вариантов, больше всего подходит метод mode, который присваивает пикселю то значение, которое чаще всего встречается среди его соседей. Результат сильно зависит от размера фильтра (кстати, значение может быть только нечетным, иначе будет ошибка). Для примера мы остановились на размере 7х7 (Neighborhood size = 7).


Модуль r.neighbors. Вкладка Required
Модуль r.neighbors. Вкладка Neighborhood


В результате получаем заметно отфильтрованное изображение:


Отфильтрованное изображение


Также, если вы хотите сравнить результаты нескольких методов, то в Name of output можно через запятую перечислить несколько имен (по количеству используемых методов), а во вкладке Neighborhood установить несколько флагов (названия присваиваются соответственно по порядку методов). Конечный результат можно, при необходимости, уже знакомыми вам методами, перевести в вектор. Подгрузив в отображение векторный слой воды мы получим практически готовую карту.

Карта

В сущности, задача решена, остается только оформить его (если это вообще требуется). GRASS не является идеальной средой для оформления макетов карт, поэтому мы использовали для этих целей открытую ГИС Quantum GIS. Подразумевается, что мы произвели предварительное преобразование результатов постклассификации в векторный формат с помощью модуля r.to.vect. Далее, из GRASS можно экспортировать (с помощью модуля v.out.ogr (File - Export Vector map - Common Export formats)) этот векторный слой конечного результата постклассификации (при экспорте нежелательно устанавливать флаг напротив 'Also export features without category…'. во вкладке Selection, только если вы знаете, что полигоны без категорий вам тоже необходимы). Если вы экспортируете результаты в ESRI SHP, то также не забудьте выгрузить и файл .prj, установив соответствующий флаг во вкладке Creation. Дополнительно мы экспортировали в .shp слой воды, который создавался на этапе маскирования.


Модуль v.out.ogr. Вкладка Required
Модуль v.out.ogr. Вкладка Selection
Модуль v.out.ogr. Вкладка Creation


Экспортированные данные оформляются в удобной для вас среде. В нашем случае была получена карта такого вида:


Файл:Grass7 landsat8 processing 74.jpg
Оформленная в QGIS карта