17

私は何百万もの地理的ポイントを持っています。これらのそれぞれについて、すべての「隣接ポイント」、つまり、数百メートルなどの半径内の他のすべてのポイントを見つけたいと考えています。

この問題に対する単純な O(N^2) ソリューションがあります。単純に、すべての点のペアの距離を計算します。ただし、適切な距離メトリック (地理的距離) を扱っているため、これを行うためのより迅速な方法があるはずです。

これをPython内で実行したいと思います。頭に浮かぶ 1 つの解決策は、何らかのデータベース (GIS 拡張機能を備えた mySQL、PostGIS) を使用し、そのようなデータベースが何らかのインデックスを使用して上記の操作を効率的に実行することを期待することです。ただし、そのようなテクノロジを構築して学習する必要がない、より単純なものを好みます。

いくつかのポイント

  • 「隣人を探す」操作を何百万回も実行します
  • データは静的なままです
  • ある意味単純な問題なので、それを解くpythonコードを出してほしいです。

Pythonコードの観点から言えば、次のようなものが必要です:

points = [(lat1, long1), (lat2, long2) ... ] # this list contains millions lat/long tuples
points_index = magical_indexer(points)
neighbors = []
for point in points:
    point_neighbors = points_index.get_points_within(point, 200) # get all points within 200 meters of point
    neighbors.append(point_neighbors)
4

2 に答える 2

7

scipy

まず最初に： kdツリーなどのようなことを行うための既存のアルゴリズムがあります。Scipyには、指定された範囲内のすべてのポイントを見つけることができるPython実装cKDtreeがあります。

二分探索

ただし、実行していることによっては、そのようなものを実装することは簡単ではない場合があります。さらに、ツリーの作成はかなり複雑であり（潜在的にかなりのオーバーヘッド）、以前に使用した単純なハックで回避できる可能性があります。

  1. データセットのPCAを計算します。最も重要な方向が最初になり、直交する（大きくない）2番目の方向が2番目になるように、データセットを回転させます。これをスキップしてXまたはYを選択することもできますが、計算コストが低く、通常は簡単に実装できます。XまたはYを選択するだけの場合は、分散が大きい方向を選択します。
  2. ポイントを主方向で並べ替えます（この方向をXと呼びます）。
  3. 特定のポイントの最近傍を見つけるには、バイナリ検索によってXに最も近いポイントのインデックスを見つけます（ポイントがすでにコレクションにある場合は、このインデックスをすでに知っているので、検索する必要はありません）。これまでの最良の一致と検索ポイントからの距離を維持しながら、次と前のポイントを繰り返し調べます。Xの差がこれまでのベストマッチまでの距離以上になると、見るのをやめることができます（実際には、通常はごくわずかなポイントです）。
  4. 指定された範囲内のすべてのポイントを見つけるには、Xの差が範囲を超えるまで停止しないことを除いて、手順3と同じように実行します。

事実上、O（N log（N））の前処理を行っており、ポイントの分布が悪い場合は、各ポイントについておおよそo（sqrt（N））-またはそれ以上です。ポイントがほぼ均一に分布している場合、Xで最も近い隣接ポイントよりも近いポイントの数は、Nの平方根のオーダーになります。多くのポイントが範囲内にある場合は効率が低下しますが、ブルートフォースよりもはるかに悪くなることはありません。

この方法の利点の1つは、すべてが非常に少ないメモリ割り当てで実行可能であり、ほとんどの場合、非常に優れたメモリローカリティで実行できることです。つまり、明らかな制限にもかかわらず、非常に優れたパフォーマンスを発揮します。

ドロネー三角形分割

別のアイデア：Delauney三角形分割が機能する可能性があります。Delauney三角形分割の場合、任意のポイントの最近傍が隣接ノードであると見なされます。直感的には、検索中に、クエリポイントからの絶対距離に基づいてヒープ（優先度付きキュー）を維持できます。最も近いポイントを選択し、範囲内にあることを確認し、範囲内にある場合は、すべての隣接ポイントを追加します。このような点を見逃すことは不可能だと思いますが、確実にもっと注意深く見る必要があります...

于 2011-06-16T11:54:41.460 に答える
7

Eamon からヒントを得て、SciPy に実装された btree を使用した簡単な解決策を思いつきました。

from scipy.spatial import cKDTree
from scipy import inf

max_distance = 0.0001 # Assuming lats and longs are in decimal degrees, this corresponds to 11.1 meters
points = [(lat1, long1), (lat2, long2) ... ]
tree = cKDTree(points)

point_neighbors_list = [] # Put the neighbors of each point here

for point in points:
    distances, indices = tree.query(point, len(points), p=2, distance_upper_bound=max_distance)
    point_neighbors = []
    for index, distance in zip(indices, distances):
        if distance == inf:
            break
        point_neighbors.append(points[index])
    point_neighbors_list.append(point_neighbors)
于 2011-06-16T13:15:14.833 に答える