title: Python calculate lots of distances quickly
body: <p>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.</p>
  
  <p>I want to keep the pair combinations which are within 50km of each other. My current set-up is as follows</p>
  
  <pre><code>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 &lt;= 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)
  </code></pre>
  
  <p>This will obviously take hours and hours. Some possibilities I was thinking of:</p>
  
  <ul>
  <li>Use numpy to <strong>vectorise</strong> these calculations rather than looping through</li>
  <li>Use some kind of <strong>hashing</strong> to get a quick rough-cut off (all stores within 100km) and then only calculate accurate distances between those stores</li>
  <li>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 <strong>geodatabase</strong></li>
  <li>I can obviously try the <strong>haversine</strong> or project and use <strong>euclidean</strong> distances, however I am interested in using the most accurate measure possible</li>
  <li>Make use of <strong>parallel</strong> processing (however I was having a bit of difficulty coming up how to cut the list to still get all the relevant pairs).</li>
  </ul>
  
  <p><strong>Edit</strong>: I think geohashing is definitely needed here - an example <a href="https://pypi.python.org/pypi/geoindex://">from</a>:</p>
  
  <pre><code>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))
  </code></pre>
  
  <p>However, I would also like to vectorise (instead of loop) the distance calculations for the stores returned by the geo-hash.</p>
  
  <p><strong>Edit2: Pouria Hadjibagheri</strong> - I tried using lambda and map:</p>
  
  <pre><code># [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 &lt;= 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)
  </code></pre>
  
  <p>And they were all around <strong>61 seconds</strong> (I restricted number of stores to 2000 from 32,000). Perhaps I used map incorrectly?</p>