从栅格中提取形状内的点答案

【问题标题】：Extract points within a shape from a raster从栅格中提取形状内的点
【发布时间】：2011-02-15 18:07:03
【问题描述】：

我有一个接近一百万点的光栅文件（基本上是二维数组）。我正在尝试从光栅中提取一个圆（以及圆内的所有点）。为此，使用 ArcGIS 非常缓慢。任何人都可以推荐任何既容易学习又功能强大且足够快的图像处理库来完成这样的事情吗？

谢谢！

【问题讨论】：

【解决方案1】：

有效地提取点子集取决于您使用的确切格式。假设您将栅格存储为一个 numpy 整数数组，您可以像这样提取点：

from numpy import *

def points_in_circle(circle, arr):
    "A generator to return all points whose indices are within given circle."
    i0,j0,r = circle
    def intceil(x):
        return int(ceil(x))
    for i in xrange(intceil(i0-r),intceil(i0+r)):
        ri = sqrt(r**2-(i-i0)**2)
        for j in xrange(intceil(j0-ri),intceil(j0+ri)):
            yield arr[i][j]

points_in_circle 将创建一个返回所有点的生成器。请注意，我使用了yield 而不是return。此函数实际上并不返回点值，而是描述如何找到所有点值。它在圆内的点值上创建一个顺序迭代器。有关yield 工作原理的更多详细信息，请参阅Python documentation。

我使用了这样一个事实，即对于 circle，我们只能在内部点上显式循环。对于更复杂的形状，您可以遍历形状范围的点，然后检查一个点是否属于它。诀窍不是检查每个点，只检查其中的一小部分。

现在举例说明如何使用points_in_circle：

# raster dimensions, 10 million points
N, M = 3200, 3200
# circle center and its radius in index space
i0, j0, r = 70, 20, 12.3

raster = fromfunction(lambda i,j: 100+10*i+j, (N, M), dtype=int)
print "raster is ready"
print raster

pts_iterator = points_in_circle((i0,j0,r), raster) # very quick, do not extract points yet
pts = array(list(pts_iterator)) # actually extract all points
print pts.size, "points extracted, sum = ", sum(pts)

在 1000 万个整数的栅格上，它非常快。

如果您需要更具体的答案，请描述文件格式或在某处放置示例。

【讨论】：

您可能希望删除对in_circle 的调用，并将intcell(j0+r) 替换为intcell(j0+ri)。另请注意，出于对称原因，您对待圆的左侧和底部与右侧和顶部不同：左侧和底部的界限应该是 int(floor(...))。
好的，谢谢。我修复了+ri 并删除了in_circle。左侧和底部（顶部）应该是int(ceil(...))，因为我们只需要内部点。对于右侧和顶部（底部），我使用相同的int(ceil(...))，因为xrange 不包含第二个参数。因此，为了进行对称选择，我在两端使用ceil。
谢谢！速度非常惊人！我不确定为什么这对于提取重叠圆圈不起作用，正如您在下面的评论中提到的那样。
再次感谢！我能够将它与将 arcGIS 栅格转换为 numpy 数组的代码结合起来。 arcGIS 2 分钟完成的工作现在在一秒钟内完成:)
为避免边缘效应（例如负数或超出范围的索引），您可以在 points_in_circle 中的 yield arr[i][j] 之前的一行添加以下条件 if (i >= 0 and i < rows) and (j>=0 and j < cols):。

【解决方案2】：

Numpy 允许你这样做，而且速度非常快：

import numpy

all_points = numpy.random.random((1000, 1000))  # Input array
# Size of your array of points all_points:
(image_size_x, image_size_y) = all_points.shape
# Disk definition:
(center_x, center_y) = (500, 500)
radius = 10

x_grid, y_grid = numpy.meshgrid(numpy.arange(image_size_x),
                                numpy.arange(image_size_y))
# Array of booleans with the disk shape
disk = ((x_grid-center_x)**2 + (y_grid-center_y)**2) <= radius**2

# You can now do all sorts of things with the mask "disk":

# For instance, the following array has only 317 points (about pi*radius**2):
points_in_disk = all_points[disk]
# You can also use masked arrays:
points_in_circle2 = numpy.ma.masked_array(all_points, ~disk)
from matplotlib import pyplot
pyplot.imshow(points_in_circle2)

【讨论】：

从巨大的栅格中提取小形状时，这种方法可能会占用更多内存。事实上，在相同的输入（10M 点，小圆圈）上，它几乎比我的脚本慢三倍，而在 100M 点（和 2GB RAM）下，它几乎慢了 10 倍。但与我的解决方案不同的是，它保留了形状（使用掩码数组时），并且在两次提取相同形状时可能会更好。（顺便说一句，你忘了定义all_points）。
@jetxee：是的，这种方法会消耗内存，并且对于小半径来说速度较慢。但是，代码非常清晰，几乎可以立即在现代机器上执行。

【解决方案3】：

您需要一个可以读取栅格的库。我不知道如何在 Python 中做到这一点，但如果你想用 Java 编程，你可以查看 geotools（尤其是一些新的栅格库集成）。如果你擅长 C，我会推荐使用 GDAL 之类的东西。

如果你想看一个桌面工具，你可以看看用 python 扩展 QGIS 来做上面的操作。

如果我没记错的话，PostGIS 的栅格扩展可能支持基于矢量裁剪栅格。这意味着您需要为数据库中的要素创建圆圈，然后导入您的栅格，但随后您也许可以使用 SQL 来提取您的值。

如果你真的只是一个带有网格数字的文本文件，那么我会接受上面的建议。

【讨论】：

感谢您的回答！我打算使用GDAL从光栅转换为numpy数组