Создаем бесшовную мозаику из космических снимков в ArcGIS Desktop

В ходе решения практической задачи по созданию единого покрытия космоснимками Landsat-8 2017 года съёмки на территорию Ямало-Ненецкого АО (ЯНАО) был получен собственный опыт решения такой задачи с использованием настольного ПО ArcGIS Desktop, которым спешим поделиться со всеми!

Начнём с постановки задачи. Создать подложку (базовую карту) территории ЯНАО по снимкам Landsat-8 за 2017 год с выполнением цветового уравнивания и построением линий сшивок для минимизаций различий между отдельными снимками, базовая карта должна работать как веб-сервис ArcGIS Server. Требуется использовать максимальное разрешение – 15 метров и минимизировать облачность на снимках.

Конечным результатом должен стать кэш из снимков в масштабах 1:18 489 298 — 1:72 224 (тайловая схема Google Maps|ArcGIS Online|Bing Maps)

Ожидаемый результат после подготовки мозаики по снимкам Landsat-8

Подобная задача комплексная и начинается она с подбора необходимых исходных данных.

У Landsat постоянные орбиты, т.е. можно выбрать оптимальное покрытие заранее по номерам орбит, что и было сделано, но…облачность – потребовалось набирать перекрытие с других орбит за разные даты съемки. Для самостоятельной оценки можно скачать шейп-файл с покрытием снимками с разных маршрутов съёмки Landsat-8 на весь земной шар: https://landsat.usgs.gov/pathrow-shapefiles  (орбита Descending only, WRS 2).

Теоретическое минимальное число снимков для полного покрытия региона и реальный результат после выбора снимков с минимальной облачностью за 2017 год

Поиск и скачивание снимков. Для поиска и скачивания был выбран сайт EarthExplorer https://earthexplorer.usgs.gov/

Причина выбора ресурса EarthExplorer:

  • возможность пакетной загрузки отобранных снимков с поддержкой опции докачки
  • удобная настройка критериев поиска
  • в файлах со снимками доступны метаданные по облачности снимков, что потребуется в дальнейшем при сортировке в наборе данных мозаики. Метаданные автоматически добавляются средствами ArcGIS в атрибутивную таблицу к растрам
  • метаданные позволят загружать в мозаику сразу пан-шарпенинг (pan-sharpening) изображения без их предварительного ручного создания

В итоге: объём загруженных снимков в архиве ~ 60 Гб, после распаковки и построения пирамид ~ 160 Гб. Загружайте снимки сразу в рабочую папку, чтобы не гонять большой массив данных по сети. Для пакетной загрузки потребуется установить программу https://earthexplorer.usgs.gov/bulk

Несколько слов о поиске снимков. Достаточно творческий процесс по выбору безоблачных данных. Поиск не по району интереса, а по отобранным Path и Row на первом шаге. Поиск только за 2017 год и только за три месяца (установлено опытным путем, что для данного региона эти месяцы оптимальны): июнь, июль, август. Для непокрытых снимками областей при первой выборке или областей с высоким процентом облачности далее итеративно добирались снимки для создания избыточного покрытия. Снимки добавлялись в общий заказ и потом скачивались пакетно.

Интерфейс поиска веб-сайта EarthExplorer

Переходим в ArcGIS Desktop – строим пирамиды и статистику. Используем инструмент «Build Pyramids And Statistics» («Построить пирамидные слои и статистику»): инструмент позволяет указать на вход корневую папку для построения пирамид и статистики для всех растров в подпапках. Работа инструмента может длиться до пары часов.

Создаем БГД и набор данных мозаики, загружаем описания. Для набора данных мозаики (mosaic dataset) выбираем проекцию Web Mercator. Загружаем описание снимков в мозаику – добавляем файлы _MTL.txt из каждой папки со снимком.

Примечание: после добавления описаний снимков в мозаику в таблице атрибутов к слою Footprints будут добавлены все необходимые метаданные о снимке: дата съёмки (Acquisition Date), облачность кадра (Cloud Cover), номер снимка (Name), номер маршрутка и кадра (WRS_PATH, WRS_ROW)

Подготовка набора данных мозаики

Перестраиваем границы снимков (Footprints). Границы снимков надо обновить после первоначального добавления исходных данных.

Удаляем лишние растровые продукты из мозаики. В мозаику автоматически добавляется Multispectral и Pansharpened растровые продукты. Для окончательной мозаики будут использоваться только Pansharpened продукт (цветные снимки с разрешением 15 метров в случае Landsat-8), поэтому требуется удалить Multispectral вариант из мозаики.

Отображение мозаики на всех масштабах. Для отображения растров на мелких масштабах открываем таблицу атрибутов для Footprints и указываем для значения MaxPS значение 10000 или любое число больше.

Максимальное число растров в мозаике. По умолчанию будет отображаться только 20 снимков из набора данных мозаики. В ArcCatalog, в свойствах набора данных мозаики указываем максимальное число отображаемых растров (по умолчанию стоит 20). Это число должно быть больше числа снимков, добавленных в мозаику.

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

Вырезание по границам региона – используем растровые функции. В свойствах мозаики через ArcCatalog переходим на вкладку функции (Functions) и добавляем функцию Clip (вырезание).

Построение линий сшивок. Важно указать сортировку по атрибуту Облачность (CloudCover), чтобы безоблачные снимки лежали выше облачных в областях их перекрытия. Также задаем значения смешивания цветов (Blend) на границе линии сшивки.

Примечание: В проекте использовались линии сшивки по методу COPY_FOOTPRINT

Пример линии сшивки по методу RADIOMETRY:

Если вести речь идет о построении линии сшивки по методу COPY_FOOTPRINT, т.е. когда берется рамка кадра и используется в качестве линии сшивки, то подход следующий и обязательно включающий ручную правку геометрии.

Модернизируем линию сшивки. После построения линии сшивки по методу COPY_FOOTPRINT потребуется немного уменьшить рамки снимков внутрь, чтобы по краям не возникало черных полос. Строим отрицательный буфер (500 – 700 метров) для слоя Seamline в наборе данных мозаики. Обновляем линии сшивки путем связывания старой линии сшивки и новой по полю RasterID.

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

Примечание: линия сшивки редактируется как обычный векторный слой. Для перестроения мозаики потребуется сначала сохранить изменения в редактировании.

Настраиваем свойства набора данных мозаики в ArcCatalog. 

  • Default Resampling Method – Cubic
  • Allowed Mosaic Methods – только Seamline (стыковка снимков по линии сшивки)
  • Default Sorting Order – Ascending
  • Default Mosaic Operator – Mean
  • Blend Width Unit – Pixels
  • Blend Width – 30 или другое значение

Цветовое уравнивание. Используем метод COLOR_GRID, который хорошо работает по большим областям.

Перед запуском цветового уравнивания снова запускаем построение пирамидных слоев и статистики с помощью инструмента «Build Pyramids And Statistics» («Построить пирамидные слои и статистику»), но на вход указываем уже набор данных мозаики. Для ускорения процесса можно пропускать больше 1 строки и столбца при построении статистики (1 — значение по умолчанию в инструменте, можно указать 100 или другое число).

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

Построение обзорных изображений. Строим обзорные изображения со стандартными настройками для ускоренного отображения мозаики на обзорных масштабах.

И наступает самая приятная часть  — результат проделанной работы.

Масштаб 1:18 489 298:

Бесшовная мозаика ЯНАО, построенная по снимкам Landsat-8 (15 метров), полученным в летний период 2017 года.

Масштаб 1:72 224:

Максимальное разрешение для построения кэша на основе снимков Landsat-8 (15 метров), полученных в летний период 2017 года.

Кэширование. Для быстрого отображения данных на ArcGIS Server запускаем кэширование мозаики средствами ArcGIS Server. Выставляем обозначенные масштабы 1:18 489 298 — 1:72 224 (тайловая схема Google Maps|ArcGIS Online|Bing Maps). Перед началом кэширования в таблице содержания выбираем изображение набора данных мозаики и в его свойствах меняем на вкладке Mosaic оператор на Blend.

Примечание: убедитесь, что исходные снимки и набор данных мозаики находятся в зарегистрированном для ArcGIS Server местоположении, иначе перед началом кэширования ваши 160Гб+ данных начнут автоматически копироваться на сервер!

Удачи в построении собственной мозаики! Надеемся, что наш опыт будет полезен.