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

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


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

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

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

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

  • Работоспособность самой актуальной сборки GRASS GIS 7 на день публикации статьи является адекватной только в операционных системах семейства Linux и в Mac OS. Поведение пакета в Windows непредсказуемо и нестабильно, большая часть из описанных в статье операций в ней не работает (например, генерируются пустые файлы или программа завершает работу с ошибкой). При работе над статьей использовались Linux Debian (в VirtualBox) и нативный Mac OS Lion, GRASS GIS 7.0svn (r64224).
  • Все (или почти все) действия в статье выполняются через графический интерфейс GRASS, несмотря на то, что часто проще и быстрее использовать командную строку. Это сделано для увеличения наглядности и уменьшения порога вхождения. GRASS GIS всегда показывает текстовую команду, соответствующую выполняемым в интерфейсе действиям, поэтому, при желании, вы можете ознакомиться с этим аспектом работы подробнее.
  • Предполагается, что вы загрузили и распаковали набор данных Landsat 8 уровня обработки "Level 1 Product" (Например с помощью EarthExplorer, краткие инструкции можно получить здесь или здесь) с одиннадцатью одноканальными изображениями и файлом метаданных.


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

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


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


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


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


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


Выбор проекции


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


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


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


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

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


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


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


Окно импорта растровых данных


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


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


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


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


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


Меню определения региона

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

Первым важным шагом при обработке данных 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 для исходных и выходных снимков. Для исходных это строка, которая должна полностью совпадать с именами загруженных вами снимков, только без цифры канала, которую модуль будет подставлять автоматически. В рассматриваемом примере: LC81740292013131LGN01_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"


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


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


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


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

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

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

Pan Sharpening

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

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

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

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

Карта