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

Эта страница является черновиком статьи.

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

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

LC81810202015048LGN00----------------------------------------------------------------------------------------------------------------------------------------------------------------------------LC81810202016083LGN00

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

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

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

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

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

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

Радиометрическая калибровка

Вопросы радиометрической калибровки достаточно полно рассмотрены в статье http://gis-lab.info/qa/landsat-data-correction.html, учитывая особенности исходных данных выбранных в качестве примера, калибровку проводить не обязательно(если не нужна разность в физических еденицах), поэтому все команды приведены для расчетов в DN. В случае использования данных разных сенсоров можно использовать готовые данные top-of-atmosfere reflectance (TOA) с EathExploer, если их нет(например для 8 канала), то можно вычислить с помощью gdal_calc.py, руководствуясь приведенной статьей. Вычисления можно проводить с исходными данными, или после приведения к общему покрытию, при этом необходимо отслеживать значения NODATA, например для данных LANDSAT 4-7 CLIMATE DATA RECORD (CDR) SURFACE REFLECTANCE это -9999.

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

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

gdalinfo rastr

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

Результат работы утилиты, будет представлен в виде текста в командной строке:

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

gdalwarp.exe -te Xmin Ymin Xmax Ymax -tr Xres Yres inputRastr outputRastr

• -te экстент обрезки
• -tr размер ячейки по x,y (для указанных данных 15 -15)
• для 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

Если имеется полигон обрезки, включающий только валидные значения во всех изображениях (например лист разграфки WRS-2) можно сразу обрезать исходные растры:

gdalwarp.exe -te Xmin Ymin Xmax Ymax -tr Xres Yres -dstnodata 0 -cutline shpfile inpRast outRast

• во избежании пересчета пикселей при трансформации выходной растр необходимо замкнуть на сетку Landsat используя параметры -te и -tr,

можно использовать значения -te для минимального экстента как указано выше или, если вычислить самостоятельно для полигона обрезки ближайшие значения x(y)min+res*n, где res-размер ячейки, n- целое число

• полигон обрезки и растры должны иметь одинаковую систему координат

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

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

gdal_calc -A earlyRast -B laterRast --outfile=outRast --calc="B.astype(int)-A.astype(int)" --type "Int32"

• преобразование входного типа данных (UInt16) в Int для вычислений необходимо для получения отрицательных значений
• для выходных данных используется тип данных 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

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

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

LC81810202015048LGN00----------------------------------------------------------------------------------------------------------------------------------------------------------------------------LC81810202016083LGN00

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

gdal_sieve.py -st value -4 inpRast outRast

• value-целое число, растровые полигоны, меньше указанного размера в пикселях, будут удалены
• -4 -флаг указывающий на число анализируемых соседних ячеек при определении связности
• флаг NODATA не наследуется от входящего файла

Преобразовать классифицированные данные в векторный формат можно с помощью gdal_polygonize.py, но прежде, для того чтобы не создавать полигоны для значений отнесенных к неинтересующему классу, удобно будет установвить флаг NODATA допонительно на значения 0 и 255 поэтапно:

gdalwarp.exe -dstnodata value inpRast outRast

• value-значение NODATA (сначала 0, затем 255)

,собственна конвертация в вектор:

gdal_polygonize.py inpRast -f "ESRI Shapefile" outShp layer

• outShp-полный путь к выходному shp-файлу
• layer-имя слоя(на усмотрение пользователя)

Векторные данные, иллюстрирующие результат были получены дополнительной обработкой вектора инструментом геообработки QGIS "Упростить геометрию".

Замечания

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

См. также

• Примеры использования инструментов GDAL
• Утилиты GDAL