Тайлы из Landsat 8

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

Упрощённая инструкция по получению и обработке снимков Landsat 8 и подготовке из них тайлов

Главное свойство данных ДЗЗ Landsat 8 — свободная лицензия, которая позволяет обрисовывать снимки для, в частности, OpenStreetMap. Эта статья не вдаётся в теорию, но содержит шаги для поиска и обработки снимков Landsat 8 и их превращения в тайлы для публикации или обрисовки в JOSM.

Получение снимка

Раньше нужно было шоркаться по EarthExplorer, вводить параметры поиска, что-то выбирать. Спасибо DevelopmentSeed, много лет работающей с Landsat, теперь выбрать нужный снимок проще. Заходите, приближайте свою область, кликайте в точку, фильтруйте и скачивайте:

http://libra.developmentseed.org

Напоминаем, что каждый снимок весит 500-800 мегабайт в архиве (2 гигабайта вне его), и для обработки потребуются ещё гигабайта четыре. Под линуксом файл распаковывается такой командой:

tar -xjvf LC*.tar.bz

В архиве вы получите 11 каналов и слой QA, на котором отмечены точки без информации о земной поверхности (например, облака). Мы будем работать с каналами 2, 3, 4 (синий, зелёный, красный) и 8 (ч/б с 15-метровым разрешением).

Паншарпенинг

У нас есть 3 канала с разрешением 30 м/пк и один — 15 м/пк. Паншарпенинг — это объединение их всех в один цветной снимок максимального разрешения. Раньше его делали длинным цепочками команд GDAL или ImageMagick, но нам повезло. Умные люди из GINA (университет Аляски) написали скрипт, и теперь достаточно одной команды. Набор аляскинских скриптов есть в большинстве дистрибутивов, поищите пакет dans-gdal-scripts. Команда такая:

gdal_landsat_pansharp -pan B8.TIF -ndv 0 -o result.tif \
-rgb B4.TIF -rgb B3.TIF -rgb B2.TIF -lum B3.TIF 0.5 -lum B4.TIF 0.5

Здесь опущены длинные идентификаторы снимков, вида «LC81850172014281LGN00». Для удобства пригодится такой скрипт:

BASE=LC81850172014281LGN00
gdal_landsat_pansharp -pan ${BASE}_B8.TIF -ndv 0 -o result.tif \
 -rgb ${BASE}_B4.TIF -rgb ${BASE}_B3.TIF -rgb ${BASE}_B2.TIF \
 -lum ${BASE}_B3.TIF 0.5 -lum ${BASE}_B4.TIF 0.5

Что тут происходит? Первые три параметра передают панхроматический канал (восьмой), значение «нет данных» и имя выходного файла. Следующая тройка с -ref — каналы для красного, зелёного и синего цветов. В Landsat 8 им соответствуют 4, 3 и 2, но можно выбрать другие, из невидимых областей спектра.

Наконец, параметры -lum корректируют наложение цвета на панхроматический канал. Описание алгоритма слишком сложное, мы взяли коэффициенты из этих функций Дерека Уоткинса.

Коррекция цвета

Автоматически

Белую дымку удалить непросто. В GDAL нет команды для обработки цвета, подобной вкладке в QGIS, но тот же аляскинском наборе содержит скрипт, выдающий приемлемый результат:

gdal_contrast_stretch -ndv 0 -percentile-range 0.05 0.95 result.tif result2.tif

Доли отбрасываемых цветов можно покрутить на ±0.03: чем ближе к 0-100, тем темнее получается картинка, но тем больше деталей сохраняется.

В графическом редакторе

Раз tif — обычная картинка с доп. полями, можно открыть его в ImageJ или ином редакторе и подкрутить контраст вручную, акцентируя нужные полосы. Чтобы сохранить геопривязку, понадобятся программы из пакета libgeotiff (иногда geotiff-bin). Экспорт геопривязки в файл:

listgeo result.tif > result.geo

Восстановление геопривязки для отредактированного файла (сохранять нужно, понятно, тоже в формате TIFF):

geotifcp -g result.geo edited.tif edited_geo.tif

Изображение

Если снимок нужен для иллюстрирования, а не для обрисовки, можно просто взять файл tif: этот формат читается всеми. Как вариант, с помощью ImageMagick можно из него сделать более удобный jpeg:

convert result.tif -quality 85 result.jpg

Тайлы

С помощью GDAL

В комплекте GDAL есть программа для создания тайлов из геопривязанного растра:

gdalwarp -t_srs EPSG:3857 result.tif result_merc.tif
gdal2tiles.py -r lanczos -z 6-13 -a 0 result_merc.tif tiles

Она работает очень медленно, сохраняет тайлы только в PNG и выдаёт странные ошибки во время работы. Альтернатива — QTiles.

В QGIS

Собранный растр можно обработать в QGIS, а не запускать дополнительные скрипты. Если не правили контраст скриптом, в свойствах растрового слоя, во вкладке «Стиль» нужно выбрать «Среднее ± отклонение × 2», нажать кнопку «Загрузить» — и «OK».

Отсюда же легко сделать тайлы с помощью модуля QTiles. Если его не видно в списке модулей, включите в настройках «также показывать экспериментальные модули». Если нет вкладки настроек — установите пакет qgis-python или подобный.

Тайлы рекомендуем делать в формате jpeg, до 13 масштаба. Демонстрационную html-страницу придётся поправить, заменив расширение тайлов на jpg. В JOSM тайлы подключаются так же, как и раньше, только поправьте расширение.

Подключение в JOSM

Чтобы подключить слой в JOSM, скопируйте строку для него из демонстрационной страницы, либо напишите сами:

tms[13]:file:///home/user/путь/к/тайлам/{zoom}/{x}/{y}.png

Для Windows путь записывается чуть по-другому:

tms[13]:file:///C:/путь/к/тайлам/{zoom}/{x}/{y}.png

Не забывайте проверять привязку по трекам.

См. также

• Сравнение каналов Landsat 7 и 8 Кристофом Хорманном
• Эксперименты с каналами
• Команды для Landsat 7
• Орторектифицирован ли снимок? (обычно да)
• Подготовка снимков в ImageMagick
• Паншарпенинг данных Landsat 8 средствами свободного ПО