Обработка разновременных снимков
Применение разновременных снимков Sentinel-2 для классификации культур и оценки неоднородности полей
Снимки спутников Sentinel-2 являются одним из наиболее оптимальных источников данных дистанционного зондаирования для сельского хозяйства, поскольку обладают следующими свойствами: • регулярной съемкой одной и той же территории; • достаточно высоким разрешением (10 метров); • большим количеством спектральных каналов, в том числе инфракрасного диапазона; • возможностью бесплатного получения.
Для того, чтобы работать с разновременными снимками в ArcGIS, их нужно выбрать и объединить в мозаику растров с поддержкой времени или в многомерный растр формата CRF.
Классификация культур по разновременным сценам Sentinel-2
Поля, которые находятся в обороте, проходят последовательные циклы посадки, роста и уборки урожая, которые в свою очередь достаточно хорошо дешифрируются при обработке разновременных сцен. Отсутствие существенной динамики растительности тоже является признаком, который позволяет отнести поле к определенному классу.
В качестве тестового участка была выбрана территория Узбекистана в районе города Вабкент, были отобраны сцены за май-октябрь 2022 года с шагом в 1 месяц (6 сцен), которые были объединены в один многомерный растр . Для каждой сцены был построен индекс NDVI.
Обучающая выборка была составлена на основе случайно размещенных точек, спектры которых были сгруппированы в 5 классов с помощью инструмента многофакторной кластеризации . В качестве спектра в данном случае выступает динамика индекса NDVI на протяжении 6 месяцев с мая по октябрь.
Кластеры с сходными показателями динамики индекса в течение 6 месяцв
Для классификации был применен новый алгоритм классификации растров на основе спектра , который появился в ArcGIS Pro 3.1. Основная идея алгоритма заключается в следующем: спектр пиксела можно представить как n-мерный вектор, где n это количество каналов изображения. Принадлежность пиксела к классу измеряется на основе расчета угла между вектором конкретного пиксела и эталонным вектором (спектром) обучающей выборки.
Результат классификации
Наземных (контрольных) данных не было, поэтому интерпретация динамики индекса NDVI основана на предположениях.
Характерные спектры классов:
Постепенный рост с апреля по август и уборка в сентябре
Уборка урожая в мае-июне и далее поле не используется до конца октября
Уборка двух урожаев за полгода, первый в мае-июне, второй в сентябре-октябре
В класс "сады" попали все остальные поля, которые не имеют четко выраженной сезонной динамики
Предварительные выводы, плюсы и возможные проблемы данного метода классификации
Динамика индекса NDVI (как вариант, можно использовать любую комбинацию каналов) дает возможность более точно разделять сельскохозяйственные культуры по сравнению с классификацией отдельного снимка. В качестве классификатора можно использовать самые разные инструменты, в т.ч. алгоритмы машинного обучения (Random Forest, SVM и другие). В ArcGIS также есть возможность использования специализированных нейросетей ( например, PSETAE ).
Для классификации по разновременных ДДЗ желательно подобрать безоблачные снимки за определенный период времени. Период времени подбирается в зависимости от задачи. Не факт, что классификатор, который обучен на определенной территории и на определенный год, можно будет потом использовать в другие годы для той же самой территории. Возможно, эту проблему смогут преодолеть нейросети архитектуры PSETAE, которые ориентируются не столько на конкретные значения пикселей, сколько на общую картину динамики показателей. В любом случае ArcGIS на данный момент содержит все необходимые инстурменты для экспериментов в этом направлении.
Статистическая обработка разновременных снимков для оценки неоднородности почвенного покрова
Пример базируется на попытке реализовать средствами ArcGIS идею расчета индекса неоднородности полей. Данный метод придуман Дмитрием Руховичем: https://istina.msu.ru/profile/RukhovichDI/ Здесь можно посмотреть его доклад на одной из конференций: https://www.youtube.com/watch?v=J5ln2L7QLfY&t=6972s
Для расчета этого показателя нужно сделать следующее: 1) Строится индекс NDVI для набора разновременных снимков. 2) Каждое поле внутри себя разбивается на две или три зоны (высокие, средние и низкие значения). 3) Подсчитывается суммарная статистика за период наблюдений.
Каждое поле делится на три зоны методом equal area (перцентили)
Если некий участок поля имеет низкие значения индекса (относительно других участков того же поля) много месяцев или много лет подряд, значит эти плохие значения не случайны. Если просуммировать эту статистику за 6 месяцев, то можно увидеть следующую картину.
Красный цвет - участки с систематически плохими значениями индекса NDVI
На примере видно, что участки полей с постоянно низким (относительно других участков поля) индексом NDVI часто располагаются по краям, некоторые пятна переходят с одного поля на другое, т.е. эти участки связаны не с посаженной в данный год с/х культурой и динамикой ее роста, а с какими-то постоянными негативными факторами. Это может быть эррозия, влияние рельефа, проблемы с поливом и так далее. Причину можно установить только на месте или путем сопоставления с другими пространственными слоями. Такая статистика в теории может помочь найти проблемные места и внести корректировки в количество вносимых удобрений. "Плохие" участки как правило хуже отзываются на высокие дозы удобрений и начиная с определенного порога, применение удобрений может стать нерентабельным, т.е. есть возможность сокращения затрат. Есть интересная статья Дмитрия Руховича и коллектива авторов, посвященная этой теме .
Как рассчитать подобную статистику в ArcGIS Pro? Можно вырезать из снимка отдельное поле, обработать его инструментом Slice модуля Spatial Analyst и потом соединить все эти поля в один итоговый растр, но это неудобно. В ArcGIS Pro 3.1 эту задачу можно быстрое решить с помощью скрипта Python. Создаем объект Raster Collection на базе многомерного растра или директории с растрами, применяем функцию map, которая обрабатывает все снимки набора, конвертируем Raster Collection в многомерный растр и сохраняем результат на диск в формате CRF. Для рассчета статистика по каждому полю потребуется растр с маской полей, где каждое поле имеет свое уникальное целочисленное значение яркости.
# временные растры будут храниться в памяти, а не в БГД по умолчанию
arcpy.env.scratchWorkspace = 'memory'
arcpy.env.overwriteOutput = True
# Инициализация raster collection, поле с временем должно называться
rc = arcpy.ia.RasterCollection('путь к многовременному растру или директории с растрами')
# Название поля с временным измерением, если raster collection на базе директории с снимками,
# то нужно это поле создать и вычислить в нем даты для каждого снимка
time_field = 'StdTime'
# Растр с полигонами (полями), внутри которых будет рассчитана статистика
# Каждый полигон должно иметь свое уникальное занчение в атрибуте Value
zone_raster = arcpy.Raster('путь к растру с зонами')
def zonal_3_parts(item):
"""Функция разбивает растровый полигон на три зоны методом equal area"""
item['Raster'].save('memory/temp')
zonal_stat33 = arcpy.ia.ZonalStatistics(zone_raster, 'Value', 'memory/temp', 'PERCENTILE',
percentile_value=33,
process_as_multidimensional='CURRENT_SLICE')
zonal_stat66 = arcpy.ia.ZonalStatistics(zone_raster, 'Value', 'memory/temp', 'PERCENTILE',
percentile_value=66,
process_as_multidimensional='CURRENT_SLICE')
good = arcpy.ia.Con('memory/temp' > zonal_stat66, 1, 0)
medium = arcpy.ia.Con(('memory/temp' >= zonal_stat33) & ('memory/temp' <= zonal_stat66), 2, 0)
bad = arcpy.ia.Con('memory/temp' < zonal_stat33, 3, 0)
# подсчет суммарной статистики
result = arcpy.ia.Sum([good, medium, bad])
return {"raster": result, "name": "good_medium_bad", "StdTime":item[time_field]}
# Применение функции zonal_3_parts к набору растров
rc_sliced = rc.map(zonal_3_parts)
# Конвертация raster collection в многомерный растр с измерением StdTime
output = rc_sliced.toMultidimensionalRaster('name', ['StdTime'])
# Сохранение результата на диск
output.save('путь к растру с расширением .crf')Если у вас возникли вопросы, присылайте их на адреса ykopin@esri-cis.com или market@esri-cis.com.
