我有一个包含 36,742 个点的输入数据，这意味着如果我想计算距离矩阵（使用 Vincenty 近似）的下三角，则需要生成 36,742*36,741*0.5 = 1,349,974,563 个距离。
我希望保留在彼此之间距离在 50 公里以内的配对组合。我的当前设置如下：
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公里内的所有商店），然后只计算这些商店之间的准确距离
- 不要将点存储在列表中，而是使用类似于四叉树的东西，但我认为这只有助于接近点的排名，而不是实际距离 -> 所以我猜测需要某种地理数据库
- 我可以尝试使用haversine或投影并使用欧几里得距离，但我有兴趣使用可能最准确的度量标准
- 利用并行处理（但我发现难以想出如何切割列表以仍然获得所有相关的对）。
编辑：我认为这里绝对需要geohashing - 例如此处的示例：
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))
然而，我希望你能将地理哈希返回的商店的距离计算向量化（而不是循环）。
Edit2: Pouria Hadjibagheri - 我尝试使用lambda和map：
# [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个限制为2,000个）。也许我使用地图不正确？