python - 查找与多边形相交的栅格单元的索引

Question

我想获取位于多边形要素内或与多边形要素相交的所有栅格单元的索引(行，列)列表。在 python 中寻找解决方案，最好使用 gdal/ogr 模块。

其他帖子建议对多边形进行栅格化，但如果可能的话，我宁愿直接访问单元索引。

Answer 1

由于您没有提供有效的示例，因此不清楚您的出发点是什么。我制作了一个包含 1 个多边形的数据集，如果您有一个包含多个多边形的数据集，但只想定位特定的多边形，则可以将 SQLStatement 或 where 添加到 gdal。光栅化调用。

示例多边形

geojson = """{"type":"FeatureCollection",
"name":"test",
"crs":{"type":"name","properties":{"name":"urn:ogc:def:crs:OGC:1.3:CRS84"}},
"features":[
{"type":"Feature","properties":{},"geometry":{"type":"MultiPolygon","coordinates":[[[[-110.254,44.915],[-114.176,37.644],[-105.729,36.41],[-105.05,43.318],[-110.254,44.915]]]]}}
]}"""

光栅化

可以使用gdal.Rasterize完成光栅化。您需要指定目标网格的属性。如果没有预定义的网格，则可以从多边形本身提取这些网格

ds = gdal.Rasterize('/vsimem/tmpfile', geojson, xRes=1, yRes=-1, allTouched=True,
                    outputBounds=[-120, 30, -100, 50], burnValues=1, 
                    outputType=gdal.GDT_Byte)
mask = ds.ReadAsArray()
ds = None
gdal.Unlink('/vsimem/tmpfile')

转换为索引

可以使用 Numpy 从栅格化多边形中检索索引:

y_ind, x_ind = np.where(mask==1)