4

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

谢谢！

4

3 回答 3

2

有效地提取点子集取决于您使用的确切格式。假设您将栅格存储为一个 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. 此函数实际上并不返回点值，而是描述如何找到所有点值。它在圆内的点值上创建一个顺序迭代器。有关如何工作的更多详细信息，请参阅Python 文档。yield

我使用了这样一个事实，即对于 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 万个整数的栅格上，它非常快。

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

于 2010-05-05T14:49:38.630 回答
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)
于 2010-05-05T15:43:51.763 回答
1

您需要一个可以读取栅格的库。我不确定如何在 Python 中做到这一点，但如果你想用 Java 编程，你可以查看 geotools（尤其是一些新的栅格库集成）。如果你对 CI 很好，建议使用 GDAL 之类的东西。

如果您想查看桌面工具，您可以查看使用 python 扩展 QGIS 以执行上述操作。

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

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

于 2010-05-07T00:38:47.663 回答