我正在寻找将栅格转换为多边形(无ArcPy)的开源python解决方案。

我确实知道将栅格转换为多边形的GDAL函数,这是手册:http://pcjericks.github.io/py-gdalogr-cookbook/raster_layers.html#polygonize-a-光栅带

尽管如此,我希望输出可以是形状多边形或任何对象,暂时存储在内存中,而不是另存为文件。是否有任何程序包或代码可解决此问题?

如果栅格是在numpy数组中处理的,则下面列出了该方法。

#1 楼

使用Sean Gillies的栅格。它可以轻松地与Fiona(读写shapefile)结合使用,并且形状相同。
在脚本rasterio_polygonize.py
中,开始是
import rasterio
from rasterio.features import shapes
mask = None
with rasterio.Env():
    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.transform)))

结果是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))

创建geopandas数据框并启用易于使用的空间连接,绘图,另存为geojson,ESRI shapefile的功能等
geoms = list(results)
import geopandas as gp
gpd_polygonized_raster  = gp.GeoDataFrame.from_features(geoms)


评论


如果栅格已作为numpy数组处理,是否可以将numpy数组转换为多边形?谢谢!

– Vicky Liau
16-4-3在18:45



从理论上讲,是的-

–基因
16年4月3日在19:17

您的示例中的mask变量和参数似乎是不必要的。但是,我建议将if> src.nodata添加到列表推导中,以利用源的nodata值并丢弃与之对应的任何形状。不知道如果没有数据的话会发生什么。 :o)

–bugmenot123
17年4月28日在18:58

同时,他们将rasterio.drivers更改为rasterio.Env,并将src.affine更改为src.transform

–狮子座
18/12/3在9:31

当我尝试这个时,我得到一个多边形,其crs为None。所以我尝试将crs设置为原始图像,但随后将其显示在完全错误的位置。有任何想法吗?

– Emtomp
20年4月1日在14:50

#2 楼

这是我的实现。

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',如果存在的话是安装问题。

评论


rasterio的结果直接是一个numpy数组,因此您不需要myarray = srcband.ReadAsArray()#或任何数组

–基因
16年4月4日在6:16

@gene我修改了笔记。此行(myarray = srcband.ReadAsArray())使用gdal导入图像。

– Vicky Liau
16年4月4日在18:21

将图像作为numpy数组导入,而rasterio将图像直接作为numpy数组导入

–基因
16年4月4日在18:51

尽管我必须为vec编制索引,因为它作为元组返回,但它对我有用。形状(vec [0])

– Narmaps
19-10-9的17:00