Tengo una entrada de 36,742 puntos, lo que significa que si quisiera calcular el triángulo inferior de una matriz de distancia (usando la aproximación de vincenty) necesitaría generar 36,742 * 36,741 * 0.5 = 1,349,974,563 distancias.

Quiero mantener las combinaciones de pares que están dentro de los 50 km el uno del otro. Mi configuración actual es la siguiente

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)

Obviamente, esto llevará horas y horas. Algunas posibilidades en las que estaba pensando:

Use numpy paravectorizar estos cálculos en lugar de recorrerlosUsa algún tipo dehashing para obtener un corte rápido (todas las tiendas dentro de los 100 km) y luego solo calcular distancias precisas entre esas tiendasEn lugar de almacenar los puntos en una lista, use algo como un árbol cuádruple, pero creo que eso solo ayuda con la clasificación de los puntos cercanos en lugar de la distancia real -> así que supongo que algún tipo degeodatabaseObviamente puedo probar elHaversine o proyectar y usareuclidiana distancias, sin embargo, estoy interesado en utilizar la medida más precisa posibleHacer uso deparalela procesamiento (sin embargo, me estaba costando un poco encontrar la forma de cortar la lista para obtener todos los pares relevantes).

Editar: Creo que definitivamente se necesita geohashing aquí - un ejemplode:

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))

Sin embargo, también me gustaría vectorizar (en lugar de bucle) los cálculos de distancia para las tiendas devueltas por el geo-hash.

Edit2: Pouria Hadjibagheri - Intenté usar lambda y mapa:

# [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)

Y estaban por todos lados61 segundos (Limité el número de tiendas a 2000 de 32,000). ¿Quizás usé el mapa incorrectamente?