Пример использования утилит GDAL для вычисления яркостных характеристик снимков (на примере выявления рубок по зимним данным Landsat)
Пример решения технических задач, возникающих при выявлении изменений отражательных свойств поверхности по разновременным материалам ДЗЗ.
Одним из способов выявления изменений между сериями разновременных снимков является проведение различных математических операций между значениями яркости пикселей. Настоящая статья создана с целью помочь начинающим пользователям освоить основы техники данного приема с использованием свободного ПО и является личным опытом автора. Для примера (результат на рис. 1) выбрана достаточно простая практическая задача и соответствующие ей исходные данные.
Используемое ПО
ОС Windows 7, командная строка (CMD.EXE), Python 2.7.*, 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. Радиометрическую калибровку можно также осуществить с помощью gdal_calc.py, руководствуясь указанной выше статьей.
Создание общей области покрытия
Для корректной работы утилиты gdal_calc.py с несколькими растрами необходимо, чтобы входящие растры имели одинаковый охват (экстент). Получить сведения об охвате можно, используя утилиту gdalinfo:
gdalinfo rastr
- здесь и далее приведен пример синтаксиса командной строки
- rastr - полный путь к набору растровых данных
Результат работы утилиты будет представлен в виде вывода командной строки:
Необходимо выбрать максимальные значения для левой и нижней и минимальные для правой и верхней границ среди охватов всех используемых изображений. Обрезка растров производится с помощью утилиты gdalwarp:
gdalwarp.exe -te Xmin Ymin Xmax Ymax -tr Xres Yres inputRastr outputRastr
- после имени исполняемого файла синтаксис следующий: "-(название флага) (параметры через пробел)"; могут быть флаги без параметров
- флаг -te экстент обрезки (параметры Xmin Ymin Xmax Ymax, определенные выше)
- -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 для минимального охвата как указано выше или вычислить самостоятельно для полигона обрезки
- полигон обрезки и растры должны иметь одинаковую систему координат
Вычисление разности каналов
Теперь можно провести вычисление собственно разности каналов:
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
Постобработка
Погрешности в геометрической привязке снимков, а также особенности снегонакопления на границе различных категорий земель могут привести к ошибочному отнесению к вырубкам участков вдоль линейных объектов(рек, дорог, лесополос) и по берегам водоемов(рисунок 2).
Уменьшить количество "нежелательных" пикселей, а также заполнить "пустоты" в полигонах можно с помощью фильтра- утилиты 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 "Упростить геометрию".
Замечания
Следует заметить, что делянки, вырубленные незадолго до более поздней съемки определяются плохо, т.к. они покрыты порубочными остатками и разница в яркости может оказаться ниже установленного порогового значения.