Python calculate lots of distances quickly

I have an input of 36,742 points which means if I wanted to calculate the lower triangle of a distance matrix (using the vincenty approximation) I would need to generate 36,742*36,741*0.5 = 1,349,974,563 distances.

I want to keep the pair combinations which are within 50km of each other. My current set-up is as follows

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)

This will obviously take hours and hours. Some possibilities I was thinking of:

  • Use numpy to vectorise these calculations rather than looping through
  • Use some kind of hashing to get a quick rough-cut off (all stores within 100km) and then only calculate accurate distances between those stores
  • Instead of storing the points in a list use something like a quad-tree but I think that only helps with the ranking of close points rather than actual distance -> so I guess some kind of geodatabase
  • I can obviously try the haversine or project and use euclidean distances, however I am interested in using the most accurate measure possible
  • Make use of parallel processing (however I was having a bit of difficulty coming up how to cut the list to still get all the relevant pairs).
  • Edit : I think geohashing is definitely needed here - an example from:

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

    However, I would also like to vectorise (instead of loop) the distance calculations for the stores returned by the geo-hash.

    Edit2: Pouria Hadjibagheri - I tried using lambda and 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)
    

    And they were all around 61 seconds (I restricted number of stores to 2000 from 32,000). Perhaps I used map incorrectly?


    This sounds like a classic use case for kD trees.

    If you first transform your points into Euclidean space then you can use the query_pairs method of scipy.spatial.cKDTree :

    from scipy.spatial import cKDTree
    
    tree = cKDTree(data)
    # where data is (nshops, ndim) containing the Euclidean coordinates of each shop
    # in units of km
    
    pairs = tree.query_pairs(50, p=2)   # 50km radius, L2 (Euclidean) norm
    

    pairs will be a set of (i, j) tuples corresponding to the row indices of pairs of shops that are ≤50km from each other.


    The output of tree.sparse_distance_matrix is a scipy.sparse.dok_matrix . Since the matrix will be symmetric and you're only interested in unique row/column pairs, you could use scipy.sparse.tril to zero out the upper triangle, giving you a scipy.sparse.coo_matrix . From there you can access the nonzero row and column indices and their corresponding distance values via the .row , .col and .data attributes:

    from scipy import sparse
    
    tree_dist = tree.sparse_distance_matrix(tree, max_distance=10000, p=2)
    udist = sparse.tril(tree_dist, k=-1)    # zero the main diagonal
    ridx = udist.row    # row indices
    cidx = udist.col    # column indices
    dist = udist.data   # distance values
    

    Have you tried mapping entire arrays and functions instead of iterating through them? An example would be as follows:

    from numpy.random import rand
    
    my_array = rand(int(5e7), 1)  # An array of 50,000,000 random numbers in double.
    

    Now what is normally done is:

    squared_list_iter = [value**2 for value in my_array]
    

    Which of course works, but is optimally invalid.

    The alternative would be to map the array with a function. This is done as follows:

    func = lambda x: x**2  # Here is what I want to do on my array.
    
    squared_list_map = map(func, test)  # Here I am doing it!
    

    Now, one might ask, how is this any different, or even better for that matter? Since now we have added a call to a function, too! Here is your answer:

    For the former solution (via iteration):

    1 loop: 1.11 minutes.
    

    Compared to the latter solution (mapping):

    500 loop, on average 560 ns. 
    

    Simultaneous conversion of a map() to list by list(map(my_list)) would increase the time by a factor of 10 to approximately 500 ms .

    You choose!


    "Use some kind of hashing to get a quick rough-cut off (all stores within 100km) and then only calculate accurate distances between those stores" I think this might be better called gridding. So first make a dict, with a set of coords as the key and put each shop in a 50km bucket near that point. then when you are calculating distances, you only look in nearby buckets, rather than iterate through each shop in the whole universe

    链接地址: http://www.djcxy.com/p/90274.html

    上一篇: 在Internet Explorer和多个监视器中的Window.screen

    下一篇: Python快速计算大量距离