Векторизация растра

[Denis Rykov]

Задача следующая. Имеется растр в формате GeoTIFF, из которого нужно векторизовать только области, удовлетворяющие определенным значениям. Например, пиксели со значениями {8, 9}.

Способ 1.
1) :

Код: Выделить всё

PATH=C:\OSGeo4W\apps\gdal-dev\bin\;%PATH%
set PYTHONPATH=C:\OSGeo4W\apps\gdal-dev\pymod
gdal_polygonize.py A2008137_Days_01.tif -f "ESRI Shapefile" Days01.shp "Active_fires" "fire_mask"

где A2008137_Days_01.tif - входной растр, Days01.shp - выходной шейп, "Active_fires" и "fire_mask" - необязательные параметры: имя выходного слоя и имя аттрибутивного поля соответственно.
Запускаем скрипт и получаем на выходе файл Days01.shp.

2) Открываем шейп в любой ГИС, например в QGIS:

polygonize1.gif (8.11 КБ) 6794 просмотра

3) Открываем таблицу атрибутов и выбираем только те полигоны, для которых поле "fire_mask" равно 8 или 9:

Код: Выделить всё

fire_mask = 8 OR fire_mask = 9

4) Сохраняем выбранные полигоны как шейп-файл ("Сохранить выделение как шейп-файл").
5) Открываем получившийся шейп-файл:

polygonize2.gif (1.93 КБ) 6793 просмотра

Таким образом задача решена.

Способ 2.
Воспользуемся скриптом classify.py.
1) Отредактируем его в соответствии с задачей:

Код: Выделить всё

import gdal
import gdalnumeric
try:
    import numpy
except:
    import Numeric as numpy


class_defs = [(8, 8,8), (9, 9, 9)]

src_ds = gdal.Open('A2008137_Days_01.tif')
xsize = src_ds.RasterXSize
ysize = src_ds.RasterYSize

src_image = gdalnumeric.LoadFile( 'A2008137_Days_01.tif' )

dst_image = numpy.zeros((ysize,xsize))

for class_info in class_defs:
    class_id = class_info[0]
    class_start = class_info[1]
    class_end = class_info[2]

    class_value = numpy.ones((ysize,xsize)) * class_id
    
    mask = numpy.bitwise_and(
        numpy.greater_equal(src_image,class_start),
        numpy.less_equal(src_image,class_end))

    dst_image = numpy.choose( mask, (dst_image,class_value) )

gdalnumeric.SaveArray( dst_image, 'classes.tif' ) 

2) Запускаем скрипт и получаем на выходе файл classes.tif с нужными пикселями. Полученный файл утратил исходную информацию о географической привязки, поэтому восстановим ее с помощью tfw-файла:

Код: Выделить всё

gdal_translate -of GTiff -co "TFW=YES" A2008137_Days_01.tif A2008137_Days_01_tfw.tif

3) Удаляем файл A2008137_Days_01_tfw.tif, а файл A2008137_Days_01_tfw.tfw переименовываем в classes.tfw.
4) Снова воспользуемся утилитой gdal_polygonize:

Код: Выделить всё

PATH=C:\OSGeo4W\apps\gdal-dev\bin\;%PATH%
set PYTHONPATH=C:\OSGeo4W\apps\gdal-dev\pymod
gdal_polygonize.py classes.tif -f "ESRI Shapefile" classes.shp "Active_fires" "fire_mask"

5) Получили требуемый шейп-файл.

На первый взгляд все хорошо, 2-й способ позволяет автоматизировать процесс векторизации растров, но есть небольшая проблема. Если в нашем случае посмотреть таблицу атрибутов шейп файла созданного первым способом и вторым - то мы увидим, что они отличаются. В таблице, полученной первым способом - 645 записей, а вторым - 649, причем эти 4 отличные записи имеют аттрибут 0. Да и сам шейп выглядит несколько странно:

polygonize3.gif (1.83 КБ) 6801 просмотр

Есть какие-нибудь предположения, откуда взялись эти нулевые полигоны. Исходный растр во вложении.

[Максим Дубинин]

В таблице, полученной первым способом - 645 записей, а вторым - 649, причем эти 4 отличные записи имеют аттрибут 0.

это острова-дырки + внешний полигон, можно использовать ogr2ogr -sql "SELECT .... чтобы выбрать то что нужно.

[Denis Rykov]

Спасибо. Таким образом, лучшим способом для выполнения задания в пакетном режиме является Способ 1, только все манипуляции с QGIS заменяются одной командой:

Код: Выделить всё

ogr2ogr -f "ESRI Shapefile" -sql "SELECT * FROM Days01 WHERE ((fire_mask=8) or (fire_mask=9))" Days01_out.shp Days01.shp