Как полигонизировать растр в стройные полигоны

14

Я ищу Python-решение с открытым исходным кодом для преобразования растра в полигон (без ArcPy).

Я знал функцию GDAL для преобразования растра в многоугольник, и вот руководство: http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-raster-band

Тем не менее, я ожидаю, что на выходе могут быть фигурные полигоны или любой объект, временно находящийся в памяти, а не сохраненный в виде файла. Есть ли пакет или код для решения этой проблемы?

Если растр был обработан в массиве NumPy, подход указан ниже.

Вики Ляу
источник

Ответы:

21

Используйте растерио Шона Джиллиса . Его можно легко комбинировать с Fiona (чтение и запись шейп-файлов) и shapely того же автора.

В скрипте rasterio_polygonize.py начало

import rasterio
from rasterio.features import shapes
mask = None
with rasterio.drivers():
    with rasterio.open('a_raster') as src:
        image = src.read(1) # first band
        results = (
        {'properties': {'raster_val': v}, 'geometry': s}
        for i, (s, v) 
        in enumerate(
            shapes(image, mask=mask, transform=src.affine)))

Результатом является генератор функций GeoJSON

 geoms = list(results)
 # first feature
 print geoms[0]
 {'geometry': {'type': 'Polygon', 'coordinates': [[(202086.577, 90534.3504440678), (202086.577, 90498.96207), (202121.96537406777, 90498.96207), (202121.96537406777, 90534.3504440678), (202086.577, 90534.3504440678)]]}, 'properties': {'raster_val': 170.52000427246094}}

Что вы можете превратить в стройную геометрию

from shapely.geometry import shape
print shape(geoms[0]['geometry'])
POLYGON ((202086.577 90534.35044406779, 202086.577 90498.96206999999, 202121.9653740678 90498.96206999999, 202121.9653740678 90534.35044406779, 202086.577 90534.35044406779))

Создайте геопанду Dataframe и включите простые в использовании функции пространственного соединения, построения графиков, сохранения в виде геоджона, шейп-файла ESRI и т. Д.

geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)
ген
источник
Если растр был обработан как массив NumPy, есть ли способ преобразовать массив NumPy в полигоны? Благодарность!
Вики Ляу
В теории да
ген
1
Переменная и параметр маски в вашем примере кажутся ненужными. Однако я бы порекомендовал добавить if value > src.nodataк списку понимание, чтобы использовать значение нодаты источника и отбросить любые формы, которые ему соответствуют. Не уверен, что произойдет, если не будет nodata vaue. : o)
bugmenot123
3
тем временем они изменили rasterio.drivers на rasterio.Env и src.affine на src.transform
Лев
4

Вот моя реализация.

from osgeo import ogr, gdal, osr
from osgeo.gdalnumeric import *  
from osgeo.gdalconst import * 
import fiona
from shapely.geometry import shape
import rasterio.features

#segimg=glob.glob('Poly.tif')[0]
#src_ds = gdal.Open(segimg, GA_ReadOnly )
#srcband=src_ds.GetRasterBand(1)
#myarray=srcband.ReadAsArray() 
#these lines use gdal to import an image. 'myarray' can be any numpy array

mypoly=[]
for vec in rasterio.features.shapes(myarray):
    mypoly.append(shape(vec))

Способ установки rasterio можно получить с помощью 'conda install -c https://conda.anaconda.org/ioos rasterio', если возникает проблема с установкой.

Вики Ляу
источник
Результатом растерио является непосредственно массив numpy, поэтому вам не нужноmyarray=srcband.ReadAsArray() #or any array
ген
@Gene Я пересмотрел записку. Эта строка (myarray = srcband.ReadAsArray ()) использует gdal для импорта изображения.
Вики Ляу
импортировать изображение в виде массива и растерио напрямую импортировать изображение в виде массива
ген
Это сработало для меня, хотя я должен был индексировать vec, потому что он возвращался как кортеж. shape (vec [0])
user2723146