Пример использования утилит GDAL для вычисления яркостных характеристик снимков (на примере выявления рубок по зимним данным Landsat)

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Эта страница является черновиком статьи.


Пример решения технических задач, возникающих при выявлении изменений отражательных свойств поверхности по разновременным материалам ДЗЗ

Одним из способов выявления изменений между сериями разновременных снимков является проведение различных математических операций между значениями яркости пикселей. Настоящая статья создана с целью помочь начинающим пользователям освоить основы техники данного приема с использованием свободного ПО и является опытом её автора. Для примера выбрана достаточно простая практическая задача и соответствующие ей исходные данные.

Используемое ПО

Командная строка Windows 7, Python27, GDAL 1.11.1

Исходные данные

Зимние безоблачные снимки Landsat 8 OLI панхроматический(_B8) канал одной сцены формат TIF, уровень обработки L1T в проекции UTM(отмечу что для моего региона в средней полосе Европейской части, удалось найти порядка по 1 такому снимку за весь период работы сенсора(с апреля 2013 г.) в периоде февраль-март за 2015 и 2016 г.).

Краткое описание действий

Участки вырубленные за период между датами съемок на более позднем снимки будут иметь значительно большие значения яркости, чем на раннем (естественно при условии наличия снежного покрова и не слишком большого временного интервала). Расчет разности яркости и последующая классификация результата по пороговому значению осуществляется с использованием утилиты gdal_calc.py.

Создание общей области покрытия

Для корректной работы утилиты gdal_calc.py с несколькими растрами необходимо, чтобы входящие растры имели одинаковый экстент. Получить сведения об экстенте можно используя утилиту gdalinfo.exe: 

gdalinfo rastr
  • здесь и далее приведен пример синтаксиса командной строки
  • rastr - полный путь к набору растровых данных

Результат работы утилиты, будет представлен в виде текста в командной строке: Вывод значений координат углов изображения(экстент) в командной строке

Необходимо определить максимальные значения для левой и нижней, и минимальные для правой и верхней границ используемых изображений. Обрезка растров производится с помощью утилиты gdalwarp.exe: 

gdalwarp.exe -te Xmin Ymin Xmax Ymax -tr 15 -15 inputRastr outputRastr
  • -te экстент обрезки
  • -tr размер ячейки по x,y
  • для gdalwarp и gdal_calc возможно использование флага -overwrite для перезаписи выхдногофайла, если он уже существует

Для получения правильного результата необходимо также, чтобы валидные значения яркости одного растра не попадали на значения NODATA (0) другого, т.е. получить область пересечения растров. Данную операцию необходимо провести для каждого из растров приведенных к одинаковому экстенту используя gdal_calc.py: 

gdal_calc.py -A inpRastr1 -B inpRastr2 --outfile=outRastr1 --calc="A*(A>0)*(B>0)" --NoDataValue=0
gdal_calc.py -A inpRastr1 -B inpRastr2 --outfile=outRastr2 --calc="B*(A>0)*(B>0)" --NoDataValue=0
  • gdal_calc.py(как и в случае с любым другим исполняемым файлом python)-сокращенное обозначение текста
[путь к файлу python.exe] [путь к gdal_calc.py]

,например: 

C:\Python27\python.exe "C:\Program Files\GDAL\gdal_calc.py"(пути с пробелами беруться в кавычки)
  • для значений 0 будет установлен флаг NODATA

Вычисление разности каналов

Теперь можно провести вычисление собственно разности каналов: 

gdal_calc -A earlyRast -B laterRast --outfile=outRast --calc="B.astype(int)-A.astype(int)" --type "Int32"
  • для выходных данных используется тип данных Int32-значения от -2 147 483 647 до +2 147 483 647
  • пиксели имевшие значение 0(NODATA) получают значение -2 147 483 647 с флагом NODATA, а значения 0 полученные в результате вычисления остаются валидными

Классификация

Для выявление непосредственно вырубок, как объектов, необходимо классифицировать полученный одноканальный растр разности на два класса- 1)пиксели с изменением яркости соответствующим вырубке и 2)все остальные: 

gdal_calc.py -A inpRast --outfile=outRast --calc="1*(A>value)" --type "Byte"
  • value - пороговое значение определяемое пользователем эмпирическим путем
  • тип выходных данных Byte (для экономии памяти)-значения от 0 до 255
  • пиксели имевшие значение -2 147 483 647(NODATA) получают значение 255 с флагом NODATA, пиксели, удовлетворяющие условию- значение 1, все остальные-0

Постобработка

Погрешности в геометрической привязке снимков, а также особенности снегонакопления на границе различных категорий земель могут привести к ошибочному отнесению к вырубкам участков вдоль линейных объектов(рек, дорог, лесополос) и по берегам водоемов. Уменьшить количество "нежелательных" пикселей, а также заполнить "пустоты" в полигонах можно с помощью фильтра- утилиты gdal_sieve.py:

Преобразовать классифицированные данные в векторный формат можно с помощью gdal_polygonize.py:

Замечания

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

Источник — https://wiki.gis-lab.info/index.php?title=Пример_использования_утилит_GDAL_для_вычисления_яркостных_характеристик_снимков_(на_примере_выявления_рубок_по_зимним_данным_Landsat)&oldid=24156