Comments 41
Использование Sen2Cor осложняется тем, что это как конь сферический в вакууме. Сцены Sentinel-2 нарезаны на тайлы, если вы хотите сделать атмосферную коррекцию для двух и более тайлов, то готовьтесь к Different results with tiles from one scene. Даже если тайлы сняты одним КА и в одно время.
Я ещё долго искал где-нибудь, какой именно алгоритм обработки реализован в Sen2Cor (там же разные есть — DOS и 6S, например), но так и не нашёл.
По ссылке в моем комментарии есть информация о проекте MAJA от ФранцКосмоса, похоже, что там более соответствующие друг-другу результаты получаются. MAJA бесплатна для некоммерческого использования, но сразу скажу, что сам не пробовал. Ну и задачи атмосферной коррекции могут быть разными, если нужно прямо вот суперточный BOA, то это одно, а если нужно, например, отслеживать динамику NDVI, то главное, чтобы результаты отражали реальный тренд.
Хотелось бы увидеть в статье не только описание того, на какие кнопки нужно нажимать, но и того зачем это нужно делать. Какую новую информацию я получу со снимка, применив тот или иной шаг коррекции? Как это поможет мне использовать этот снимок? Пример в шапке не показателен, т.к. похожий результат достигается десятком кликов в фотошопе.
Использую космоснимки для анализа рельефа при подготовке пешеходных маршрутов в малопосещаемых районах. Было бы интересно почитать о применении подобного ПО для анализа снимков, а не просто о том, какие кнопки тыкать.
Пример в шапке действительно не очень показателен, потому что это просто картинка, к тому же полноцветная. Если сравнивать отдельные каналы по отражению в заданных точках, то разница будет видна, и будет понятно, что это не Фотошоп. Скорректированные каналы можно, например, использовать для расчёта вегетационных индексов — считается, что в результате атмосферной коррекции значения должны больше соответствовать действительности. То есть если мы посчитаем индекс до коррекции и индекс после коррекции и сравним, например, с результатами полевых измерений, то в теории индекс после коррекции должен быть ближе к результатам полевых измерений, чем индекс до коррекции.
На самом деле, разрешение 10 м/пиксель — это довольно неплохая точность, если сравнивать с другими доступными проектами. Google Maps, например, предлагает снимки 15-30 м/пиксель, и то на отдельные районы.
DigitalGlobe как-то попытались использовать, но для некоммерческой группы там запредельные цены. Дальше запроса стоимости дело не пошло.
OpenTopoMap для построения рельефа использует данные SRTM с разрешением ~90м, что не обеспечивает достаточной точности. К тому же, в приполярных областях данные SRTM вообще отсутствуют. Неплохие результаты даёт ALOS, но только для некоторых районов. Большие области или не отсняты, или данные получены со значительной погрешностью. Существенные проблемы приносит сглаживание участков с резким изменением крутизны склонов. То, что выглядит как 40° склон, может запросто оказаться двумя участками по 20° и 70°. Такая стенка наверняка станет неприятным сюрпризом.
Ещё используем снимки Landsat 7,8 в видимом диапазоне — для выявления форм рельефа по характерным признакам, освещённости и т.д. — и иных — например, различить в ИК лёд и снег, которые в видимом диапазоне выглядят одинаковым белым пятном. Всё это, естественно, тоже вносит кучу погрешностей при ортотрансформации, склейке снимков и т.д. В итоге приходится полагаться только на невеликий накопленный опыт, чтобы свести всё воедино и попытаться угадать что же там такое на самом деле, поэтому и интересно получить максимальное количество информации с таких снимков.
ArcGIS World Imagery
www.arcgis.com/home/item.html?id=10df2279f9684e4a9f6a7f08febac2a9
Например, на горы Северо-Восточного Алтая там разрешение 1.2 метра на пиксель, а точность привязки снимков 8.5 м (данные показываются после клика левой кнопкой на веб-карте). Я пользуюсь через SasPlanet.
Мне кажется вы путаете тональную балансировку снимков и атмосферную коррекцию — Atmospheric Correction. Для первой задачи подходит растровый редактор, вторая — очень сложная научная задача.
Я иногда занимаюсь ручной дешифровкой степных пожаров в Даурии. Раньше я выкачивал 150 снимков на десктоп, и настраивал цвета нажимая кучу кнопок в QGIS. А теперь я за городом, и не могу выкачать с сервера 150 16-битных снимков. Я могу сконвертировать их в RGB, но перед этим надо с ними сделать что-то вроде Auto Colors.
Честно говоря не знаю. Можно посмотреть или OrfeoToolbox или, кажется, видел это в ScanEx Image Processor. Orfeo хорош тем, что там есть API, в т.ч. для python.
Но общая проблема в том, что из-за облачности очень тяжело найти тайлы примерно совпадающие по периоду съемки. Это приведет к тому, что с одного края у вас будут тайлы июльские, а с другого майские, а такие, сильно не совпадающие по времени тайлы, практически невозможно отбалансировать, ну только если это не пустыня. Так что для моих задач тональная балансировка не актуальна.
Хотя, сам долго жил в Даурии, солнца там очень много, может для Даурии и актуально делать бесшовные тонально сбалансированные мозаики.
Вот мозаика по снимкам за 3 года — все отлично с цветами, на мой взгляд. Облачные пиксели исключены из обработки фильтрацией. Делается бесплатно в Google Earth Engine.
Ну средняя за 3 года конечно получится, даже по СЗФО, внутри одного года с высокой вероятностью не получится, точнее зависит от региона. В Даурии получится, там под 300 солнечных дней в году, в других регионах по 60 солнечных дней в году.
А Google Earth Engine потом в как WMTS у себя в проекте можно использовать?
Атмосферная коррекция «обывателю» не нужна вообще. Эта процедура нужна только тогда, когда вы собираетесь производить фотограмметрические изменения.
A, ну если вам нужна лицензионная чистота — Sentinel вам подойдет, потому что на снимках QuickBird/WorldView/GeoEye вы лицензионно-чистую схему не нарисуете.
Слышали ли вы что-нибудь про воспроизводимость исследований?.. Есть же готовый датасет Sentinel-2 Surface Reflectance (c атмосферной коррекцией, коррекцией угла съемки, положения Солнца,..) https://developers.google.com/earth-engine/datasets/catalog/COPERNICUS_S2_SR Берите данные отсюда и ваша работа будет воспроизводима — может, вам же и пригодится в дальнейшем.
P.S. А вообще лучше бы обойтись минимумом коррекций — например, посмотрите TOA (Top of Atmosphere Reflectance). Surface Reflectance снимков с перекрытием заметно отличается, что физически некорректно, так что на практике необходимо использовать серии снимков для создания композитного изображения, при этом фильтровать по маске облачности (приведенные в статье снимки явно очень облачные). Вот как пример, я как раз сегодня интересную статью про композит Sentinel-2 TOA на Google Earth Engine нашел и воспроизвел: https://github.com/mobigroup/gis-snippets/tree/master/GEE
Dataset Availability
2017-03-28T00:00:00 — Present
Как указано — все снимки с весны 2017 доступны. Сравните с «Sentinel-2B is a European optical imaging satellite that was launched on 7 March 2017.» Спутник Sentinel 2A с лета 2015 летает, то есть примерно полтора года архива по нему не обработаны. В любом случае, для существенно большей истории наблюдений смотрите Landsat 8 или предыдущие — по ним тоже есть датасеты TOA и SR.
А не в курсе ли вы на счет подвижек по гармонизации данных Landsat и Sentinel? Проблема в том, что в СЗФО всегда очень облачно и приходится использовать часть данных Landsat, а часть Sentinel-2. Из-за отсутствия гармонизации, что TOA, что BOA будут разными. А если у вас поле, у которого часть на S2, а часть на L8, то получается полная жесть, про анализ NDVI можно забыть.
Помню, что где-то мелькал проект по гармонизации этих данных, но как-то перестал следить и где-то потерялось все.
Разные спутники с разным оборудованием и разными спектральными частотами каналов и самими каналами — универсально никак нельзя объединить. А вот для конкретной задачи и территории — можно сравнить значения NDVI для перерывающейся области разных снимков и объединить их. Кстати, вы точно не отдельными сценами оперируете, а набором фильтрованных от облаков сцен с медианной фильтрацией значений каждого пиксела? Для Sentinel-2 известно:
With two satellites all regular areas indicated above will be revisited every five days under the same viewing conditions.
Даже если вам нужно как можно чаще анализ делать — можно маскировать облачные пиксели, получая частичное покрытие NDVI.
Спектральные диапазоны каналов в значительной степени пересекаются у разных КА, поэтому можно попытаться найти общий знаменатель. Вопрос в том, насколько это применимо. Думаю для задач анализа тренда, например NDVI, — более чем. Тут попытались Harmonized Landsat Sentinel-2, признаюсь, однако, что результат не проверял.
Кстати, вы точно не отдельными сценами оперируете, а набором фильтрованных от облаков сцен с медианной фильтрацией значений каждого пиксела?
Вы про композитный результат на основе серии снимков (я так понял revisit time вы для этого привели)? Нет, я так не делаю. Причин несколько:
- Возьмите некоторые области СЗФО, там облачность такая, что будет один безоблачный снимок в мае, один в июле, а следующий в декабре-январе. И revisit time даже чаще чем 5 дней, бывает даже через день, но увы… облака.
- На данный момент мои интересы в стороне от сельского хозяйства, мне не особо сейчас нужны атмосферные коррекции и фильтрации. Меня больше интересует лес и обработка радаров Sentinel-1.
Но да, про фильтрацию я согласен, если есть набор данных, то делать можно.
Все же, гармонизация Landsat и Sentinel важный вопрос в районах с частой облачностью. Ситуация когда полполя в Landsat, а полполя в Sentinel — просто рядовая.
Возьмите некоторые области СЗФО, там облачность такая, что будет один безоблачный снимок...
Не нужно искать безоблачные снимки! Возьмите облачные и маскируйте облачные пиксели — по соответствующему спектральному каналу и, если хотите, дополнительные маски добавьте. Лично я часто сталкиваюсь с около экваториальными территориями — там безоблачных снимков не бывает вообще, так что только фильтрация безоблачных пикселов с разных сцен спасает (и последующая медианная фильтрация попиксельно).
Ситуация когда полполя в Landsat, а полполя в Sentinel — просто рядовая.
Так и анализируйте раздельно — используя потом разные пороговые коэффициенты для NDVI по Landsat и Sentinel-2. Точнее будет, чем пытаться привести результат к единому NDVI и по нему анализировать.
Возьмите облачные и маскируйте облачные пиксели — по соответствующему спектральному каналу и, если хотите, дополнительные маски добавьте. Лично я часто сталкиваюсь с около экваториальными территориями — там безоблачных снимков не бывает вообще, так что только фильтрация безоблачных пикселов с разных сцен спасает (и последующая медианная фильтрация попиксельно).
Извините, но никак в толк не возьму, что маскировать если облачность 100%? Наверное на экваторе нет смены времен года и такой подход оправдан. Но если у вас один приемлемый снимок в начале мая, потом облачность непролазная, потом один в июле, следующий, может быть осенью, а потом зимой в морозы и больше ничего. Правильно ли я понял, что надо майские пиксели смешивать с июльскими? Что даст маскирование? Ну уберется 75% снимка, как это поможет при условии, что больше на эту территорию снимков нет за интересующий период нет? Бывает необходимость анализа внутри определенного периода, например только во время цветения определенного растения, так запросто в цветение может не быть никаких снимков вообще, только 100% облачность.
Ну и облака бывают мелкими барашками и солнце низко, но яркое, тогда будут четкие тени от облаков на земле, которые пока непонятно как автоматически находить и маскировать. В СЗФО не получается даже подобрать ближайший Landsat 8 под "гребенку" Landsat 7.
Могли бы вы привести пример того о чем вы говорите, но не за 3 года, а внутри сезона, желательно лета, например, для средней полосы европейской части РФ?
Возьмем скрипт примера из GEE и добавим фильтр облачности (не более 20% облаков) и фильтр по дате (май 2019) и отфильтруем по заданной территории (прямоугольник 29.5, 59.5, 31.5, 60.5 — Санкт-Петербург примерно в центре).
Получаем композитное изображение с полным покрытием:
На мой взгляд, тут есть некоторые небольшие артефакты, их можно вычистить более аккуратной фильтрацией облаков. Как видим, ситуация на деле далека от "не быть никаких снимков вообще".
Вот и сам скрипт GEE:
// This example uses the Sentinel-2 QA band to cloud mask
// the collection. The Sentinel-2 cloud flags are less
// selective, so the collection is also pre-filtered by the
// CLOUDY_PIXEL_PERCENTAGE flag, to use only relatively
// cloud-free granule.
// Function to mask clouds using the Sentinel-2 QA band.
function maskS2clouds(image) {
var qa = image.select('QA60')
// Bits 10 and 11 are clouds and cirrus, respectively.
var cloudBitMask = 1 << 10;
var cirrusBitMask = 1 << 11;
// Both flags should be set to zero, indicating clear conditions.
var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
qa.bitwiseAnd(cirrusBitMask).eq(0))
// Return the masked and scaled data, without the QA bands.
return image.updateMask(mask).divide(10000)
.select("B.*")
.copyProperties(image, ["system:time_start"])
}
// Set map area
var area = ee.Geometry.Rectangle(29.5, 59.5, 31.5, 60.5);
Map.centerObject(area, 8);
// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var images = ee.ImageCollection('COPERNICUS/S2_SR')
.filterDate('2019-05-01', '2019-05-31')
// Pre-filter to get less cloudy granules.
.filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
.map(function(image) { return image.clip(area); });
var composite = images
.map(maskS2clouds)
.median()
// Display the results.
Map.addLayer(composite, {bands: ['B4', 'B3', 'B2'], min: 0, max: 0.3}, 'RGB')Да, но май не интересен:
- почти всегда в начале мая есть немалое количество солнечных дней;
- в мае только начинается рост растений, для с/х не подойдет.
Возьмите, например, северную часть Рыбинского водохранилища и Череповец, тайл 37VDF, и возьмите июль. Скорее всего ничего не получится, а даже если получится, то будут четкие тени от маскированных облаков — барашков. Это в том случае если облака замаскируются хорошо. Перистые регулярно пролетают через все фильтры и маски.
Хорошо, смотрим на ту же территорию июль:
И здесь явно речь не идет про озвученную вами 100% облачность! Притом, что мы использовали простейший скрипт обработки.
Upd. Дополним еще данными Sentinel-2
Тоже взят скрипт из примеров GEE для нужного района и времени с ограничением не более 20% облачности на сцене.
А теперь возьмем минимальное значение вместо медианы — облака это точно отфильтрует:
Итого — ваша задача превратилась из не решаемой в элементарную, верно?:) Последней картинки вам может вполне хватить, трудозатраты на это — запуск официального примера и замена функции фильтрации медианы на минимум. Куда уж проще.
Еще можно посмотреть другие спутники, можно подобрать более подходящий алгоритм фильтрации облаков (обычно для летних снимков оценивают наиболее зеленый и наиболее синий пиксели, чтобы выбрать наилучшее значение для суши и воды, вместо медианной и минимальной фильтрации).
:-) собственно мы обсуждаем метод, но не определили цель, что немного наоборот :-). Конечно, ваш метод применим, особенно для получения визуальных мозаик, собственно я с этим никогда не спорил. Да, метод применим даже для более глубокого анализа когда есть серия близкорасположенных по времени снимков. Но:
- у вас в июле много облачности, не факт, что интересующий вас объект будет виден под облачностью
- попробуйте собрать мозаику на экстент тайла 37VDF за июль.
Собственно весь вопрос: какую цель мы хотим достичь? Ну а от цели можно и метод выбрать и да, ваш метод важен и нужен и вполне может быть применен. Но есть области, где он работать просто не будет и тогда нужны другие методы. Собственно к радарам обратился именно из-за нехватки оптических данных. За всю осень может не быть ничего кроме "молока"!
Тем же методом и радарные данные обрабатываются :) Экстент назовите, который вас интересует, посмотрим — широта_мин, долгота_мин, широта_макс, долгота_макс.
Экстент (широта_мин, долгота_мин, широта_макс, долгота_макс) 58.5407855054331 37.2814375455109 59.5382353517154 39.1725811832316. За июль, без облаков, без резких теней от "барашков", пожалуйста :-)
Ну с радарами отдельная пестня, там погода влияет, но не как в случае с облаками. Если вам надо найти вырубку, и вы даже знаете квартал, то далеко не всегда она будет видна на радаре. Может быть видна сегодня, а через 3 дня не видна. И вырубка быстроизменяющийся объект, для них подходят или одна сцена или ближайшие друг к другу сцены.
Пожалуйте:
А вырубки вы, надеюсь, не глазами ищете? Сделайте оконный спектральный анализ и свертку с гауссовой функцией — увидите, что получится. Понимаете, на определенном масштабе кроны достаточно симметричны (круговая симметрия) и это легко выделить.
Извините, но вы широту с долготой перепутали, в Ашхабаде наверное все хорошо с облачностью :-) Вы же сами попросили
Экстент назовите, который вас интересует, посмотрим — широта_мин, долгота_мин, широта_макс, долгота_макс.
Я ответил
Экстент (широта_мин, долгота_мин, широта_макс, долгота_макс) 58.5407855054331 37.2814375455109 59.5382353517154 39.1725811832316.
А я еще подумал почему долгота с широтой в экстенте наоборот находятся :-)
Хех, ок. А я еще подумал, что это вам в теплые края захотелось ;) Смотрим Sentinel-2 SR:
Это вовсе не 100% облачность, согласны?
Это июль одного года? Какого? Визуально выглядит хорошо.
Да, итоговая облачность далеко не 100% но она явно есть. Вопрос в том, что здесь явно будут присутствовать резкие тени от барашков и сервис на таких данных опасно построить. Если останется даже 20% облачности, то это означает, что из 10 заказов 2 выполнить не удастся, и заранее неизвестно будет облачно или нет. В остальном вопросов нет, для визуального представления все отлично.
Июль 2019 года, в скрипте (см. выше) месяц поменял. Резкие тени фильтруются тоже, это не проблема. Для своих целей я еще удаляю шумовую компоненту пространственного спектра (масштаб одного пикселя) — кажется крохоборством, но это существенно, и полученный композит отлично годится для дальнейшей обработки.
А как вы находите резкие тени? Вы их потом тоже маскируете и меняете значения?
Вы такие композиты классифицируете или у вас только для визуальных целей? Если классифицируете то насколько точно у вас получается находить требуемые объекты?
Сами понимаете, не стоит вместо медианы брать минимальное значение. Если данных мало — можно первый квартиль взять. Если совсем мало данных — анализируйте в два прохода, на первом оцениваете гистограмму и находите валидный интервал значений, на втором фильтруете по нему. Есть и отдельные алгоритмы фильтрации тени, но если у вас площади небольшие — подход на основе одной гистограммы для всей площади (точнее, одна гистограмма для суши и другая для воды, если нужно) работает отлично. Весь композит я использую только для визуализации, а работаю с данными по каналам — выделяю высокочастотные компоненты для детализации рельефа и потом гравики, что мне нужно для решения обратной задачи гравики и построения 3д модели плотности среды для геологоразведочных проектов. На Linkedin у меня много про это постов и статей, от теории до практики.
Сделайте оконный спектральный анализ и свертку с гауссовой функцией — увидите, что получится.
Если не сложно, можете ссылку на метод? Это к какой поляризации относится? У меня подозрение, что вырубка, при определенных условиях ничем не отличается от невырубки на радаре. Отражение что у вырубки, что у окружающего леса одинаковое.
Метод имеет кучу названий — кольцевое преобразование Хафа, преобразование Радона, фокальный анализ,… Вот можете на Linkedin посты глянуть:
https://www.linkedin.com/posts/alexey-pechnikov_circlet-transform-robust-activity-6634335533017456640-u4BP/
https://www.linkedin.com/posts/alexey-pechnikov_geophysical-applications-circular-activity-6634334529068859392-Pt7X/
https://www.linkedin.com/posts/alexey-pechnikov_radon-transform-methodology-activity-6629799196331470848-7-X-/
Это именно выделение кольцевых форм, так что можно отличить одно большое дерево от более мелкой растительности. Ну а фильтрация в спектральной области понятно зачем нужна — чтобы работать только с объектами интересующих масштабов.
Как обрабатывать спутниковые снимки с помощью Sen2Cor