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

Question

我有一些相对较大的 GeoTiff 文件（10980 x 10980 像素），它们都对应于相同的地理区域（并且具有相同的坐标参考系统），并且我有大量的多边形（100,000+ ) 对应地块，我想从每个图像文件中提取与每个多边形对应的像素。目前，我这样做的方法是使用匀称的多边形和 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)

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

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]
            ...

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

Get pixel coordinates in raster from polygon boundary in coordinate reference system

python

gis

raster

rasterio