title: (Spatial) Efficient way of finding all points within X meters of a point?
body: <p>I have a large spatial dataset (12M rows). The geometries are points on a map. For each row in the dataset, I'd like to find all the points that are within 500 meters of that point. </p>
  
  <p>In r, using sf, I've been trying to do this by parallel looping through each row and running st_buffer and st_intersects, then saving the result as a list in a key-value format (the key being the origin point, the values being the neighbors). </p>
  
  <p>The issue is that the dataset is too large. Even when parallelizing to upwards of 60 cores, the operation takes too long (>1 week and usually crashes). </p>
  
  <p>What are the alternatives to this brute-force approach? Is it possible to build indexes using sf? Perhaps push the operation to an external database? </p>
  
  <p>Reprex:</p>
  
  <pre><code>library(sf)
  library(tidyverse)
  library(parallel)
  library(foreach)
  
  
  # example data, convert to decimal:
  nc &lt;- st_read(system.file("shape/nc.shp", package="sf")) %&gt;% st_transform(32618)
  # expand the data a a bit to make the example more interesting:
  nc &lt;- rbind(nc,nc,nc)
  nc &lt;- nc %&gt;% mutate(Id = row_number())
  
  
  ## can run in parallel if desired:
  # num_cores &lt;- parallel::detectCores()-2
  # cl &lt;- makeSOCKcluster(num_cores)
  # registerDoSNOW(cl)
  
  # or just run in sequence:
  registerDoSEQ()
  
  neighbors &lt;- foreach(ii = 1:nrow(nc)
                        , .verbose = FALSE
                        , .errorhandling = "pass") %dopar% {
  
                          l = 500 # 500 meters
  
                          # isolate the row as the origin point:
                          row_interest &lt;- filter(nc, row_number()==ii)
  
                          # create the buffer:
                          buffer &lt;- row_interest %&gt;% st_buffer(dist = l)
  
                          # extract the row numbers of the neighbors
                          comps_idx &lt;- suppressMessages(st_intersects(buffer, nc))[[1]]
  
                          # get all the neighbors:
                          comps &lt;- nc %&gt;% filter(row_number() %in% comps_idx)
  
                          # remove the geometry:
                          comps &lt;- comps %&gt;% st_set_geometry(NULL)
  
                          # flow control in case there are no neibors:
                          if(nrow(comps)&gt;0) {
                            comps$Origin_Key &lt;- row_interest$Id
                          } else {
                            comps &lt;- data_frame("lat" = NA_integer_,"lon" = NA_integer_, "bbl" = row_interest$bbl)
                            comps$Origin_Key &lt;- row_interest$Id
                          }
  
  
                          return(comps)
                        }
  
  closeAllConnections()
  
  length(neighbors)==nrow(nc)
  [1] TRUE
  </code></pre>