Извлечение точек в фигуре из растра

Question

У меня есть растровый файл (в основном, 2D-массив) с почти миллионом точек. Я пытаюсь извлечь круг из растра (и все точки, которые лежат в круге). Использование ArcGIS чрезвычайно медленно для этого. Кто-нибудь может предложить какую-нибудь библиотеку обработки изображений, которая была бы простой в освоении, мощной и достаточно быстрой для чего-то подобного?

Спасибо!

Answer 1

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 2

Эффективное извлечение поднабора точек зависит от того, какой именно формат вы используете. Предполагая, что вы храните свой растр как массив целых чисел, вы можете извлечь точки следующим образом:

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)

На растре из 10 миллионов целых чисел это довольно быстро.

Пожалуйста, опишите формат файла или поместите образец где-нибудь, если вам нужен более конкретный ответ.

Answer 3

Вам нужна библиотека, которая может читать ваш растр. Я не уверен, как это сделать в Python, но вы могли бы взглянуть на геоинструменты (особенно с некоторыми из новой интеграции растровой библиотеки), если вы хотите программировать на Java. Если вы хорошо разбираетесь в C, я бы рекомендовал использовать что-то вроде GDAL.

Если вы хотите взглянуть на настольный инструмент, вы можете посмотреть на расширение QGIS с помощью python для выполнения описанной выше операции.

Если я правильно помню, расширение Raster для PostGIS может поддерживать растровые отсечения на основе векторов. Это означает, что вам нужно будет создать круги для объектов в БД, а затем импортировать растр, но тогда вы сможете использовать SQL для извлечения ваших значений.

Если вы на самом деле просто текстовый файл с числами в сетке, то я бы пошел с предложениями выше.