查找与多边形相交的栅格像元的索引

Question

我想获取落入多边形特征或与多边形特征相交的所有栅格像元的索引列表（行、列）。在 python 中寻找解决方案，最好使用 gdal/ogr 个模块。

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

Answer 1

如果您想手动执行此操作，您需要测试每个单元格：方形 v 多边形相交和方 v 线交点。

如果您将每个正方形视为一个 2d 点，这会变得更容易 - 它现在是点 v 多边形问题。在游戏开发论坛中查看碰撞算法。

祝你好运！

Answer 2

显然，Rutger 的上述解决方案是解决此问题的方法，但我将保留我的解决方案。我开发了一个脚本来完成我需要的内容：

获取我要检查的每个矢量特征的边界框
使用边界框限制计算window（确定栅格的哪一部分可能有交点）
遍历这部分栅格内的像元并为每个像元构建多边形几何图形
使用ogr.Geometry.Intersects()检查单元格是否与多边形要素相交

请注意，我只定义了方法，但我认为实现应该非常清楚——只需使用适当的参数（ogr.Geometry 对象和地理变换矩阵）调用 match_cells。代码如下：

from osgeo import ogr

# Convert projected coordinates to raster cell indices
def parse_coords(x,y,gt):
    row,col = None,None
    if x:
        col = int((x - gt[0]) // gt[1])
        # If only x coordinate is provided, return column index
        if not y:
            return col
    if y:
        row = int((y - gt[3]) // gt[5])
        # If only x coordinate is provided, return column index
        if not x:
            return row
    return (row,col)

# Construct polygon geometry from raster cell
def build_cell((row,col),gt):
    xres,yres = gt[1],gt[5]
    x_0,y_0 = gt[0],gt[3]
    top = (yres*row) + y_0
    bottom = (yres*(row+1)) + y_0
    right = (xres*col) + x_0
    left = (xres*(col+1)) + x_0
    # Create ring topology
    ring = ogr.Geometry(ogr.wkbLinearRing)
    ring.AddPoint(left,bottom)
    ring.AddPoint(right,bottom)
    ring.AddPoint(right,top)
    ring.AddPoint(left,top)
    ring.AddPoint(left,bottom)
    # Create polygon
    box = ogr.Geometry(ogr.wkbPolygon)
    box.AddGeometry(ring)
    return box

# Iterate over feature geometries & check for intersection
def match_cells(inputGeometry,gt):
    matched_cells = []
    for f,feature in enumerate(inputGeometry):
        geom = feature.GetGeometryRef()
        bbox = geom.GetEnvelope()
        xmin,xmax = [parse_coords(x,None,gt) for x in bbox[:2]]
        ymin,ymax = [parse_coords(None,y,gt) for y in bbox[2:]]
        for cell_row in range(ymax,ymin+1):
            for cell_col in range(xmin,xmax+1):
                cell_box = build_cell((cell_row,cell_col),gt)
                if cell_box.Intersects(geom):
                    matched_cells += [[(cell_row,cell_col)]]
    return matched_cells

Answer 3

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

多边形示例

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)

查找与多边形相交的栅格像元的索引

Find indices of raster cells that intersect with a polygon

python

gdal

ogr

多边形示例

栅格化

转换为指数