Обработка и интерпретация данных Landsat 8 (OLI) средствами GRASS GIS 7
по адресу http://gis-lab.info/qa/grass7-landsat8-processing.html
Использование открытой геоинформационной системы 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. Конечная редакция скриншотов и текстов команд выполнялась в Windows 8 x64.
- Все (или почти все) действия в статье выполняются через графический интерфейс GRASS, несмотря на то, что часто проще и быстрее использовать командную строку. Это сделано для увеличения наглядности и уменьшения порога вхождения. Основные операции сопровождаются текстом соответствующей команды.
- Предполагается, что вы загрузили и распаковали набор данных Landsat 8 уровня обработки "Level 1 Product" (Например с помощью EarthExplorer, краткие инструкции можно получить здесь или здесь) с одиннадцатью одноканальными изображениями и файлом метаданных.
- Не гарантируется корректная работа модулей 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 и d.rast:
d.mon wx0 d.rast LC80110462013131LGN01_B1
Далее необходимо задать регион с помощью модуля g.region. Для этого обращаемся в меню Settings - Region - Set Region.
Поскольку все изображения одного набора данных Landsat точно соответствуют друг другу геометрически, в качестве региона можно выбрать любой из стандартных каналов (один) в поле [multiple] Set Region to match raster map(s) ("установить регион в соответствии с растровой картой"). Оговоримся, что на этом этапе следует выбрать канал с разрешением 30 метров (т.е. любой, кроме 8-го, панхроматического).
g.region raster=LC80110462013131LGN01_B2
Шаг с выбором региона отдельной функцией можно расценивать как избыточный, так как на этапе импорта данных в r.in.gdal можно установить флаг 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.
Во вкладке required ("обязательные") задаем basename ("базовое имя") для исходных и выходных снимков. Для исходных это строка, которая должна полностью совпадать с именами загруженных вами снимков, только без цифры канала, которую модуль будет подставлять автоматически. В рассматриваемом примере: LC80110462013131LGN01_B. Для выходных — произвольная понятная вам строка. Разумеется, желательно, чтобы по ней был понятен характер хранимой информации. Например, в нашем примере это toar_refl_B (результирующие растры будут иметь имена toar_refl_B1, toar_refl_B2 и т.д.).
Во вкладке metadata достаточно указать путь к файлу метаданных (в примере LC80110462013131LGN01_MTL.txt). Параметры сенсора и прочие параметры модуль автоматически прочитает из него.
Во вкладке optional ключевой флаг отвечает за тип результатов. По умолчанию модуль рассчитывает отражательную способность (reflectance), для расчета же энергетической яркости (radiance) нужно отметить флаг Output at-sensor radiance isntead of reflectance for all bands ("рассчитать яркость вместо отражательной способности для всех каналов").
Примечание: если у вас возникает необходимость перезаписать какие-либо из уже существующих файлов, практически в каждом модуле (во вкладке Optional или где-нибудь еще) доступен флаг "Allow output files overwrite existing files" (в командной строке за это отвечает опция "--overwrite" или её короткий вариант "--o"). Ещё примечание: флаг verbose module output (подробный вывод), также доступный во вкладке Optional, позволяет получать детализированный отчет о ходе выполнения операции. В командной строке этой цели служит опция "--verbose" или её короткий вариант "--v".
i.landsat.toar --overwrite --verbose input=LC80110462013131LGN01_B output=toar_refl_B metfile=G:\data\image_proc\Cuba\input\LC80110462013131LGN01_MTL.txt
После запуска программа будет выполнять последовательное преобразование всех растров до последнего существующего в рабочей области канала (в нашем случае включительно до 8-го панхроматического) .
Если вы не импортировали какие-то каналы на предыдущих этапах, то модуль сообщит об ошибке: "Unable to open header file for raster map..." Это не должно вас напугать, всё что мы добавили — благополучно обработано.
Обратим также внимание на параметр Esun. В документации к i.landsat.toar указано, что результаты без параметра Esun дают приблизительно такие же значения, как и с ним, и сообщение вида "WARNING: ESUN evaluated from REFLECTANCE_MAXIMUM_BAND", как и вышеупомянутая ошибка, не стоит ваших переживаний.
Сейчас мы произвели расчет отражательной способности (reflectance). Для расчета спектральной яркости (radiance) достаточно переименовать выходные слои в поле "Prefix for output raster maps" (toar_radiance_B) и на вкладке Optional поставить выше упоминавшийся флаг Output at-sensor radiance instead of reflectance for all bands (-r). Запускаем модуль еще раз, получаем искомое.
i.landsat.toar -r --overwrite --verbose input=LC80110462013131LGN01_B output=toar_radiance_B metfile=G:\data\image_proc\Cuba\input\LC80110462013131LGN01_MTL.txt
Теперь обещанные развлечения с калькулятором растров. Согласно документации, radiance рассчитывается по следующей формуле:
Lλ = ((Lmax - Lmin)/(Qcalmax - Qcalmin))*(Qcal - Qcalmin) + Lmin
где Lλ — количество приходящего излучения; Lmin — количество приходящего излучения, которое после масштабирования становится Qmin; Lmax — количество приходящего излучения, которое после масштабирования становится Qmax; Qcalmin — минимальное калиброванное значение DN (0 или 1); Qcalmax — максимальное калиброванное значение DN; 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.
r.mapcalc expression='mapc_radian_B2=(762.78229-(-62.99081))/(65535-1)*(LC80110462013131LGN01_B2-1)+(-62.99081)'
Больше ничего настраивать не нужно. Запускаем и получаем новый растр. Добавим 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). Импортируем их как обычные растровые слои. В примере использованы наборы данных srtm_21_08 и srtm_21_09.
Оценку состояния атмосферы можно посмотреть здесь и здесь. В частности, нам нужен параметр оптической толщины аэрозолей для диапазона 550 нм (только число, наборы данных не нужны). В нашем случае значение параметра 0.120.
В качестве одного из параметров необходимо будет указать среднюю высоту (помимо самой ЦМР). Рассчитаем её с помощью модуля r.univar. Запускаем его через меню Raster - Reports and statistics - Univariate raster statistics или через командную строку во вкладке Command Console:
В появившемся окне выбираем слой с рельефом (ЦМР).
r.univar map=Cuba_SRTM
Запускаем модуль и отслеживаем результат работы во вкладке Command output. Запоминаем (записываем) значение средней высоты (mean).
Также нам понадобятся координаты центра снимка. Если набор исходных данных был загружен с сайта EarthExplorer, то при выборе снимка указывались географические координаты центра сцены. Можете быстро найти свою сцену по имени.
Создаем файл с параметрами для 6S (например atcorr.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. Для этого воспользуемся меню (Imagery - Satellite images tools - Atmospheric correction) или вкладкой Command Console:
В поле указания входного растрового слоя будем поочередно записывать откалиброванные ранее снимки (рекомендуется использовать reflectance). Указываем путь к файлу atcorr.txt, нажимаем кнопку "Load", указываем соответствующий индекс канала (загруженный файл конфигурации можно редактировать прямо в окне i.atcorr). Это необходимо проделать со всеми каналами, которые в дальнейшем мы собираемся использовать для построения индексов и проведения классификаций (в рассматриваемом примере со 2 по 7).
Во вкладке Input необходимо обратить внимание на параметр "Input range". По умолчанию он установлен в интервале от 0 до 255. Это минимальное и максимальное значения оптической плотности по всему растру (в нашем случае значение reflectance). Дело в том, что алгоритм 6S принимает в качестве входного параметра ограниченный диапазон значений, который в дальнейшем будет приведен к значениям [0..1]. Укажем его, используя метаданные соответствующего канала. Для этого в списке слоев в контекстном меню слоя выбираем Metadata, вызывая модуль r.info.
r.info map=toar_refl_B2
Находим строку "Range of data".
Обычно диапазон значений для reflectance близок к единице. Обратим внимание, что в случае, когда максимальное значение меньше 1, GRASS автоматически приведет входной интервал i.atcorr к виду "от 0 до 255".
В нашем случае подобная ситуация происходит сразу с несколькими каналами.
Опытным путем было установлено, что модуль осуществляет автоматические отбрасывание знаков после запятой, поэтому мы имеем полное право ввести сразу округленные целочисленные границы. Также не забываем отметить флаг атрибута (-r), Input raster map converted to reflectance ("исходные данные конвертированы в отражательную способность"). Соответственно, запускаем модуль без этого флага, если корректируем radiance.
i.atcorr -r --overwrite --verbose input=toar_refl_B2 range=0,1 elevation=Cuba_SRTM parameters=G:\data\image_proc\Cuba\atcorr.txt output=atcorr_refl_B2 rescale=0,255
Стоит отметить, что на выходе получаются значения, приведенные к диапазону, указанному в atcorr output.
Если сравнить результаты до и после коррекции, они будут значительно отличаться. Это необходимо учитывать, и каждый канал приводить к одинаковым диапазонам.
Опытным путём было установлено, что результаты расчетов индексов и других производных не зависят от выбранного диапазона.
Аналогично следует поступить со всеми остальными каналами, проявив внимательность при указании диапазонов и создании конфигурационного файла 6S. Атмосферную коррекцию можно считать законченной.
Расчет индекса NDVI
О том, что такое NDVI и зачем нужен этот индекс, замечательно написано в этой статье.
Существуют разные взгляды на то, как правильнее строить вегетационные индексы: по radiance или по reflectance (главное не в DN!). Вы можете поэкспериментировать. В примере мы воспользуемся полученными на прошлом шаге скорректированными с учётом атмосферы данными reflectance.
Для расчета индекса снова откроем растровый калькулятор r.mapcalc, в котором, для удобства используем разделы 'insert existing raster map' и 'operators' и вводим выражение, не забыв указать имя выходного слоя (в нашем случае "ndvi_refl").
r.mapcalc expression='ndvi_refl=(atcorr_refl_B5 - atcorr_refl_B4)/(atcorr_refl_B5 + atcorr_refl_B4)'
В сущности, результат получен. Для того, чтобы раскрасить растр NDVI в принятую для этого индекса цветовую схему, воспользуемся модулем r.colors (запуск через Command Console, либо через контекстное меню слоя, пункт Set Color Table (Задать цветовую таблицу)). Убедитесь, что в основной вкладке указаны нужные наборы данных, т.е. рассчитанные NDVI.
Собственно цветовая шкала задается на вкладке Define. Поскольку индекс NDVI является классическим продуктом, в GRASS есть заранее заготовленная шкала ("ndvi").
r.colors ndvi_refl,ndvi_FLAASH color=ndvi
В результате наше посчитанное по reflectance индексное изображение выглядит следующим образом:
Сравнение результатов с расчетами в ENVI после FLAASH
Для проверки результатов мы посчитали индекс NDVI также в среде ENVI 5.1, предварительно применив к набору исходных данных алгоритм атмосферной коррекции FLAASH. В GRASS существует несколько возможностей сравнить растровые изображения, начиная от простого визуального или попиксельного сравнения и заканчивая расчетом корреляции через модуль d.correlate. Мы воспользуемся модулем r.what, который позволяет простым способом осуществить попиксельное сравнение. Запустить его, помимо консольной команды, можно через меню: Raster - Query Raster Maps - Query values by coordinates.
Выбираем слои для сравнения (соответственно NDVI, рассчитанный в GRASS по reflectance, и в ENVI).
r.what map=ndvi_refl,ndvi_FLAASH points=Cuba_check_pts
Результат выдается в консольном виде. Мы обойдемся визуальной оценкой, хотя понятно, что по двум рядам значений вполне можно посчитать статистику. Видно, что в целом результаты очень близкие (третья колонка - данные после FLAASH, четвертая - после ATCORR).
Паншарпенинг
Паншарпенинг (pan-sharpening, сокращение от "panchromatic sharpening") — процедура увеличения разрешения стандартных каналов (для Landsat 8 это каналы с разрешением 30 метров) по данным панхроматического канала (15 метров). Обычно используется для получения более детального изображения, синтезированного в естественных цветах. При работе с результатами паншарпенинга следует учитывать, что процедура существенно изменяет свойства изображения с точки зрения спектральных характеристик снимка, что, в свою очередь, ощутимо повлияет на результаты классификации. Обычно алгоритм используются для улучшения качества визуального дешифрирования, более точного создания обучающих выборок для классификаций, а также для иллюстрирования различных материалов.
Для начала, чтобы получить возможность визуально оценить эффективность этой и следующих процедур, создадим синтезированное изображение в естественных цветах для исходных 30-метровых каналов. Чтобы обеспечить более выразительную цветовую гамму, можно перед синтезом воспользоваться модулем i.colors.enhance, выполняющим автоматическую балансировку цветов для последующего синтеза в естественных цветах. При работе с ним можно поэкспериментировать с флагами и опциями (особое внимание следует уделить параметру strength (рекомендуемые значения от 90 до 99)). Модуль i.colors.enhance не создает новых наборов данных, модифицируя исходные, заданные в параметрах red, green и blue.
i.colors.enhance red=LC80110462013131LGN01_B4 green=LC80110462013131LGN01_B3 blue=LC80110462013131LGN01_B2 strength=98
Собственно синтез осуществляется с применением модуля r.composite (Raster - Manage colors - Create RGB), в котором достаточно указать растры соответствующих каналов (B2, B3, B4) и имя выходного слоя.
r.composite red=LC80110462013131LGN01_B4 green=LC80110462013131LGN01_B3 blue=LC80110462013131LGN01_B2 output=rgb
В результате синтеза исходных каналов, сбалансированных с помощью i.color.enhance, получается такое изображение:
Теперь приступим к работе с панхроматическим снимком. Для этого необходимо, в первую очередь, изменить регион нашего рабочего пространства, так как в нём устанавливается не только пространственный охват, но и размеры изображения (для 8-го канала оно вдвое больше, чем у остальных стандартных). Для этого с помощью уже знакомого модуля g.region (в меню Region - Set Region) нужно выбрать новый источник данных (B8). Либо просто воспользоваться командой:
g.region raster=LC80110462013131LGN01_B8
Будем работать с исходными файлами (автоматически сбалансированными, но без калибровки и коррекции). После того, как мы сменили регион, нужно изменить диапазон значений у всех предполагаемых к использованию файлов, т.е. у стандартных каналов (например, для синтеза в естественных цветах — B2,B3,B4), а также у канала B8, на единый диапазон [0;255]. Это связано с особенностями работы модуля i.pansharpen, с помощью которого мы и будем осуществлять увеличение детализации. Согласно документации, в качестве исходных данных для него предполагается использование только 8-битных растров. Мы поэкспериментировали и убедились в том, что без предварительной процедуры изменения диапазона действительно ничего не выйдет. Поэтому запускаем модуль r.rescale через меню Raster - Change category values and labels - rescale и на вкладке Required указываем преобразуемые слои и диапазон.
r.rescale input=LC80110462013131LGN01_B2 output=255_B2 to=0,255
Далее запускаем модуль i.pansharpen (в меню Imagery - Pan sharpening). На вкладке Required указываем растры, которые будут использованы как красный, синий и зелёный каналы для результирующего изображения (в нашем случае B4, B3, B2), имя для выходных растров (их будет три), а также метод для паншарненинга. На вкладке Optional можно установить флаг Rebalance, что может дать интересные результаты, но в примере обойдемся без него. После нескольких экспериментов более приемлемым нам показался результат, полученный с помощью метода IHS.
i.pansharpen red=255_B4 green=255_B3 blue=255_B2 pan=255_B8 output=pansharp method=ihs
После выполнения команды будут созданы три одноканальных изображения, имеющих имена с ключами _red, _blue, _green. Теперь нужно их объединить для синтезирования детального изображения в естественных цветах с помощью уже знакомого модуля r.composite.
r.composite red=pansharp_red green=pansharp_green blue=pansharp_blue output=panrgb
Создание детализированного изображение окончено.
Сравним фрагмент исходного сбалансированного RGB с результатом процедуры паншарпенинга.
Создание масок воды и облаков
Процедура выделения воды и облаков может показаться ненужной или избыточной, однако подобная задача возникает часто. Такие маски помогают проводить процедуру классификации, оценивать площади недоступных из-за облаков территорий и территорий водоёмов. Существуют различные подходы к решению этой задачи, мы познакомимся с одним из них. Для выделения воды рассчитаем индекс "Water Ratio Index" (WRI) с помощью уже знакомого нам растрового калькулятора. Индекс WRI рассчитывается по формуле:
(GREEN + RED) / (NIR + SWIR)
Воспользуемся растровым калькулятором "r.mapcalc", а в качестве исходных данных используем откорректированные ранее методом 6S растры со значениями reflectance.
r.mapcalc expression=wri=(atcorr_refl_B3 + atcorr_refl_B4)/(atcorr_refl_B5 + atcorr_refl_B6)
Оценив полученный растр, выбираем диапазон значений, в который попадает вода. Индекс WRI очень явно выделяет её, в нашем случае это значения больше либо равные 1.0. Для выделения пикселей только со значением большим или равным 1.0 снова воспользуемся калькулятором растров.
r.mapcalc expression='wri_water=if(wri>=1.0,1,null())'
Также выделим в воду значения "No Data" с любого изображения, полученного после атмосферной коррекции (затем, при создании карты, мы почистим полученные векторные маски). В "No Data" после атмосферной коррекции обратились все пиксели, которые не попали на сушу и, соответственно, не имели значения рельефа из SRTM. Для снимков, целиком охваченных ЦМР, такой проблемы не возникнет. Осуществляем задуманное снова с помощью калькулятора растров.
r.mapcalc expression='wri_nodata=if(isnull(atcorr_refl_B2),1,null())'
Для создания маски облаков воспользуемся изображением с ключом _BQA (в нашем случае LC80110462013131LGN01_BQA). Найдем в таблице стандартных значений пикселей для Landsat 8 нужные значения для "Cloud" и "Cirrus", которые встречаются на нашем снимке, и создадим маску облаков через растровый калькулятор. Выражение выглядит следующим образом:
r.mapcalc expression='cloudmask=if( LC80110462013131LGN01_BQA==53248,1,null() ) ||| if( LC80110462013131LGN01_BQA==28672,1,null() ) ||| if( LC80110462013131LGN01_BQA==28704,1,null() )'
Примечание: при работе в командной строке в среде Windows оператор ||| может вызывать конфликт интерпретации. В таком случае вы можете воспользоваться ещё одним способом запуска команд в GRASS: создайте текстовый файл (например, "cloud.txt") и запишите в него команду:
expression=cloudmask=if( LC80110462013131LGN01_BQA==53248,1,null() )||| if( LC80110462013131LGN01_BQA==28672,1,null() )||| if( LC80110462013131LGN01_BQA==28704,1,null() )
Затем из командной строки GRASS выполните команду, указав вместо "cloud.txt" абсолютный путь к созданному текстовому файлу:
r.mapcalc file=cloud.txt
Такой подход часто бывает полезным.
Итак, мы получили маски воды и облаков в растровом виде и отдельно друг от друга. Маска воды нам пригодится в векторном виде, поэтому следует конвертировать её в вектор. Затем объединим полученные слои и зададим результат в качестве векторной маски.
Для преобразования растра в вектор воспользуемся модулем r.to.vect (Raster - Map type conversions - Raster to vector), задав соответствующие параметры. Проделывая эту процедуру дважды, получаем слои "water_mask" и "cloud_mask".
r.to.vect input=wri_water output=water_mask type=area
Теперь объединяем эти векторные слои с помощью модуля v.patch (Vector - Overlay Vector maps - Patch Vector maps) в слой mask.
v.patch input=water_mask,cloudmask,wri_nodata output=cw_mask
Полученный слой:
Если предполагается сразу перейти к процедуре классификации, следует установить полученный слой как маску через модуль r.mask ("Raster - Mask"). Нужно выбрать вкладку Vector, указать слой mask в графе "Name of vector map to use as mask", а также на вкладке Create установить флаг Create inverse mask ("инвертировать маску") (-i), иначе нам останется доступным только то, что попадет под маску, хотя мы добиваемся обратного эффекта.
r.mask -i vector=cw_mask
Классификация без обучения
Переходя к задаче классификации, начнём с метода неконтролируемой классификации как наиболее простого. Во-первых, необходимо создать группу и подгруппу снимков, по которым будет проводиться классификация. Запускаем модуль i.group (Imagery - Develop images and groups - Create/edit group). Создаем группу "atcorr_group", включая туда все растры, которые мы получили после атмосферной коррекции, а также любые другие, какие вам могут пригодиться в дальнейшей работе с группами. Нажимаем "Apply". Теперь создаем подгруппу, где выбираем (в нашем случае) все корректированные каналы со 2 по 7 — они будут использоваться как факторы классификации.
i.group group=atcorr_group subgroup=atcorr_subgroup input=atcorr_refl_B2,atcorr_refl_B3,atcorr_refl_B4,atcorr_refl_B5,atcorr_refl_B6,atcorr_refl_B7
Сама неконтролируемая классификация выполняется в два этапа. Сначала необходимо сгенерировать так называемые сигнатуры с помощью модуля i.cluster ("Imagery - classify image - clustering input for unsupervised classification"). Задаем необходимые параметры на вкладках Required и Settings. Здесь все зависит от характера данных, многие параметры подбираются только через многочисленные эксперименты. Даже сами создатели модуля говорят о существовании "магических чисел" ("magic numbers"), которые подбираются чисто эмпирически (хотя, безусловно, они имеют конкретный алгоритмический смысл). Основной параметр — количество классов ("Initial number of classes").
i.cluster group=atcorr_group subgroup=atcorr_subgroup signaturefile=cluster_sign classes=6
Далее запускаем модуль i.maxlik (Imagery - Classify image - Maximum likelihood classification (MLC)). Настраивать здесь ничего не нужно, просто указываем имя файла сигнатур, сгенерированного на предыдущем шаге.
i.maxlik group=atcorr_group subgroup=atcorr_subgroup signaturefile=cluster_sign output=clusters_6
В результате получается классифицированное изображение, степень интерпретируемости которого сильно зависит от характера местности и подобранных параметров.
Классификация с обучением
Классификация с обучением — более сложный метод, предполагающий создание обучающих выборок. Оговоримся, что качественная классификация предполагает широкое использование дополнительной информации: картографических материалов, данных полевых наблюдений, ландшафтных описаний и так далее. Даже создание базового набора обучающих выборок на совершенно незнакомую территорию вряд ли увенчается успехом. Кроме того, любые результаты перед использованием необходимо верифицировать. В примере демонстрируется техническая составляющая процесса, авторы не ставили перед собой задачи получения детального и высококачественного результата.
Итак, для проведения классификации с обучением нам понадобятся обучающие полигоны. Их можно как создать в GRASS, так и импортировать уже существующие. Главное, чтобы тип данных поля, отвечающего за разделение полигонов на классы, был числовым. Импортируются векторные данные аналогично растровым. Рассмотрим создание собственных полигонов в среде GRASS. Создадим векторный файл с помощью меню Vector - Develop Vector maps - Create new vector map. В открывшемся окне введем название слоя, атрибут менять не будем, для нашей задачи подойдут категории.
Созданный слой добавился в список слоев, отображаемых на экране.
Для удобства создания обучающих полигонов подгрузим полученный ранее растр с результатами паншарпенинга. Откроем режим Digitize в окне отображения и слева в маленьком окне выберем наш слой. Обратите внимание на наличие в этом окне пункта "New vector layer", использование которого является альтернативным способом создать векторный слой (откроется диалог, аналогичный вышеприведенному).
До начала векторизации зайдем в настройки. Нам нужно, чтобы полигоны замыкались по окончанию (флаг Close boundary на вкладке General), также нужно настроить присваивание категории созданному объекту (Category mode - manual entry на вкладке Attributes).
При таком подходе удобно векторизовать по классам: устанавливаем в поле "category number" (вкладка Attributes) значение категории '1' и создадим столько, сколько нужно, полигонов, относящихся к данному классу. Затем изменяем category number на 2, векторизуем полигоны второго класса, и так далее. С количеством классов необходимо определиться заранее, руководствуясь объемами имеющейся вспомогательной информации и общим характером местности. Для примера мы выбрали всего пять общих категорий: урбанизированные территории, растительность, сельскохозяйственные угодья, земли без растительного покрова, не классифицируемые без вспомогательных данных территории.
По окончанию нажмем кнопку выхода из режима векторизации, согласившись сохранить изменения, и получаем векторный слой.
Модуль классификации работает с обучающими полигонами в растровом формате. Поэтому нам нужно совершить преобразование с помощью модуля v.to.rast (Vector - Map type conversions - Vector to raster). Никаких особенных настроек он не предполагает:
v.to.rast input=classes output=classes use=cat
Далее запускаем модуль i.gensigset (Imagery - Classify Image - Input for supervised SMAP), который позволит сгенерировать сигнатуры для классификации с обучением. Особое внимание обратите на параметр maximum number of sub-signtures in any class на вкладке Optional, изменение которого существенно влияет на результат. Первоначальную классификацию можно провести с настройками по умолчанию, задав только опции на вкладке Required, а именно выбор слоя с обучающей выборкой, группы и подгруппы наборов классифицируемых данных и имени выходного файла. После получения первичного результата можно будет начать эксперименты с другими настройками.
i.gensigset trainingmap=classes group=atcorr_group subgroup=atcorr_subgroup signaturefile=smap_sign
Теперь запускаем модуль i.smap (Imagery - Classify image - Sequential maximum a posteriori classification (SMAP)). Задаем знакомые параметры на вкладке General. На вкладке Optional можно установить флаг "Use maximum likelihood estimation". Возможно, для ваших данных это будет решающим моментом (в любом случае поэкспериментируйте с флагом и без него).
i.smap group=atcorr_group subgroup=atcorr_subgroup signaturefile=smap_sign output=smap_classify
Как только мы получаем удовлетворительный результат (путем всевозможных комбинаций доступных параметров), переходим к постклассификации. В процессе, возможно, вы поймете, что следует модифицировать обучающие выборки.
Постклассификация
Постклассификация — это во многом субъективный, но важный этап. Мы воспользуемся методом "скользящего окна" для того, чтобы избавиться от всевозможных "шумов" (отдельных пикселей) и генерализовать полученный результат. Запускаем модуль r.neighbors (Raster - Neighborhood analysis - Moving window). Заполняем обязательные поля на вкладке Required (растр для фильтрации, имя для результирующего растра) и переходим к выбору фильтра (вкладка Neighborhood). В нашем случае, когда каждый пиксель имеет всего 5-7 (может быть меньше или больше, в зависимости от количества классов) возможных значений, определенных однозначно, и не может быть никаких промежуточных вариантов, больше всего подходит метод mode, который присваивает пикселю значение, которое чаще всего встречается среди его соседей. Результат сильно зависит от размера фильтра (кстати, значение может быть только нечетным, иначе будет ошибка). Для примера мы остановились на размере 7х7 (Neighborhood size = 7).
r.neighbors input=smap_classify output=map_1 method=mode size=7
В результате получаем заметно отфильтрованное изображение:
Также, если вы хотите сравнить результаты нескольких методов, то в Name of output можно через запятую перечислить несколько имен (по количеству используемых методов), а во вкладке Neighborhood установить несколько флагов (названия присваиваются соответственно по порядку методов).
Конечный результат можно, при необходимости, перевести в вектор уже знакомыми вам методами,. Подгрузив на карту векторный слой воды, мы получим практически готовую карту.
Карта
В сущности, задача решена, остается только оформить карту на выходе (если это требуется, конечно). GRASS не является идеальной средой для оформления макетов карт, поэтому мы использовали для этих целей открытую ГИС QGIS. В общем-то, можно использовать любой удобный для вас инструмент, вплоть до обыкновенного векторного редактора. Подразумевается, что мы произвели предварительное преобразование результатов постклассификации в векторный формат с помощью уже использовавшегося модуля r.to.vect. Далее из GRASS можно экспортировать (с помощью модуля v.out.ogr (File - Export Vector map - Common Export formats) этот векторный слой конечного результата постклассификации. При экспорте нежелательно устанавливать флаг напротив Also export features without category… (также экспортировать объекты без категории) на вкладке Selection, только если вы знаете, что полигоны без категорий вам тоже необходимы). Если вы экспортируете результаты в шейп-файл (ESRI SHP), то также не забудьте выгрузить и файл .prj, установив соответствующий флаг во вкладке Creation. Дополнительно мы экспортировали в шейп-файл слой воды, который создавался на этапе маскирования.
v.out.ogr -e input=cuba_map output=cuba format=ESRI_Shapefile
Экспортированные данные оформляются в удобной для вас среде. В нашем случае в QGIS была получена карта такого вида:
В случае использования для оформления именно QGIS можно пропустить шаг с экспортом данных модулем v.out.ogr и воспользоваться штатным модулем поддержки GRASS (доступен в модулях QGIS по умолчанию). Он позволяет открыть набор геоданных и добавлять на карту векторные и растровые слои.
