В ходе решения практической задачи по созданию единого покрытия космоснимками 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 постоянные орбиты, т.е. можно выбрать оптимальное покрытие заранее по номерам орбит, что и было сделано, но…облачность – потребовалось набирать перекрытие с других орбит за разные даты съемки. Для самостоятельной оценки можно скачать шейп-файл с покрытием снимками с разных маршрутов съёмки Landsat-8 на весь земной шар: https://landsat.usgs.gov/pathrow-shapefiles (орбита Descending only, WRS 2).
Поиск и скачивание снимков. Для поиска и скачивания был выбран сайт EarthExplorer https://earthexplorer.usgs.gov/
Причина выбора ресурса EarthExplorer:
- возможность пакетной загрузки отобранных снимков с поддержкой опции докачки
- удобная настройка критериев поиска
- в файлах со снимками доступны метаданные по облачности снимков, что потребуется в дальнейшем при сортировке в наборе данных мозаики. Метаданные автоматически добавляются средствами ArcGIS в атрибутивную таблицу к растрам
- метаданные позволят загружать в мозаику сразу пан-шарпенинг (pan-sharpening) изображения без их предварительного ручного создания
В итоге: объём загруженных снимков в архиве ~ 60 Гб, после распаковки и построения пирамид ~ 160 Гб. Загружайте снимки сразу в рабочую папку, чтобы не гонять большой массив данных по сети. Для пакетной загрузки потребуется установить программу https://earthexplorer.usgs.gov/bulk
Несколько слов о поиске снимков. Достаточно творческий процесс по выбору безоблачных данных. Поиск не по району интереса, а по отобранным Path и Row на первом шаге. Поиск только за 2017 год и только за три месяца (установлено опытным путем, что для данного региона эти месяцы оптимальны): июнь, июль, август. Для непокрытых снимками областей при первой выборке или областей с высоким процентом облачности далее итеративно добирались снимки для создания избыточного покрытия. Снимки добавлялись в общий заказ и потом скачивались пакетно.
Переходим в 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:
Масштаб 1:72 224:
Кэширование. Для быстрого отображения данных на ArcGIS Server запускаем кэширование мозаики средствами ArcGIS Server. Выставляем обозначенные масштабы 1:18 489 298 — 1:72 224 (тайловая схема Google Maps|ArcGIS Online|Bing Maps). Перед началом кэширования в таблице содержания выбираем изображение набора данных мозаики и в его свойствах меняем на вкладке Mosaic оператор на Blend.
Примечание: убедитесь, что исходные снимки и набор данных мозаики находятся в зарегистрированном для ArcGIS Server местоположении, иначе перед началом кэширования ваши 160Гб+ данных начнут автоматически копироваться на сервер!
Удачи в построении собственной мозаики! Надеемся, что наш опыт будет полезен.
