python - 500 万行の座標データの処理を高速化

Question

500 万行を超える地理位置情報データを含む 2 つの列 (緯度、経度) を持つ csv ファイルがあります。リスト内の他のポイントから 5 マイル以内にないポイントを特定し、すべてを別の CSV に出力する必要があります。この CSV には、別のポイントが 5 マイル以内にある場合と、そうでない場合に追加の列 ( CloseToAnotherPoint) があります。だ。TrueFalse

これが私の現在の解決策ですgeopy（Web呼び出しを行わず、関数を使用して距離を計算するだけです）：

from geopy.point import Point
from geopy.distance import vincenty
import csv


class CustomGeoPoint(object):
    def __init__(self, latitude, longitude):
        self.location = Point(latitude, longitude)
        self.close_to_another_point = False


try:
    output = open('output.csv','w')
    writer = csv.writer(output, delimiter = ',', quoting=csv.QUOTE_ALL)
    writer.writerow(['Latitude', 'Longitude', 'CloseToAnotherPoint'])

    # 5 miles
    close_limit = 5
    geo_points = []

    with open('geo_input.csv', newline='') as geo_csv:
        reader = csv.reader(geo_csv)
        next(reader, None) # skip the headers
        for row in reader:
            geo_points.append(CustomGeoPoint(row[0], row[1]))

    # for every point, look at every point until one is found within 5 miles
    for geo_point in geo_points:
        for geo_point2 in geo_points:
            dist = vincenty(geo_point.location, geo_point2.location).miles
            if 0 < dist <= close_limit: # (0,close_limit]
                geo_point.close_to_another_point = True
                break
        writer.writerow([geo_point.location.latitude, geo_point.location.longitude,
                         geo_point.close_to_another_point])

finally:
    output.close()

これを見ればわかるかもしれませんが、このソリューションは非常に遅いです。実際、とても遅いので、3 日間実行しましたが、まだ完了しませんでした!

内部ループが他のすべてのポイントを見る必要がないように、データをチャンク (複数の CSV ファイルなど) に分割しようと考えましたが、境界線を確認する方法を理解する必要があります。隣接するセクションの境界に対してチェックされた各セクションの、そしてそれはあまりにも複雑に思えます。

これをより速くする方法についての指針はありますか？

Answer 1

Python calculate lot of distancesへの回答がすぐに指摘するように、これは kD ツリーの古典的な使用例です。

Python を使用して同様の座標を一致させるにはどうすればよいですか?への回答に示されているように、スイープ ライン アルゴリズムを使用することもできます。

これが、質問に合わせたスイープ ライン アルゴリズムです。私のラップトップでは、5M のランダム ポイントを実行するのに 5 分未満かかります。

import itertools as it
import operator as op
import sortedcontainers     # handy library on Pypi
import time

from collections import namedtuple
from math import cos, degrees, pi, radians, sqrt
from random import sample, uniform

Point = namedtuple("Point", "lat long has_close_neighbor")

miles_per_degree = 69

number_of_points = 5000000
data = [Point(uniform( -88.0,  88.0),     # lat
              uniform(-180.0, 180.0),     # long
              True
             )
        for _ in range(number_of_points)
       ]

start = time.time()
# Note: lat is first in Point, so data is sorted by .lat then .long.
data.sort()

print(time.time() - start)

# Parameter that determines the size of a sliding lattitude window
# and therefore how close two points need to be to be to get flagged.
threshold = 5.0  # miles
lat_span = threshold / miles_per_degree
coarse_threshold = (.98 * threshold)**2

# Sliding lattitude window.  Within the window, observations are
# ordered by longitude.
window = sortedcontainers.SortedListWithKey(key=op.attrgetter('long'))

# lag_pt is the 'southernmost' point within the sliding window.
point = iter(data)
lag_pt = next(point)

milepost = len(data)//10

# lead_pt is the 'northernmost' point in the sliding window.
for i, lead_pt in enumerate(data):
    if i == milepost:
        print('.', end=' ')
        milepost += len(data)//10

    # Dec of lead_obs represents the leading edge of window.
    window.add(lead_pt)

    # Remove observations further than the trailing edge of window.
    while lead_pt.lat - lag_pt.lat > lat_span:
        window.discard(lag_pt)
        lag_pt = next(point)

    # Calculate 'east-west' width of window_size at dec of lead_obs
    long_span = lat_span / cos(radians(lead_pt.lat))
    east_long = lead_pt.long + long_span
    west_long = lead_pt.long - long_span

    # Check all observations in the sliding window within
    # long_span of lead_pt.
    for other_pt in window.irange_key(west_long, east_long):

        if other_pt != lead_pt:
            # lead_pt is at the top center of a box 2 * long_span wide by 
            # 1 * long_span tall.  other_pt is is in that box. If desired, 
            # put additional fine-grained 'closeness' tests here. 

            # coarse check if any pts within 80% of threshold distance
            # then don't need to check distance to any more neighbors
            average_lat = (other_pt.lat + lead_pt.lat) / 2
            delta_lat   = other_pt.lat - lead_pt.lat
            delta_long  = (other_pt.long - lead_pt.long)/cos(radians(average_lat))

            if delta_lat**2 + delta_long**2 <= coarse_threshold:
                break

            # put vincenty test here
            #if 0 < vincenty(lead_pt, other_pt).miles <= close_limit:
            #    break

    else:
        data[i] = data[i]._replace(has_close_neighbor=False)

print()      
print(time.time() - start)

Answer 2

あなたが何をしているのか見てみましょう。

1. すべてのポイントを という名前のリストに読み込みますgeo_points。

   では、リストがソートされているかどうか教えていただけますか? それがソートされた場合、私たちは間違いなくそれを知りたいからです。並べ替えは貴重な情報であり、特に 500 万件ものものを扱っている場合はそうです。

2. すべてをループしますgeo_points。あなたによると、それは500万です。

3. 外側のループ内で、500 万個すべてを再度geo_pointsループします。

4. 2 つのループ アイテム間の距離をマイル単位で計算します。

5. 距離がしきい値よりも小さい場合は、最初のポイントでその情報を記録し、内側のループを停止します。

6. 内側のループが停止したら、外側のループ項目に関する情報を CSV ファイルに書き込みます。

いくつかのことに注意してください。まず、外側のループで 500 万回ループしています。そして、内側のループで 500 万回ループしています。

これが O(n²) の意味です。

次に誰かが「これは O(log n) ですが、それ以外は O(n log n) です」と話しているのを見かけたら、この経験を思い出してください。この場合、n が 5,000,000 である n² アルゴリズムを実行していることになります。クソ、ダニット？

とにかく、あなたにはいくつかの問題があります。

問題 1: 最終的には、すべてのポイントをそれ自体と比較することになります。これは、距離がゼロであるべきです。つまり、それらはすべて、距離のしきい値内にあるとマークされます。プログラムが終了すると、すべてのセルが True とマークされます。

問題 2: ポイント #1 をポイント #12345 と比較し、それらが互いにしきい値距離内にある場合、ポイント #1 に関する情報を記録していることになります。しかし、他の点について同じ情報を記録することはありません。ポイント #12345 (geo_point2) が反射的にポイント #1 のしきい値内にあることはわかっていますが、それを書き留めていません。つまり、500 万回以上の比較をスキップするチャンスを逃しています。

問題 3: ポイント #1 とポイント #2 を比較し、それらがしきい値距離内にない場合、ポイント #2 とポイント #1 を比較するとどうなりますか? 内部ループは毎回リストの先頭から開始していますが、リストの先頭とリストの末尾を既に比較していることはわかっています。i in range(0, 5million)外側のループを goにし、内側のループを go にするだけで、問題のスペースを半分に減らすことができますj in range(i+1, 5million)。

答えは？

平面上の緯度と経度を考慮してください。5 マイル以内にポイントがあるかどうかを知りたいとします。ポイント 1 を中心とする 10 マイルの正方形について考えてみましょう。これは、(X1, Y1) を中心とする正方形で、左上隅が (X1 - 5 マイル、Y1 + 5 マイル) で、右下隅が (X1 + 5 マイル、Y1 - 5 マイル) です。ポイントがその正方形内にある場合、ポイント #1 から 5 マイル以内にない可能性があります。しかし、その広場の外にある場合は、5 マイル以上離れていることは間違いありません。

@SeverinPappadeaux が指摘しているように、地球のような回転楕円体での距離は、平面での距離とまったく同じではありません。しかし、だから何？違いを考慮して正方形を少し大きく設定し、先に進みます!

ソート済みリスト

これが、ソートが重要な理由です。すべてのポイントが X、次に Y (または Y、次に X など) でソートされていて、それを知っていれば、本当にスピードアップできます。X (または Y) 座標が大きくなりすぎたときにスキャンを停止するだけで、500 万ポイントを通過する必要がないからです。

それはどのように機能しますか？前と同じ方法ですが、内側のループには次のようなチェックがあります。

five_miles = ... # Whatever math, plus an error allowance!
list_len = len(geo_points) # Don't call this 5 million times

for i, pi in enumerate(geo_points):

    if pi.close_to_another_point:
        continue   # Remember if close to an earlier point

    pi0max = pi[0] + five_miles
    pi1min = pi[1] - five_miles
    pi1max = pi[1] + five_miles

    for j in range(i+1, list_len):
        pj = geo_points[j]
        # Assumes geo_points is sorted on [0] then [1]
        if pj[0] > pi0max:
            # Can't possibly be close enough, nor any later points
            break
        if pj[1] < pi1min or pj[1] > pi1max:
            # Can't be close enough, but a later point might be
            continue

        # Now do "real" comparison using accurate functions.
        if ...:
            pi.close_to_another_point = True
            pj.close_to_another_point = True
            break

私はそこで何をしているのですか？まず、いくつかの数値をローカル変数に入れています。次に、値と外側の点への参照をenumerate与えるために使用しています。（あなたが呼んだもの）。次に、このポイントが別のポイントに近いことをすでに知っているかどうかをすばやく確認しています。igeo_point

そうでない場合は、スキャンする必要があります。したがって、リスト内の「後の」ポイントのみをスキャンしています。これは、外側のループが初期のポイントをスキャンすることを知っているためです。ポイントをそれ自体と比較したくないことは間違いありません。いくつかの一時変数を使用して、外側のループを含む計算の結果をキャッシュしています。内側のループ内で、一時オブジェクトに対していくつかのばかげた比較を行います。2 つの点が互いに近いかどうかはわかりませんが、確実に近くないかどうかを確認してスキップすることはできます。

最後に、単純なチェックに合格した場合は、先に進んで高価なチェックを行います。実際にチェックに合格した場合は、後で 2 番目のポイントをスキップできるように、必ず両方のポイントで結果を記録してください。

ソートされていないリスト

しかし、リストがソートされていない場合はどうなるでしょうか?

@RootTwo は、kD ツリーを示します (ここで、D は「次元」を表し、この場合の k は「2」です)。二分探索木について既に知っている場合、アイデアは非常に単純です。次元を循環し、ツリーの偶数レベルで X を比較し、奇数レベルで Y を比較します (またはその逆)。アイデアは次のようになります。

def insert_node(node, treenode, depth=0):
    dimension = depth % 2  # even/odd -> lat/long
    dn = node.coord[dimension]
    dt = treenode.coord[dimension]

    if dn < dt:
        # go left
        if treenode.left is None:
            treenode.left = node
        else:
            insert_node(node, treenode.left, depth+1)
    else:
        # go right
        if treenode.right is None:
            treenode.right = node
        else:
            insert_node(node, treenode.right, depth+1)

これは何をしますか？これにより、O（log n）時間でポイントを挿入できる検索可能なツリーが得られます。これは、リスト全体で O(n log n) を意味し、n の 2 乗よりも優れています。(500 万の対数底 2 は基本的に 23 です。したがって、n log n は 500 万 x 23 であり、500 万 x 500 万です!)

また、ターゲットを絞った検索を実行できることも意味します。ツリーは順序付けられているため、「近い」ポイントを探すのはかなり簡単です (@RootTwo からのウィキペディア リンクでアルゴリズムが提供されます)。

アドバイス

私のアドバイスは、必要に応じてリストをソートするコードを書くことです。書きやすく、手で確認しやすく、一度だけ作成する必要がある別のパスです。

リストを並べ替えたら、上で示した方法を試してください。それはあなたがしていたことに近く、理解しやすく、コーディングしやすいはずです。

Answer 3

リストを緯度 (n log(n)) で並べ替え、ポイントがほぼ均等に分布している場合、ポイントごとに 5 マイル以内に約 1000 ポイントになります (ナプキン計算、正確ではありません)。緯度が近い地点だけを見ると、実行時間は n^2 から n*log(n)+.0004n^2 になります。うまくいけば、これで十分に高速化されます。

Answer 4

この問題はVP ツリーで解決できます。これらにより、三角形の不等式に従うメトリックである距離を使用してデータをクエリできます。

kD ツリーに対する VP ツリーの大きな利点は、地理データを適切な 2D 空間に投影することを心配することなく、世界中のどこにでも地理データを盲目的に適用できることです。さらに、真の測地線距離を使用できます (測地線距離と投影の距離の違いについて心配する必要はありません)。

私のテストは次のとおりです。世界中でランダムかつ均一に 500 万ポイントを生成します。これらを VP ツリーに入れます。

すべてのポイントをループし、VP ツリーにクエリを実行して (0km, 10km] 離れた距離にある任意のネイバーを検索します (クエリ ポイントが見つからないように、このセットには 0km は含まれません)。そのようなネイバーがないポイントの数をカウントします。 (私の場合は 229573 です)。

VP ツリーの設定コスト = 5000000 * 20 距離の計算。

クエリのコスト = 5000000 * 23 距離計算。

セットアップとクエリの時間は 5 分 7 秒です。

私は距離を計算するためにGeographicLibで C++ を使用していますが、アルゴリズムはもちろん任意の言語で実装できます。これは GeographicLib の Python バージョンです。

補遺: このアプローチを実装する C++ コードをここに示します。

Answer 5

3 つのステップでアルゴリズムをやり直します。

1. 大円距離を使用し、誤差を 1% と仮定して、制限を 1.01*制限に等しくします。

2. 大圏距離をインライン関数としてコーディングします。このテストは高速です

3. いくつかの偽陽性が得られますが、vincenty でさらにテストできます。

Answer 6

これは最初のパスにすぎませんが、great_circle()の代わりに を使用vincinty()し、他のいくつかをクリーンアップすることで、これまでのところ半分にスピードアップしました。違いはここで説明されており、精度の損失は約0.17%です：

from geopy.point import Point
from geopy.distance import great_circle
import csv


class CustomGeoPoint(Point):
    def __init__(self, latitude, longitude):
        super(CustomGeoPoint, self).__init__(latitude, longitude)
        self.close_to_another_point = False


def isCloseToAnother(pointA, points):
    for pointB in points:
        dist = great_circle(pointA, pointB).miles
        if 0 < dist <= CLOSE_LIMIT:  # (0, close_limit]
            return True

    return False


    with open('geo_input.csv', 'r') as geo_csv:
        reader = csv.reader(geo_csv)
        next(reader, None)  # skip the headers

        geo_points = sorted(map(lambda x: CustomGeoPoint(x[0], x[1]), reader))

    with open('output.csv', 'w') as output:
        writer = csv.writer(output, delimiter=',', quoting=csv.QUOTE_ALL)
        writer.writerow(['Latitude', 'Longitude', 'CloseToAnotherPoint'])

        # for every point, look at every point until one is found within a mile
        for point in geo_points:
            point.close_to_another_point = isCloseToAnother(point, geo_points)
            writer.writerow([point.latitude, point.longitude,
                             point.close_to_another_point])

これをさらに改良していきます。

前：

$ time python geo.py

real    0m5.765s
user    0m5.675s
sys     0m0.048s

後：

$ time python geo.py

real    0m2.816s
user    0m2.716s
sys     0m0.041s