4

百万点近くのラスターファイル（基本的には2D配列）があります。ラスター（および円内にあるすべてのポイント）から円を抽出しようとしています。このため、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。

円の場合、内側の点のみを明示的にループできるという事実を使用しました。より複雑な形状の場合は、形状の範囲のポイントをループして、ポイントがその形状に属しているかどうかを確認できます。秘訣は、すべてのポイントをチェックするのではなく、それらの狭いサブセットのみをチェックすることです。

使用方法の例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のラスター拡張機能はベクトルに基づくクリッピングラスターをサポートしている可能性があります。つまり、DB内のフィーチャに円を作成してからラスターをインポートする必要がありますが、SQLを使用して値を抽出できる場合があります。

あなたが本当にグリッド内の数字を含む単なるテキストファイルであるなら、私は上記の提案に行きます。

于 2010-05-07T00:38:47.663 に答える