Faster st_is_within_distance for points to lines

I am looking for a fast way to run a `st_is_within_distance` (or equivalent) in R, to find all the lines (y) within a given distance of each point (x). (This would be an equivalent of ArcGIS's "neartable", for which we do not seem to have an easy equivalent in R's GIS packages, although there are packages for points-to-points, see e.g. Is there an R equivalent of NEAR in ESRI ArcGIS?) The R help for `?st_is_within_distance` contains this note:

For most predicates, a spatial index is built on argument x; see http://r-spatial.org/r/2017/06/22/spatial-index.html. Specifically, st_intersects, st_disjoint, st_touches st_crosses, st_within, st_contains, st_contains_properly, st_overlaps, st_equals, st_covers and st_covered_by all build spatial indexes for more efficient geometry calculations. st_relate, st_equals_exact, and st_is_within_distance do not.

So my interpretation is that `st_is_within_distance` is going to be no faster than calculating all possible pairs of distances and taking the ones within the given cutoff. From my testing, I find this is true.

So I had thought, as a way to improve on the speed of `st_is_within_distance`, to use an approach based on a spatial buffer: buffer the lines out to the distance I want (15km in my case), then use `st_within` to see if the points fall into the buffers (and pick which buffers). Conceptually this is similar to what `st_is_within_distance` is doing, although this gives slightly different answers, probably because the buffer is only an approximation to a "true buffer". Worrying about this turns out to be not important though, because, much to my surprise -- and despite using a spatial index -- this approach is SLOWER than calculating all possible distances! (about 1.7 times slower in my example)

Questions: (1) Is there a way (in theory and/or in R implementation) to use a spatial index to obtain a speed improvement for finding all lines within a given distance of a point? Or is calculating all possible distances the only sure way to accomplish this? (Note that in sf, `st_nearest_feature` uses a spatial index and is very fast, but that only obtains the closest line to the point, not all lines within a given radius of the point. This suggests to me that a spatial index could perhaps speed up the `st_is_within_distance` call.)

(2) Why, despite using a spatial index, is the `st_within` approach slower than calculating all possible distances? I am looking for any intuition to help understand why I was wrong.

• How bout buffering the points and then use st_intersects? I am assuming you have more points than lines. Sep 24 '20 at 17:09
• @TimSalabim wow, your approach is way faster (14 seconds from buffer points and st_intersects vs 200 seconds from st_is_within_distance). yes, I'm working with variable number of points (20000 now, up to 2 million later), fixed number of lines (7000 lines). is the speed diff because st_within was preparing the polygons rather than the points when I was buffering the lines? Sep 24 '20 at 20:59
• I believe so yes. Sep 25 '20 at 5:58