36,742 ポイントの入力があります。これは、距離行列の下三角を (ビンセンティ近似を使用して) 計算したい場合、36,742*36,741*0.5 = 1,349,974,563 距離を生成する必要があることを意味します。
50km以内のペアの組み合わせをキープしたいです。私の現在のセットアップは次のとおりです
shops= [[id,lat,lon]...]
def lower_triangle_mat(points):
for i in range(len(shops)-1):
for j in range(i+1,len(shops)):
yield [shops[i],shops[j]]
def return_stores_cutoff(points,cutoff_km=0):
below_cut = []
counter = 0
for x in lower_triangle_mat(points):
dist_km = vincenty(x[0][1:3],x[1][1:3]).km
counter += 1
if counter % 1000000 == 0:
print("%d out of %d" % (counter,(len(shops)*len(shops)-1*0.5)))
if dist_km <= cutoff_km:
below_cut.append([x[0][0],x[1][0],dist_km])
return below_cut
start = time.clock()
stores = return_stores_cutoff(points=shops,cutoff_km=50)
print(time.clock() - start)
これには明らかに何時間もかかります。私が考えていたいくつかの可能性：
- numpy を使用して、ループするのではなく、これらの計算をベクトル化します
- ある種のハッシングを使用して簡単な大まかなカットオフ (100 km 以内のすべての店舗) を取得し、それらの店舗間の正確な距離のみを計算します。
- ポイントをリストに保存する代わりに、四分木のようなものを使用しますが、実際の距離ではなく、近いポイントのランキングにのみ役立つと思います->だから、ある種のジオデータベースを推測します
- 私は明らかにハバーシンまたはプロジェクトを試してユークリッド距離を使用できますが、可能な限り最も正確な測定値を使用することに興味があります
- 並列処理を利用します（ただし、関連するすべてのペアを取得するためにリストを切り取る方法を考え出すのに少し苦労しました）。
編集：ここではジオハッシュが間違いなく必要だと思います-例：
from geoindex import GeoGridIndex, GeoPoint
geo_index = GeoGridIndex()
for _ in range(10000):
lat = random.random()*180 - 90
lng = random.random()*360 - 180
index.add_point(GeoPoint(lat, lng))
center_point = GeoPoint(37.7772448, -122.3955118)
for distance, point in index.get_nearest_points(center_point, 10, 'km'):
print("We found {0} in {1} km".format(point, distance))
ただし、geo-hash によって返された店舗の距離計算を (ループするのではなく) ベクトル化したいと考えています。
Edit2: Pouria Hadjibagheri - ラムダとマップを使ってみました:
# [B]: Mapping approach
lwr_tr_mat = ((shops[i],shops[j]) for i in range(len(shops)-1) for j in range(i+1,len(shops)))
func = lambda x: (x[0][0],x[1][0],vincenty(x[0],x[1]).km)
# Trying to see if conditional statements slow this down
func_cond = lambda x: (x[0][0],x[1][0],vincenty(x[0],x[1]).km) if vincenty(x[0],x[1]).km <= 50 else None
start = time.clock()
out_dist = list(map(func,lwr_tr_mat))
print(time.clock() - start)
start = time.clock()
out_dist = list(map(func_cond,lwr_tr_mat))
print(time.clock() - start)
そして、それらはすべて約61 秒でした (店舗数を 32,000 から 2000 に制限しました)。地図の使い方が間違っていたのでしょうか。