从坐标参考系中的多边形边界获取光栅中的像素坐标

Question

我有一些相对较大（10980 x 10980像素）的GeoTiff文件，它们都对应于相同的地理区域（并且具有相同的坐标参考系），我有大量对应于地块的多边形（100000+），我想从每个图像文件中提取对应于每个多边形的像素。目前，我这样做的方式是使用shapely多边形和rasterio.mask.mask函数，如下所示：

for filename in image_files:
    with rasterio.open(filename) as src:
        for shape in shapes:
            data, _ = rasterio.mask.mask(src, [shape], crop=True)

这在经验上相当缓慢。如果我预先计算了遮罩索引，那么我只需要读取每个图像的整个数据一次，然后使用预先计算的索引来提取每个多边形的相关像素（我不需要它们处于正确的二维配置中，我只需要值），这非常快。但我不知道是否有一种快速的方法来获取这些像素索引。我知道我可以使用rasterio的raster_geometry_mask函数来获得整个图像大小的遮罩，然后使用numpy来获得多边形内部元素的索引，但是这样就不需要为每个多边形构造10980x10980数组来制作遮罩，而且速度非常慢

Answer 1

我最后做的是，当我打开第一个图像，然后对每个多边形

• 使用图像变换将多边形转换为像素坐标，并在整数像素坐标中查找包含多边形的矩形边界框
• 若要确定边界框中的哪些像素实际位于多边形中，请为每个像素构造形状良好的多边形，并使用.intersects()方法（如果只希望包含完全位于多边形内部的像素，可以使用.contains()）。（我不确定这是否会很慢，但事实证明并非如此。）
• 保存每个多边形中所有像素的坐标对列表

然后，对于您打开的每个新图像，您只需读取整个图像数据并为每个多边形的部分建立索引，因为您已经有了像素索引

代码大致如下所示：

import math
import numpy
import pyproj
import rasterio.mask
from shapely.geometry import Polygon

shape_pixels = None
for filename in image_files:
    with rasterio.open(filename) as src:
        if shape_pixels is None:
            projector = pyproj.Proj(src.crs)
            pixelcoord_shapes = [
                Polygon(zip(*(~src.transform * numpy.array(projector(*zip(*shape.boundary.coords))))))
                for shape in shapes
            ]

            pixels_per_shape = []
            for shape in shapes:
                xmin = max(0, math.floor(shape.bounds[0]))
                ymin = max(0, math.floor(shape.bounds[1]))
                xmax = math.ceil(shape.bounds[2])
                ymax = math.ceil(shape.bounds[3])
                pixel_squares = {}
                for j in range(xmin, xmax+1):
                    for i in range(ymin, ymax+1):
                        pixel_squares[(i, j)] = Polygon.from_bounds(j, i, j+1, i+1)
                pixels_per_shape.append([
                    coords for (coords, pixel) in pixel_squares.items()
                    if shape.intersects(pixel)
                ])

        whole_data = src.read()

        for pixels in pixels_per_shape:
            ivals, jvals = zip(*pixels)
            shape_data = whole_data[0, ivals, jvals]
            ...