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

Материал из GIS-Lab
Перейти к навигации Перейти к поиску
Нет описания правки
 
(не показано 15 промежуточных версий 2 участников)
Строка 1: Строка 1:
{{Статья|Черновик}}
{{Статья|Опубликована|gdal_spectral_characteristics}}
{{Аннотация|Пример решения технических задач, возникающих при выявлении изменений отражательных свойств поверхности по разновременным материалам ДЗЗ}}
{{Аннотация|Пример решения технических задач, возникающих при выявлении изменений отражательных свойств поверхности по разновременным материалам ДЗЗ.}}


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


[[Файл:compare2015_2016.jpg|центр|50|frame|LC81810202015048LGN00----------------------------------------------------------------------------------------------------------------------------------------------------------------------------LC81810202016083LGN00]]
[[Файл:compare2015_2016.jpg|600px|thumb|center|Рис.1 Слева фрагмент сцены LC81810202015048LGN00(февраль 2015 г., см. описание исходных данных), справа тот же участок на сцене LC81810202016083LGN00(март 2016 г.), красный контур-полигональный слой участков леса вырубленных за период между сценами, полученный с помощью автоматизированного алгоритма, рассмотренного в данной статье.]]
==Используемое ПО==
==Используемое ПО==
Командная строка Windows 7, Python27, GDAL 1.11.1? QGIS 2.14.2
ОС Windows 7, Python 2.7.*, GDAL 1.11.1, QGIS 2.14.2.


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


==Краткое описание действий==
==Краткое описание действий==
Участки вырубленные за период между датами съемок на более позднем снимки будут иметь значительно большие значения яркости, чем на раннем (естественно при условии наличия снежного покрова и не слишком большого временного интервала). Расчет разности яркости и последующая классификация результата по пороговому значению осуществляется с использованием утилиты gdal_calc.py.
Участки, вырубленные за период между датами съемок, на более позднем снимке будут иметь значительно большие значения яркости, чем на раннем (при условии наличия снежного покрова и не слишком большого временного интервала). Расчет разности яркости и последующая классификация результата по пороговому значению осуществляется с использованием утилиты [http://www.gdal.org/gdal_calc.html 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.
Вопросы радиометрической калибровки достаточно полно рассмотрены в статье http://gis-lab.info/qa/landsat-data-correction.html, учитывая особенности исходных данных выбранных в качестве примера, калибровку проводить не обязательно(если не нужна разность в физических единицах), поэтому все команды приведены для расчетов в DN. Радиометрическую калибровку можно также осуществить с помощью gdal_calc.py, руководствуясь указанной выше статьей.
 
==Создание общей области покрытия==
==Создание общей области покрытия==
Для корректной работы утилиты gdal_calc.py с несколькими растрами необходимо, чтобы входящие растры имели одинаковый экстент. Получить сведения об охвате можно используя утилиту gdalinfo.exe:
Для корректной работы утилиты '''gdal_calc.py''' с несколькими растрами необходимо, чтобы входящие растры имели одинаковый охват (экстент). Получить сведения об охвате можно, используя утилиту '''gdalinfo''':


<pre>gdalinfo rastr</pre>
<pre>gdalinfo rastr</pre>
Строка 24: Строка 25:
*rastr - полный путь к набору растровых данных
*rastr - полный путь к набору растровых данных


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


Необходимо определить максимальные значения для левой и нижней, и минимальные для правой и верхней границ используемых изображений.
Результат работы утилиты будет представлен в виде вывода командной строки:
Обрезка растров производится с помощью утилиты gdalwarp.exe:
[[Файл:Ginfextent.jpg|center|Вывод значений координат углов изображения (охват) в командной строке]]
 
 
Необходимо выбрать максимальные значения для левой и нижней границ и минимальные значения для правой и верхней границ среди охватов всех используемых изображений.
Обрезка растров производится с помощью утилиты '''gdalwarp''':


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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


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

Текущая версия от 10:03, 17 ноября 2016

Эта страница опубликована в основном списке статей сайта
по адресу http://gis-lab.info/qa/gdal_spectral_characteristics.html


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

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

Рис.1 Слева фрагмент сцены LC81810202015048LGN00(февраль 2015 г., см. описание исходных данных), справа тот же участок на сцене LC81810202016083LGN00(март 2016 г.), красный контур-полигональный слой участков леса вырубленных за период между сценами, полученный с помощью автоматизированного алгоритма, рассмотренного в данной статье.

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

ОС Windows 7, 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
  • для значений 0 будет установлен флаг NODATA

Замечание: Для запуска утилит GDAL может понадобится указание полного пути к интерпретатору Python, например:

C:\Python27\python.exe "C:\Program Files\GDAL\gdal_calc.py" (пути с пробелами берутся в кавычки)

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

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

Во избежание пересчета пикселей при трансформации выходной растр необходимо замкнуть на сетку Landsat, используя параметры gdalwarp -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).

Рис.2 Ошибочные результаты классификации (полигональный слой с красным контуром), вызванные геометрическими погрешностями снимков (слева снимок 2015 г., справа 2016 г.)

Уменьшить количество "нежелательных" пикселей, а также заполнить "пустоты" в полигонах можно с помощью фильтра-утилиты 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 "Упростить геометрию".

Замечания

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

См. также