从栅格中提取形状的点

Question

我有一个栅格文件（基本上是2D数组），近一百万分。我正在尝试从栅格中提取一个圆（以及圆圈内的所有点）。使用ArcGI极为慢。谁能提出任何既易于学习的图像处理库，又足够快地来实现这样的东西？

谢谢！

Answer 1

有效提取点子集取决于您正在使用的确切格式。假设您将光栅存储为一个数字整数阵列，则可以提取这样的点：

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. 。此函数实际上并未返回点值，而是描述了如何找到所有功能。它在圆圈内的点值中创建了一个顺序迭代器。看 Python文档 有关更多详细信息 yield 作品。

我使用的事实是，对于圆圈，我们只能在内部点上明确循环。对于更复杂的形状，您可以在形状的范围上循环，然后检查一个点是否属于它。诀窍不是检查 每一个 点，只有一个狭窄的子集。

现在是如何使用的示例 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万个整数的栅格中，这很快。

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

Answer 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)

Answer 3

您需要一个可以阅读栅格的库。我不确定如何在Python中做到这一点，但是如果您想在Java进行编程，您可以查看Geotools（尤其是一些新的栅格库集成）。如果您对CI良好，则可以使用GDAL之类的东西进行推荐。

如果您想查看桌面工具，则可以考虑使用Python扩展QGIS以进行上面的操作。

如果我没记错的话，PostGIS的栅格扩展可能会根据向量支持剪切栅格。这意味着您需要在DB中创建圆圈，然后导入光栅，但是您可以使用SQL来提取值。

如果您实际上只是一个文本文件，其中包含网格中的数字，那么我会处理上面的建议。