как часто можно получать свежие снимки sentinel
Как обрабатывать спутниковые снимки с помощью Sen2Cor
Sen2Cor — программа для обработки снимков, сделанных со спутника Sentinel-2. В статье рассказывается, как установить, запустить и настроить её.
В магистратуре мне понадобилось сделать атмосферную коррекцию снимков со спутника Sentinel-2. Магистратура по экологии, поэтому можно было не углубляться в физику и математику, а использовать готовый инструмент. Такой инструмент нашёлся — Sen2Cor.
Проблема в том, что понятной инструкции для Sen2Cor нет. Авторы учебников по геоинформатике не лезут в такие дебри, как работа с конкретными программами, а авторы научных статей, напротив, полагают, что коррекция снимков — слишком простая часть методики, поэтому не описывают её подробно. На сайтах «для чайников» (и не совсем для чайников) о Sen2Cor тоже не пишут. Официальная документация запутанная: достаточно сказать, что раздел «Установка» на 26-й странице «Руководства пользователя» сразу отсылает к 60-й странице Release Notes. Чтобы разобраться, мне потребовалось время.
Я решил написать понятное руководство по установке, запуску и настройке Sen2Cor. Руководство ориентировано на биологов, экологов, специалистов по сельскому хозяйству и всех тех, кому может понадобиться обработка спутниковых снимков, но для кого она не является основной задачей. По сути, я делюсь собственным опытом. Абсолютную правильность и полноту не гарантирую. Специалисты по геоинформатике, думаю, сами мне что-то посоветуют.
Строго говоря, Sen2Cor делает не только атмосферную коррекцию, но и поправку на угол восхождения Солнца и рельеф местности, а также тематическую классификацию (scene classification), но для меня важнее всего была атмосферная коррекция. Поэтому для простоты я иногда вместо «обработка» пишу «атмосферная коррекция», хотя это неточно.
Sen2Cor работает в трёх режимах: как отдельное консольное приложение, как плагин в составе Sentinel-2 Toolbox и как консольное приложение, используемое в PDGS (Payload Data Ground Segment). PDGS — это наземный центр обработки данных, получаемых со спутников. Находится он в Италии в исследовательском институте. Теоретически можно запустить Sen2Cor в том же режиме, в котором его запускают в центре обработки данных, но для запуска понадобятся данные, которые просто так в интернете не скачать, поэтому мы данный вариант рассматривать не будем. Sentinel-2 Tolobox тоже обойдём стороной, потому что с этой программой я не работал. Она с графическим интерфейсом, поэтому полагаю, что разобраться в её работе несложно. Итак, в данной статье описано, как использовать Sen2Cor в режиме консольной программы.
Sen2Cor работает в 64-битных Windows, MacOS и Linux. У меня Linux, поэтому команды буду показывать на примере Linux. Впрочем, разница только в установщиках и в путях к файлам.
Sen2Cor требует 4 Гб оперативной памяти. Если меньше, то он запустится, но может аварийно завершиться на середине обработки. Расход оперативной памяти зависит от разрешения производимых снимков, поэтому если памяти мало и программа из-за этого выдаёт ошибку, то можно запустить обработку с меньшим разрешением. О том, как это сделать, написано далее. Аналогичных жёстких требований к процессору нет, но чем слабее процессор, тем больше времени займёт обработка.
Установка
Нужно скачать установщик и запустить его. Sen2Cor написан на Python, однако, интерпретатор и все необходимые пакеты находятся в установщике, поэтому самому устанавливать Python и настраивать окружение не нужно.
Далее открываем консоль, перемещаемся в каталог spaceshots и запускаем установщик.
Если ошибок нет, то в консоли должна появиться надпись «Congratulations, Installation successful. » и предложение протестировать установку. Сделаем тестовый запуск:
В MacOS программа устанавливается так же, как в Linux.
Обработка снимка
После нажатия на Enter можно отдохнуть, потому что обработка занимает довольно много времени (на Яндекс.Облаке с Intel Cascade Lake и гарантированной долей 5% CPU — около часа). Sen2Cor выводит в консоль информацию о ходе работы, так что вы будете примерно понимать, что в данный момент делает программа.
Когда программа завершит работу, в /home/user/spaceshots рядом с S2A_MSIL1C_20180727T071621_N0206_R006_T40VFJ_20180727T092607.SAFE/ должна появиться ещё одна папка с очень похожим названием — S2A_MSIL2A_20180727T071621_N9999_R006_T40VFJ_20200507T083341.SAFE
Теперь разберёмся, что делает программа и что за папку мы получили. Для этого нужно рассмотреть структуру снимков Sentinel-2.
Sen2Cor получает в качестве исходных данных папку *.SAFE — именно она указывается в качестве единственного обязательного параметра при запуске программы. Если вместо папки *.SAFE указать что-то другое, например, исходный архив со снимком или конкретный тайл, то программа не запустится.
Дело в том, что для атмосферной коррекции используются не только сами тайлы, но и метаданные, а они находятся в файлах внутри папки *.SAFE, а не в тайлах.
Результатом работы программы также является папка *.SAFE с похожей структурой: внутри есть папка GRANULE, внутри неё — ещё одна папка, в которой находится папка IMG_DATA, а в ней скорректированные тайлы. Тайлы сгруппированы по пространственному разрешению: в папке R10m находятся те, у которых пространственное разрешение составляет 10 метров, в R20m — те, где разрешение 20 метров, в R60m — 60 метров. При определённых настройках какие-то разрешения не обрабатываются и соответствующие им папки не создаются — об этом чуть ниже.
Название папки *.SAFE, получающейся в результате, отличается от названия исходной только одним: в начале вместо _MSIL1C_ указано _MSIL2A_. L1C и L2A — это коды уровня обработки. Уровень L1C означает, что снимок обработан до уровня Top-of-Atmosphere reflectance, то есть значения пикселей на снимке соответствуют отражению, зафиксированному у верхней границы атмосферы — там, где летает спутник. Уровень L2A — это обработка до уровня Bottom-of-Atmoshpere reflectance, то есть до значений отражения у поверхности Земли — так, как будто атмосфера абсолютно прозрачна. В реальности это, конечно, не так, и атмосферная коррекция позволяет убрать только часть искажений, возникающих при прохождении света через атмосферу, но в целом значения Bottom-of-Atmosphere reflectance больше соответствуют реальному отражению от земной поверхности, чем значения Top-of-Atmosphere reflectance.
Помимо скорректированных тайлов, Sen2Cor производит ещё несколько изображений: полноцветное (TCI), карты аэрозольной оптической плотности атмосферы (AOT — Aerosol Optical Thickness) и содержания водяного пара (WV — Water Vapour), карту тематической классификации типов поверхности (SC — Scene Classification). Они также находятся в папке IMG_DATA в каталогах для соответствующего пространственного разрешения.
При запуске с такими параметрами Sen2Cor обработает снимок и создаст скорректированные тайлы с разрешением 60 метров. В отличие от стандартного запуска, обработка занимает намного меньше времени — чуть больше четырёх минут. Расход памяти тоже меньше, поэтому если запуск со стандартным разрешением закончился ошибкой из-за нехватки памяти, то можно указать разрешение 60 метров. В папке IMG_DATA будет находиться только папка R60m, а внутри неё — тайлы с разрешением 60 метров и некоторые другие изображения: полноцветное, тематическая карта и карта содержания водяного пара.
Конфигурационный файл
В конфигурационном файле много разных параметров. Рассмотрим некоторые из них.
Log_Level — параметр, который определяет, насколько подробным будет журнал работы программы. Лог отображается в консоли во время работы и дополнительно сохраняется в папке для логов (в нашем случае — /home/user/sen2cor/2.8/log ). Значение по умолчанию — INFO. Обычно нет необходимости менять эту настройку, если вы не столкнулись с ошибками в работе программы.
Учёт рельефа и угла восхождения Солнца
DEM_Reference — URL, с которого нужно скачать цифровую модель рельефа, если её нет в папке, указанной в предыдущем пункте. Разработчики Sen2Cor предлагают использовать — это ссылка на STRM DEM, свободно распространяемую цифровую карту рельефа с 90-метровым пространственным разрешением. Учитывая, что масштаб ненамного мельче, чем у самих спутниковых снимков, эта модель должна давать хороший результат.
Опция Generate_DEM_Output позволяет получить в результате работы программы отдельный тайл с цифровой картой рельефа. По умолчанию FALSE.
Опция DEM_Terrain_Correction частично отключает использование цифровой карты рельефа: рельеф по-прежнему будет учитываться при тематической классификации (SC) и построении карты AOT, но не при корректировке значений отражения от поверхности.
Если DEM не используется, то следует указать параметр Altitude — это средняя высота над уровнем моря в районе, запечатлённом на снимке. Высота указывается в километрах.
Красивая картинка
Параметр Generate_TCI_Output включает и выключает создание полноцветного изображения. По умолчанию TRUE, но если красивая картинка не нужна, то можно выбрать FALSE.
Учёт состояния атмосферы
Удаление облачности
Как увидеть на снимке лес? Наш опыт сегментации снимков Sentinel-2
Перед капитальной застройкой большой территории необходимо её детально исследовать. В зависимости от вида участка серьёзно варьируется стоимость строительства, предварительной обработки местности и многих других сопутствующих работ.
Чтобы серьезно минимизировать издержки на строительство объектов, некоторые компании используют спутниковые снимки с размеченной территорией. Такой метод исследования участка в разы дешевле проведения экспедиций, с помощью которых и в настоящее время часто решают эту задачу.
Но как получить полностью размеченные по необходимым критериям спутниковые снимки?
Мы в компании «Цифровое проектирование» решаем различные задачи сегментации. На каждой подробно остановимся в наших следующих материалах. В этом же тексте расскажем про наш опыт нахождения леса на спутниковых снимках, как пытаемся построить ML модель по разметке леса и про то, что уже удалось добиться на данном этапе, а к чему мы еще только стремимся.
Какая перед нами стоит задача?
Нам необходимо разработать алгоритм автоматической разметки леса довольно больших по площади областей для последующего анализа территории.
Помимо того, что нам необходимо разметить, где находится лес, нужно ещё классифицировать его на хвойный и лиственный. В дальнейшем планируется проводить оценку биомассы, для чего также будут нужны снимки с размеченным лесом. Важное условие на текущем этапе — распознавание леса должно проводиться по открытым данным, которые можно найти в интернете.
Мы использовали многоспектральные снимки спутников дистанционного зондирования Земли Sentinel-2, которые созданы Европейским космическим агентством. Спутники были запущены в 2015 году и продолжают функционировать до сих пор.
Почему именно Sentinel-2? В этом решении мы руководствовались несколькими, важными для нас, факторами.
Новые снимки Sentinel появляются достаточно часто. Для одной и той же территории снимок может обновиться уже через 5 дней. К тому же, они распространяются свободно, что для нас, как вы уже поняли, было критически важно. Также есть интерфейс, с помощью которого можно скачать нужный снимок. Правда, не очень удобный, но главное — он есть. Можно найти и API под Python и CLI, который позволяет легко автоматизировать запросы даже к архивным данным Sentinel.
Большое преимущество Sentinel-2 — они видят 12 спектральных каналов. Спутники видят не только каналы видимого спектра (синий, зеленый и красный), но также фиксируют невидимые глазу спектральные диапазоны (например, NIR и SWIR), что в контексте распознавания леса может быть довольно критично.
Ещё одним преимуществом снимков Sentinel-2 является их высокое разрешение в рамках основных RGB каналов и ближнего инфракрасного диапазона — 10 метров на пиксель. По остальным каналам, правда, разрешение хуже — 20 или 60 метров на пиксель, но для нас они были уже не столь важны. Для сравнения — у некоторых других агентств снимки по каналам RGB + ближний инфракрасный идут с разрешением 30 метров на пиксель, что серьезно усложняет задачу построения точных контуров леса. Коммерческие продукты, такие как WorldView-3, могут предложить существенно лучшее разрешение (
30см панхроматический канал, 1.24 м — RGB). Но, как уже упоминалось, такие продукты на текущем этапе проекта не используются.
Как решаются подобные задачи?
Классический подход заключается в расчёте множества различных индексов и сравнении их с некоторыми эмпирически подобранными пороговыми значениями.
В нашей работе мы использовали один из самых известных индексов для определения леса и прочей растительности — нормализованный разностный вегетационный индекс (NDVI — normalized difference vegetation index). Это числовой показатель качества и количества растительности на участке.
NIR — ближний инфракрасный, Red — красный канал
А почему вегетационный? Потому что значения индекса определяются тем, как растения отражают и поглощают световые волны разной длины. Если мы посмотрим на график отражения листа дерева в зависимости от длины волны, то в красной области спектра будет максимум поглощения солнечной энергии хлорофиллом, а в ближнем инфракрасном спектре, который уже не виден человеческим глазом, находится максимум отражения солнечного света клеточными структурами листа.
Поэтому, чем больше разница между NIR и Red, тем большее значение количества хлорофилла находится в соответствующем пикселе. Взаимосвязь этих каналов друг к другу позволяет четко отделять растительность от иных природных объектов.
Этот индекс настолько широко и так давно используется, что про него даже можно найти мемы.
Ищем данные
Для решения поставленной задачи мы использовали методы машинного обучения. Нам необходимо было создать две модели. Одна из них по снимку определяла бы, где находится лес, а другая сегментировала его на хвойный и лиственный.
Благодаря такому подходу нам не нужно самостоятельно подбирать сложные правила взаимодействия всех индексов — модель это сделает за нас, на основании данных. Данных, которые предстояло найти…
Мы занялись поисками территории, для которой уже известно — где лес, а где не лес. Первые данные, которые мы обнаружили, нас сильно воодушевили — это была полностью размеченная карта Словении.
Сразу же стал очевиден и первый камень преткновения. Словения не очень большая страна, особенно в сравнении с Россией или даже отдельными субъектами РФ. Однако, как оказалось позднее — это совсем не главная проблема.
По найденным данным была обучена модель, и мы получили очень высокие значения метрики intersection-over-union (метрика степени пересечения между двумя множествами). Грубо говоря — точность определения леса этой моделью на той части разметки, которая не участвовала в процессе обучения. И она составила 91%!
Не буду долго тянуть и скажу сразу — в дальнейшем значение этой метрики у нас только падало. Однако в процессе работы результаты становились осмысленней.
Важно отметить, что мы использовали несколько аугментаций (техника искусственного увеличения исходного набора данных) датасетов — горизонтальное / вертикальное отражение и случайный поворот на 90 градусов. С гаммой, яркостью и прочими классическими аугментациями на начальном этапе было принято решение не работать.
«Естественной» аугментацией также послужили снимки за разные даты для одной и той же территории, которые мы использовали в работе. Какие-то искажения это явно дало, но они максимально естественные (листва желтеет, появляется снег и т.п.).
В качестве модели мы взяли сверточную нейронную сеть, потому что она позволяет «из коробки» учитывать как взаимное пространственное расположение пикселей, так и расположение «в глубину» (все спектральные каналы). Основная архитектура модели — U-Net с кодировщиком ResNet-34.
Первые результаты давали большие надежды — в отдельных случаях только на основе разметки модель могла делать предсказания лучше, чем сама разметка. Хотя последняя проходила много стадий валидации, в том числе проводились и экспедиции.
По изображениям вы можете увидеть, что модель сразу же начала отбивать правильно те места, где леса нет.
Ground Truth маска
Предсказание модели
Но в применении к участку, которого, напомним, в тренировочной выборке не было, получился совершенно бессмысленный результат. Если ещё поиграться с чувствительностью модели либо нормализацией кадров, то появлялись некоторые пятна. В основном же модель показывала, что либо везде находится лес, либо леса нигде нет. Важно также подчеркнуть — дело не в несбалансированности целевых классов, мы проверили это, регулируя соответствующие веса целевой функции.
Основная проблема в том, что леса в Словении совсем не похожи на леса в России (а именно — Сибири). Там растут другие деревья, на других почвах и их характеристики в разных спектральных каналах не похожи на то, что можно наблюдать в Сибири.
Могли возникнуть предположения, что спутник снимал под другим углом и были искажения, из-за чего и появилась разница. Однако у Sentinel серьёзная процедура нормировки данных — и по уровню яркости, и по углам. Есть понятные умозрительные представления, почему критическим является именно состав леса и состав почвы. Достоверно известно, что спектральные свойства разных почв (разной влажности, разного состава) отличаются.
Никаких размеченных данных по Сибири нам найти не удалось, поэтому было решено добавить к имеющейся разметке Словении кусочки ручной разметки (т.е. выполненной своими силами) по регионам, которые похожи на Сибирь. К чему мы с новым вдохновением и приступили.
Размечаем вручную
Чтобы улучшить предсказания по Сибири и сохранить масштабируемость, для ручной разметки мы выбрали участки в Иркутской и Архангельской областях. Области были определены из соображений простоты ручной разметки и наибольшего сходства с интересующим нас регионом. Мы решили провести независимую разметку тремя людьми, что должно было помочь компенсировать систематические ошибки в работе каждого.
Разметка Иркутской области
Разметка Архангельской области
На этом этапе у нас возникло несколько дилемм, которые необходимо было решить.
Считать ли лесом отдельно стоящее дерево? Поскольку дальнейшим этапом работы модели будет оценка биомассы, функционально мы решили считать лесом более-менее плотно стоящие друг к другу деревья. Поэтому, при разметке отдельное дерево мы могли пропустить.
Как добиться того, чтобы все участники процесса трактовали снимок одинаково? Практически невозможно прийти к универсальному гайдлайну «что такое лес, а что такое не лес». Поэтому наша идея заключалась в диверсификации подходов. Попытки по-разному комбинировать наборы данных показали, что наши систематические ошибки могут хорошо друг друга компенсировать.
Чтобы исключить вариант обучения модели работать с конкретным кадром и конкретным регионом, мы не стали брать для ручной разметки именно тот регион, с которым в дальнейшем будем работать. В таком случае сохранялась возможность для последующей обобщаемости нашей модели.
Поскольку те небольшие участки, что мы разметили, нельзя сравнить по объему с имеющимся датасетом Словении, на этом этапе параллельно мы разрабатывали модели на основе Random Forest, которые можно использовать для классификации каждого отдельного пикселя без учета «оконных признаков». Даже при небольшой ручной разметке «набегает» приличное количество размеченных пикселей, которые складываются в солидный для обучения RF датасет.
И снова первые результаты были очень воодушевляющие! Модель, обученная только по Словении, давала бессмысленные результаты. Но модель, обученная по совокупности выборок, оказалась способна исправить некоторые ошибки ручной разметки! Общая информация о том, что такое лес из Словении и специфичный небольшой кусочек информации из Архангельской и Иркутской областей в сумме позволили получить такой результат. Для визуализации полученных результатов мы использовали бинарную маску и карту вероятностей.
Снимок
Разметка
Предсказание модели: бинарная маска
Предсказание модели: карта вероятностей
Казалось бы, все отлично — наша нейросеть уже видит лучше нас! Ученик превзошел учителя! Однако, на деле это не совсем так.
К сожалению, было много участков, где смысл результатов был весьма сомнителен. И на одном из примеров (кусок Иркутской области) можно продемонстрировать те трудности, с которыми мы столкнулись.
Предсказание модели
Ручная разметка
Явно выраженной разницы между ручной разметкой и предсказанием модели нет. Важно понимать, что на подложке здесь не снимок Sentinel с 10-метровым разрешением, а снимок из Яндекса, в который мы иногда подглядывали при разметке.
Мы снова принялись искать способы улучшения нашей модели.
Добавляем еще данных
В процессе последующих поисков были найдены новые данные с разметкой леса Европы (от Великобритании и Португалии до Финляндии и Турции). К тому же — на хвойный и лиственный. Показалось, что они в корне изменят всю дальнейшую работу.
Для работы мы выбрали достаточно обширные территории Финляндии и Норвегии, леса которых похожи на сибирские.
Однако, и попытка распознавания этих данных нейросетями или Random Forest как «лес/не лес», и попытка распознавания как «хвойный/лиственный», хоть и давали какие-то, порой осмысленные, результаты, но в итоге все они были отбракованы. Несмотря на весьма неплохое качество, подходящее разрешение и более высокие метрики моделей, данные снова оказались неприменимы для качественного решения поставленных задач в целевом регионе.
Из чего мы пришли к разочаровывающему для нас итогу — леса Европы не похожи на леса Сибири. Это объясняется и почвой, и прочими геоботаническими соображениями. После множественных просмотров мы заключили, что ни данные Словении, ни данные Финляндии и Норвегии не делают предсказания модели лучше. Модели, обученные только на нашей ручной разметке, давали гораздо более осмысленный результат.
Мы приняли решение отложить данные европейских стран и сделать дополнительную ручную разметку уже по Сибири. Для неё отбирались регионы таким образом, чтобы учесть самые типичные ошибки нашей модели. Была выдвинута гипотеза, что имеет смысл выделить наиболее очевидные участки хвойного и лиственного лесов для разметки. Потому что даже на снимке с разрешением 1 метр на пиксель далеко не всегда понятно, какой лес перед тобой — хвойный или лиственный.
Предсказания модели до использования точечной разметки
Предсказания модели после использования точечной разметки
Идея с точечной разметкой, которая указывает модели её типичные ошибки, в результате показала себя эффективной. Но очевидно также и то, что остаются зоны, на которых нет леса, но обе модели размечают его там.
Для того, чтобы решить эту проблему, мы использовали второй снимок. Как я и сказал в начале, вегетационные индексы наибольшую разницу между разными классами поверхности дают именно в динамике, и мы решили применить это способ.
Однако в свободном доступе есть только инструмент, который не позволяет загрузить сразу много кадров с приемлемым уровнем облачности за необходимые даты и требует большое время ожидания для выгрузки архивных данных, кроме того, сервис попросту иногда глючит. Разумеется, есть и продукт без подобных ограничений, но за определенный прайс, а мы, как вы помните, работаем только с открытыми данными.
Поэтому, мы не могли использовать длинный временной ряд снимков, однако, как выяснилось, всего одного дополнительного снимка за другой вегетационный период оказалось вполне достаточно. При одновременном использовании сразу двух снимков результаты действительно намного улучшились. На тех проблемных участках, которые раньше распознавались как лес, новая модель уже не делала подобных ошибок. Мы избавились от ложноположительных срабатываний довольно серьёзно.
После применения второго снимка результаты серьезно улучшились
На этом этапе мы стали использовать летние и осенние снимки, потому что наиболее заметная динамика в разных участках спектра происходит именно на этом переходе.
На ранних этапах у нас еще было предположение использовать зимние снимки для того, чтобы размечать лес на «хвойный/лиственный». Для этого предполагалось использовать не вегетационный, а снежный индекс и смотреть на снежное покрытие и соотносить его с маской леса.
После проверки ряда гипотез выделить надежный критерий, по которому было бы удобно размечать разные типы леса с использованием зимних снимков, не получилось. Во многом это связано с тем, что снег — явление не постоянное. Его может сдуть, он ложится не только на лес и не всегда лежит равномерно. Поэтому, от этой идеи мы решили пока что отказаться.
Снимок
Результаты предсказаний моделей
Еще один пример. На снимке желтым показаны результаты работы последней модели, фиолетовым — промежуточной, красным — то, что отбивалось как лес первой моделью. Как видите, тут тоже довольно серьёзные улучшения.
Промежуточные итоги
Наша работа ещё далека от завершения, но на данном этапе мы уже пришли к нескольким важным для нас выводам.
В настоящий момент мы работаем над внедрением в пайплайны данных цифровой модели рельефа (DEM) и ожидаем, что это даст нам большой прирост точности. Также мы улучшаем наши практики проведения экспериментов (логирование процесса обучения, документирование моделей, разработка библиотеки для удобной работы с различными форматами геоданных в контексте задач ML).
Обо всем, конечно же, расскажем в следующих материалах.
Задавайте вопросы и делитесь своим опытом в решении подобных задач — будет интересно продолжить обсуждение в комментариях.