从栅格中提取形状内的点

Question

我有一个光栅文件（基本上就是一个二维数组），里面有将近一百万个点。我想从这个光栅中提取一个圆形区域，以及所有在这个圆形区域内的点。用ArcGIS来做这个实在是太慢了。有没有人能推荐一些既简单易学又强大、速度够快的图像处理库，适合做这样的事情？

谢谢！

Answer 1

你需要一个可以读取栅格数据的库。我不太确定在Python中怎么做，但如果你想用Java编程，可以看看geotools（特别是它的一些新栅格库集成）。如果你对C语言比较熟悉，我建议使用GDAL这样的工具。

如果你想找一个桌面工具，可以考虑用Python扩展QGIS来完成上面的操作。

如果我没记错的话，PostGIS的栅格扩展可能支持根据矢量数据来裁剪栅格。这意味着你需要在数据库中创建你的圆形特征，然后导入你的栅格数据，之后你就可以用SQL提取你的数值了。

如果你只是一个包含数字的文本文件，那我建议你参考上面的建议。

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

高效提取一部分点的方式，取决于你使用的具体格式。假设你把栅格数据存储为一个整数的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文档。

我利用了一个事实：对于圆形，我们可以明确地只遍历内部的点。对于更复杂的形状，你可以遍历形状的边界点，然后检查某个点是否属于这个形状。关键是不要检查每一个点，而是只检查一小部分。

下面是如何使用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万个整数的栅格上，这个过程相当快。

如果你需要更具体的答案，请描述一下文件格式或者提供一个样本。