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

{{Статья|Черновик}}
{{Аннотация|Пример решения технических задач, возникающих при выявлении изменений отражательных свойств поверхности по разновременным материалам ДЗЗ.}}

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

[[Файл:compare2015_2016.jpg|центр|50|frame|Рис.1 Слева фрагмент сцены LC81810202015048LGN00(февраль 2015 г.,см. описание исходных данных),справа тот же участок на сцене LC81810202016083LGN00(март 2016 г.), красный контур-полигональный слой участков леса вырубленных за период между сценами, полученный с помощью автоматизированного алгоритма, рассмотренного в данной статье.]]
==Используемое ПО==
ОС 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 г.).

==Краткое описание действий==
Участки, вырубленные за период между датами съемок, на более позднем снимке будут иметь '''значительно большие значения яркости''', чем на раннем (при условии наличия снежного покрова и не слишком большого временного интервала). Расчет разности яркости и последующая классификация результата по пороговому значению осуществляется с использованием утилиты [http://www.gdal.org/gdal_calc.html gdal_calc.py]. Естественно, в результат попадут все пиксели, значительно увеличившие свою яркость по разным причинам, но '''в подавляющем большинстве''' (если речь идет о лесной зоне) это будут именно изменения, связанные с исчезновением лесного покрова (в частности, сплошные рубки, приводящие к полному исчезновению древостоя).

==Радиометрическая калибровка==
Вопросы радиометрической калибровки достаточно полно рассмотрены в статье http://gis-lab.info/qa/landsat-data-correction.html, учитывая особенности исходных данных выбранных в качестве примера, калибровку проводить не обязательно(если не нужна разность в физических единицах), поэтому все команды приведены для расчетов в DN. Радиометрическую калибровку можно также осуществить с помощью gdal_calc.py, руководствуясь указанной выше статьей.

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

<pre>gdalinfo rastr</pre>

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

Результат работы утилиты будет представлен в виде вывода командной строки:
[[Файл:Ginfextent.jpg|лево|300|Вывод значений координат углов изображения (охват) в командной строке]]

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

<pre>gdalwarp.exe -te Xmin Ymin Xmax Ymax -tr Xres Yres inputRastr outputRastr</pre>
*после имени исполняемого файла синтаксис следующий: "-(название флага) (параметры через пробел)"; могут быть флаги без параметров 
*флаг -te экстент обрезки (параметры Xmin Ymin Xmax Ymax, определенные выше)
*-tr размер ячейки по x,y (для указанных данных 15 -15)
*для gdalwarp и gdal_calc возможно использование флага -overwrite для перезаписи выходного файла, если он уже существует

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

<pre>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</pre>
*для значений 0 будет установлен флаг NODATA

''Замечание:''
Для запуска утилит GDAL может понадобится указание полного пути к интерпретатору Python, например:
<pre>C:\Python27\python.exe "C:\Program Files\GDAL\gdal_calc.py" (пути с пробелами берутся в кавычки)</pre>

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

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

Во избежание пересчета пикселей при трансформации выходной растр необходимо замкнуть на сетку Landsat, используя параметры gdalwarp ''-te'' и ''-tr''. Можно использовать значения ''-te'' для минимального охвата, как указано выше, или вычислить самостоятельно для полигона обрезки.

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

==Вычисление разности каналов==
Теперь можно провести вычисление собственно разности каналов:
<pre>gdal_calc -A earlyRast -B laterRast --outfile=outRast --calc="B.astype(int)-A.astype(int)" --type "Int32"</pre>

* Преобразование входного типа данных (UInt16) в Int для вычислений необходимо для получения отрицательных значений.
* Для выходных данных используется тип данных Int32 (значения от -2 147 483 647 до +2 147 483 647).
* Пиксели, имевшие значение 0 (NODATA), получают значение -2 147 483 647 с флагом NODATA, а значения 0, полученные в результате вычисления, остаются корректными.

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

<pre>gdal_calc.py -A inpRast --outfile=outRast --calc="1*(A>value)" --type "Byte"</pre>

*value - пороговое значение определяемое пользователем эмпирическим путем
*тип выходных данных Byte (для экономии памяти)-значения от 0 до 255
*пиксели имевшие значение -2 147 483 647(NODATA) получают значение 255 с флагом NODATA, пиксели, удовлетворяющие условию- значение 1, все остальные-0

==Постобработка==
Погрешности в геометрической привязке снимков, а также особенности снегонакопления на границе различных категорий земель могут привести к ошибочному отнесению к вырубкам участков вдоль линейных объектов(рек, дорог, лесополос) и по берегам водоемов(рисунок 2).
[[Файл:compare2015_2016_lineError.jpg|центр|100|frame|Рис.2 Ошибочные результаты классификации(полигональный слой с красным контуром) вызванные геометрическими погрешностями снимков(слева 2015 г.,справа 2016 г.)]]

Уменьшить количество "нежелательных" пикселей, а также заполнить "пустоты" в полигонах можно с помощью фильтра- утилиты gdal_sieve.py:
<pre>gdal_sieve.py -st value -4 inpRast outRast</pre>
*value-целое число, растровые полигоны, меньше указанного размера в пикселях, будут удалены
*-4 -флаг указывающий на число анализируемых соседних ячеек при определении связности
*флаг NODATA не наследуется от входящего файла

Преобразовать классифицированные данные в векторный формат можно с помощью gdal_polygonize.py, но прежде, для того чтобы не создавать полигоны для значений отнесенных к не интересующему классу, удобно будет установить флаг NODATA дополнительно на значения 0 и 255 поэтапно:
<pre>gdalwarp.exe -dstnodata value inpRast outRast</pre>
*value-значение NODATA (сначала 0, затем 255)

собственно конвертация в вектор:
<pre>gdal_polygonize.py inpRast -f "ESRI Shapefile" outShp layer</pre>
*outShp-полный путь к выходному shp-файлу
*layer-имя слоя(на усмотрение пользователя)
Векторные данные, иллюстрирующие результат были получены дополнительной обработкой вектора инструментом геообработки QGIS "Упростить геометрию".

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

== См. также ==
* [http://gis-lab.info/qa/gdal-examples.html Примеры использования инструментов GDAL]
* [http://www.gdal.org/pages.html Утилиты  GDAL]