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

Ответить
Аватара пользователя
Denis Rykov
Гуру
Сообщения: 3376
Зарегистрирован: 11 апр 2008, 21:09
Статьи: 33
Проекты: 9
Репутация: 526
Ваше звание: Author
Контактная информация:

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

Сообщение Denis Rykov » 15 авг 2009, 02:40

Задача следующая. Имеется растр в формате 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
polygonize1.gif (8.11 КБ) 4522 просмотра
3) Открываем таблицу атрибутов и выбираем только те полигоны, для которых поле "fire_mask" равно 8 или 9:

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

fire_mask = 8 OR fire_mask = 9
4) Сохраняем выбранные полигоны как шейп-файл ("Сохранить выделение как шейп-файл").
5) Открываем получившийся шейп-файл:
polygonize2.gif
polygonize2.gif (1.93 КБ) 4521 просмотр
Таким образом задача решена.

Способ 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
polygonize3.gif (1.83 КБ) 4529 просмотров
Есть какие-нибудь предположения, откуда взялись эти нулевые полигоны. Исходный растр во вложении.
Вложения
A2008137_Days_01.7z
(99.63 КБ) 313 скачиваний
Spatial is now, more than ever, just another column- The Geometry Column.

Аватара пользователя
Максим Дубинин
MindingMyOwnBusiness
Сообщения: 9037
Зарегистрирован: 06 окт 2003, 20:20
Статьи: 231
Проекты: 12/6
Репутация: 713
Ваше звание: NextGIS
Откуда: Москва
Контактная информация:

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

Сообщение Максим Дубинин » 16 авг 2009, 07:38

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

Аватара пользователя
Denis Rykov
Гуру
Сообщения: 3376
Зарегистрирован: 11 апр 2008, 21:09
Статьи: 33
Проекты: 9
Репутация: 526
Ваше звание: Author
Контактная информация:

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

Сообщение Denis Rykov » 16 авг 2009, 08:43

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

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

ogr2ogr -f "ESRI Shapefile" -sql "SELECT * FROM Days01 WHERE ((fire_mask=8) or (fire_mask=9))" Days01_out.shp Days01.shp
Spatial is now, more than ever, just another column- The Geometry Column.

Ответить

Вернуться в «GDAL/OGR»

Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей и 1 гость