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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
Нет описания правки
Строка 183: Строка 183:
Получаемые съемочной системой данные о подстилающей поверхности и объектах на ней, очевидно, являются искаженными из-за влияния множества факторов, среди которых основной - атмосфера, сложная и разнообразная. Ясно, что для корректной обработки необходимо максимально подробно учесть эти искажения, восстановив "истинное лицо" объектов на земной поверхности, именно для этого применяется так называемая '''Атмосферная коррекция'''. Чтобы учесть атмосферу, её нужно как-то смоделировать, или, выражаясь корректнее, смоделировать поведение отраженной солнечной радиации на пути от объекта до съемочной аппаратуры. В GRASS (модуль i.atcorr) реализована довольно сложная и подробная модель 6S (Second Simulation of a Satellite Signal in the Solar Spectrum). [http://6s.ltdri.org/pages/manual.html Документация модели]. Ей мы и попробуем воспользоваться.
Получаемые съемочной системой данные о подстилающей поверхности и объектах на ней, очевидно, являются искаженными из-за влияния множества факторов, среди которых основной - атмосфера, сложная и разнообразная. Ясно, что для корректной обработки необходимо максимально подробно учесть эти искажения, восстановив "истинное лицо" объектов на земной поверхности, именно для этого применяется так называемая '''Атмосферная коррекция'''. Чтобы учесть атмосферу, её нужно как-то смоделировать, или, выражаясь корректнее, смоделировать поведение отраженной солнечной радиации на пути от объекта до съемочной аппаратуры. В GRASS (модуль i.atcorr) реализована довольно сложная и подробная модель 6S (Second Simulation of a Satellite Signal in the Solar Spectrum). [http://6s.ltdri.org/pages/manual.html Документация модели]. Ей мы и попробуем воспользоваться.


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


Хотелось бы отметить, что проведение процедуры атмосферной коррекции не является гарантом получения достоверных данных об объектах. Откровенно говоря, никто не может обещать, что результат не окажется ещё более искаженным, чем исходный. Но всё же, при аккуратной и внимательной настройке модуля атмосферной коррекции, ''обычно'' результат оказывается более приемлемым.
Хотелось бы отметить, что проведение процедуры атмосферной коррекции не является гарантом получения достоверных данных об объектах. Откровенно говоря, никто не может обещать, что результат не окажется ещё более искаженным, чем исходный. Но всё же, при аккуратной и внимательной настройке модуля атмосферной коррекции, ''обычно'' результат оказывается более приемлемым.
В соответствии с [http://grass.osgeo.org/grass64/manuals/i.atcorr.html документацией] к модулю i.atcorr для начала необходимо создать конфигурационный файл для метода 6S. Для улучшения результатов предварительно добудем данные по рельефу и атмосферным условиям на исследуемую территорию (модель 6S во многом опирается на эти данные, хотя они и являются опциональными).
Рельеф можно брать [http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp отсюда] (SRTM). Импортируем его как обычный растровый набор данных.
Состояние атмосферы [http://aeronet.gsfc.nasa.gov/cgi-bin/type_piece_of_map_opera_v2_new здесь] и [http://gdata1.sci.gsfc.nasa.gov/daac-bin/G3/gui.cgi?instance_id=aerosol_daily здесь]. В частности, нам отсюда нужен параметр концентрации аэрозолей для диапазона 550 нм (только число, наборы данных отсюда не нужны).
В качестве одного из параметров необходимо будет указать среднюю высоту (помимо самой ЦМР). Рассчитаем её с помощью модуля r.univar. Запускаем его через вкладку Command Console:
[[Файл:Grass7_landsat8_processing_22.jpg|center|thumb|500px|Запуск модуля r.univar]]
В появившемся окне выбираем слой с ЦМР
[[Файл:Grass7_landsat8_processing_23.jpg|center|thumb|500px|Модуль r.univar. Выбор растрового слоя с ЦМР]]
Исполняем модуль и отслеживаем результат работы во вкладке command output. Запоминаем (записываем) значение средней высоты (mean).
[[Файл:Grass7_landsat8_processing_24.jpg|center|thumb|500px|Результат работы модуля r.univar]]
Также нам понадобятся координаты центра снимка. Если набор исходных данных был загружен с сайта EarthExplorer, то при выборе снимка указывались географические координаты центра сцены. Можете быстро найти свою сцену по имени.
[[Файл:Grass7_landsat8_processing_25.jpg|center|thumb|500px|EarthExplorer. Координаты центра сцены]]
Создаем файл с параметрами для 6S (например atcorr_pars.txt). Это самый обыкновенный текстовый документ, где одна строка - один параметр (в жестком порядке). Разные каналы корректируются отдельно (то есть процедуру придется повторить несколько раз). Рассмотрим файл конфигурации для нашей Кубинской сцены, синий канал (B2).
<pre>
18 - читаем следующую строку как положение спутника Landsat8
05 11 15.48 -75.722 20.187 - месяц| день| время GMT| долгота | широта
2 - атмосферная модель
1 - аэрозольная модель
0 - не указываем видимость
0.014 - параметр концентрации аэрозолей для диапазона 550 нм
-0.244 - средняя высота над уровнем моря в км(с обратным знаком)
-1000 - означает, что сенсор находится на борту спутника
116 - индекс соответствующего канала landsat 8 (здесь - синий)
</pre>
** обратите внимание, что индекс панхроматического канала 119 между red и nir
Первая строка это "Geometrical conditions". В сущности здесь - разные сенсоры. Landsat 7 (ETM+) имеет код 8, а используемый нами Landsat 8 (OLI) - код 18. В зависимости от указанного здесь кода формат второй строки будет отличаться (см. таблицу в документации).
Во второй строке (Для Landsat 8) указывается время (GMT) и координаты центра сцены в десятичных долях с соответствующим знаком для широты и долготы. Время извлекается из файла метаданных (строка SCENE_CENTER_TIME):
<pre>
SCENE_CENTER_TIME = 15:28:49.4315057Z
</pre>
Координаты - указанным выше способом.
Далее следуют важные параметры, требующие от вас некоторой смекалки. Это атмосферная модель, выбираемая из семи стандартных вариантов:
<pre>
0 Поглощение отсутствует
1 Тропики
2 Средние широты (лето)
3 Средние широты (зима)
4 Субарктические широты (лето)
5 Субарктические широты (зима)
6 Стандарт US 62
</pre>
либо можно задать собственную модель на основе, к примеру, измерений радиозондом.
Следующая строка - аэрозольная модель. Её также можно выбрать из семи стандартных вариантов:
<pre>
0 Нет аэрозолей
1 Континентальная модель
2 Приморская модель
3 Модель для урбанизированных территорий
4 Модель для пустыни (Shettle)
5 Горение биомассы (biomass burning)
6 Стратосферная модель
</pre>
Следующий параметр
Подробно обо всех этих параметрах можно прочесть в документации i.atcorr.





Версия от 00:01, 2 марта 2015

Эта страница является черновиком статьи.


Использование открытой геоинформационной системы 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"


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


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


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


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

Получаемые съемочной системой данные о подстилающей поверхности и объектах на ней, очевидно, являются искаженными из-за влияния множества факторов, среди которых основной - атмосфера, сложная и разнообразная. Ясно, что для корректной обработки необходимо максимально подробно учесть эти искажения, восстановив "истинное лицо" объектов на земной поверхности, именно для этого применяется так называемая Атмосферная коррекция. Чтобы учесть атмосферу, её нужно как-то смоделировать, или, выражаясь корректнее, смоделировать поведение отраженной солнечной радиации на пути от объекта до съемочной аппаратуры. В 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.014					- параметр концентрации аэрозолей для диапазона 550 нм
-0.244					- средняя высота над уровнем моря в км(с обратным знаком)
-1000					- означает, что сенсор находится на борту спутника
116					- индекс соответствующего канала landsat 8 (здесь - синий)
    • обратите внимание, что индекс панхроматического канала 119 между red и nir


Первая строка это "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	Стратосферная модель	 

Следующий параметр


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


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

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

Pan Sharpening

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

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

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

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

Карта