In 2018 Lorenzo Busetto published Speeding up spatial analyses by integrating sf and data.table: a test case on r-bloggers. The post addresses a common task: for each point in a set, count the other points within a given distance. It measured two methods up to 475,000 points and projected the 12-million-point case at 85 hours. We reproduced the post’s measurements on one machine, then computed the same counts on the same points with rgSpatial. The exact count, at 12 million points and 2,000 m, took 19.4 minutes.
The question appears across the field: exposure within a radius, neighbourhood density, distance-band weights for spatial regression, sample points that share a footprint. The answer is a single number per point, and the work behind it grows with the square of the point density, since every point must be tested against every point that might be near it.
The post starts from the standard sf formulation — buffer every point by
d, then st_intersects the buffers against the points — measures
it, and then improves on it with a function that partitions the extent into tiles
and processes each tile with data.table, keeping the buffered sets small and allowing
the tiles to run in parallel. It times both methods up to 475,000 points, at 500 m
and at 2,000 m, and extrapolates to 12 million points.
The post’s projections for 12 million points are 15 hours at 500 m and 85 hours at 2,000 m for the baseline, 2.5 and 15.2 hours for the tiled function on one core, and 0.57 and 3.18 hours on six cores. The post is a clear piece of work and its tiled function is a real improvement over the baseline. We used it as our reference.
# buffer every point by maxdist, then intersect the # buffers with the points: one list of neighbours per point pts_buf <- sf::st_buffer(pts, maxdist) int <- sf::st_intersects(pts_buf, pts)
# sizes: 25,000 · 100,000 · 475,000 points # radii: 500 m and 2,000 m # points: uniform random on a 100 km square # then: extrapolated to 12,000,000 points, # serial and on six cores
points_in_distance(), the post’s contribution as published, unchangedpoints_in_distance <- function(in_pts, maxdist, ncuts = 10) { require(data.table) require(sf) pts <- data.table(in_pts) pts <- pts[, or_id := 1:dim(in_pts)[1]] range_x <- range(pts$x) limits_x <-(range_x[1] + (0:ncuts)*(range_x[2] - range_x[1])/ncuts) range_y <- range(pts$y) limits_y <- range_y[1] + (0:ncuts)*(range_y[2] - range_y[1])/ncuts pts[, `:=`(xcut = as.integer(cut(x, ncuts, labels = 1:ncuts)), ycut = as.integer(cut(y, ncuts, labels = 1:ncuts)))] %>% setkey(x) results <- list() count <- 0 for (cutx in seq_len(ncuts)) { min_x_comp <- ifelse(cutx == 1, limits_x[cutx], (limits_x[cutx] - maxdist)) max_x_comp <- ifelse(cutx == ncuts, limits_x[cutx + 1], (limits_x[cutx + 1] + maxdist)) subpts_x <- pts[x >= min_x_comp & x < max_x_comp] %>% setkey(y) for (cuty in seq_len(ncuts)) { count <- count + 1 min_y_comp <- ifelse(cuty == 1, limits_y[cuty], (limits_y[cuty] - maxdist)) max_y_comp <- ifelse(cuty == ncuts, limits_y[cuty + 1], (limits_y[cuty + 1] + maxdist)) subpts_comp <- subpts_x[y >= min_y_comp & y < max_y_comp] subpts_buf <- subpts_comp[ycut == cuty & xcut == cutx] %>% sf::st_as_sf() %>% sf::st_buffer(maxdist) subpts_comp <- sf::st_as_sf(subpts_comp) inters <- sf::st_intersects(subpts_buf, subpts_comp) results[[count]] <- data.table( id = subpts_buf$or_id, int_ids = lapply(inters, FUN = function(x) subpts_comp$or_id[x])) } } data.table::rbindlist(results) }
st_intersects call sees one tile and its margin, a small fraction of the points, and buffer geometry exists only for that tile at a time.The function returns the neighbour lists, as the baseline does. Both methods therefore hold every pair in memory at once; section 06 returns to this. The answer is the buffer-polygon answer in both cases (section 05).
We ran both of the post’s methods as published — the buffer-and-intersect baseline and the tiled data.table function — on seeded uniform random points over a 100 km square, generated the way the post generates them. Then the same points, the same radius, through one rgSpatial verb.
Bars to scale. Each method counts the same 279,309,308 neighbour pairs. rgSpatial from R is 114× faster than the post’s function and 339× faster than the baseline.
The task is a distance join followed by a count per point. rgSpatial exposes it as one verb that performs the join and the aggregate in a single pass. The R and Python bindings call the same engine and return the same counts.
# every point, with how many points lie within 2 km of it pts <- rg_points_from_xy("x", "y", data = rg_read("pts.parquet")) counts <- rg_join_agg(pts, pts, n = c("id", "count"), join = "dwithin", distance = 2000)
# every point, with how many points lie within 2 km of it tbl = rg.read("pts.parquet") pts = rg.GeoDataFrame({"id": tbl["id"]}, geometry=rg.points_from_xy(tbl["x"], tbl["y"])) counts = pts.sjoin_agg(pts, predicate="dwithin", distance=2000, n=("id", "count"))
No tiling, buffering or chunk loop is written by the user, and no pair table is produced. The result is a frame with one row per point and its count, carrying the CRS.
The post stopped measuring at 475,000 points and extrapolated to 12 million; section 06 explains why the larger case is out of reach for any method that materialises the pairs. We ran it at both radii, from R and from Python. The post’s projections are reproduced as published and marked as estimates.
Bars to scale within each chart. Relative to the post’s six-core projection, rgSpatial from R is 24× faster at 500 m and 9.8× at 2,000 m, on one desktop.
“Points within 2,000 m” means points inside a circle of radius 2,000 m. A buffer is
a polygon with a fixed number of segments per quarter turn, thirty by default in sf,
inscribed in that circle. A neighbour in the outer 0.03% of the radius falls outside
the polygon and is not counted. This is a property of the idiom, not of the post; the
idiom is common because sf’s exact predicate, st_is_within_distance, is
too slow to use at this scale. Both of the post’s methods are built on the buffer and
both inherit the approximation.
rgSpatial’s dwithin evaluates the distance itself. We checked it
against st_is_within_distance, point for point, at the sizes where that
call completes in reasonable time.
At 475,000 points and 2 km, one pair in 2,200 is missed by the buffer method. The effect is small enough to pass unnoticed and large enough to matter where a count feeds a threshold, a model or a regulatory test.
The dataframe libraries answer a distance join by building an index, finding the candidate pairs, and returning the list of pairs for the user to aggregate. At 12 million points and 2 km that list has 178 billion rows, which is why the post could extrapolate the case but not run it. rgSpatial computes the aggregate inside the join.
A spatial index is built over the points in one pass. It is the last large allocation the computation makes.
Each point queries the index for candidates within its radius and tests the exact distance to each. Native points are tested on their coordinates without decoding or copying.
Each hit updates a running accumulator for its point (a count, a sum, a minimum) and is discarded. The accumulators are the result; no pair is stored, sorted or grouped.
Points are processed in chunks across all cores, each chunk with its own accumulators, merged once at the end. Memory use is independent of the number of pairs.
This is the plan a query optimiser produces for a GROUP BY over a join, the
aggregate pushed into the join’s probe loop, offered here as a dataframe verb. The same
fused aggregate serves rgSpatial’s point-in-polygon joins. When the pairs themselves are
wanted, the ordinary join returns them.
Compute seconds on one machine, reads excluded for every method. Blank cells are sizes the post did not run; italic cells are the post’s published extrapolations.
| points | radius | post’s baseline | post’s function | sf exact | rgSpatial R | rgSpatial Python | neighbour pairs |
|---|---|---|---|---|---|---|---|
| 25,000 | 500 m | 0.95 s | 1.30 s | 11.7 s | 0.03 s | 0.007 s | 73,782 |
| 25,000 | 2,000 m | 1.08 s | 1.36 s | 11.9 s | 0.05 s | 0.005 s | 794,930 |
| 100,000 | 500 m | 4.9 s | 4.5 s | 359 s | 0.05 s | 0.021 s | 880,046 |
| 100,000 | 2,000 m | 8.1 s | 6.0 s | 352 s | 0.05 s | 0.029 s | 12,457,986 |
| 475,000 | 500 m | 31.9 s | 23.6 s | 0.17 s | 0.14 s | 18,127,930 | |
| 475,000 | 2,000 m | 156 s | 52.3 s | 0.46 s | 0.42 s | 279,309,308 | |
| 2,000,000 | 500 m | 1.85 s | 1.71 s | 314,866,686 | |||
| 2,000,000 | 2,000 m | 13.9 s | 13.9 s | 4,943,197,648 | |||
| 12,000,000 | 500 m | 15 h | 2.5 h · 0.57 h on 6 cores | 86 s | 88 s | 11,273,747,962 | |
| 12,000,000 | 2,000 m | 85 h | 15.2 h · 3.18 h on 6 cores | 19.4 min | 22.0 min | 177,896,999,236 |
rgSpatial’s counts equal sf’s exact st_is_within_distance on every point at
25,000 and 100,000; above that the exact call was not run. The post’s methods return the
buffer approximation at every size (section 05).
Data. Uniform random points on a 100 km × 100 km square, generated the way the post
generates them, with one change: a fixed seed, so every method sees the same points.
Sizes 25,000 · 100,000 · 475,000 · 2,000,000 · 12,000,000; radii 500 m and 2,000 m.
The post’s code. Lorenzo Busetto, Speeding up spatial analyses by
integrating sf and data.table: a test case, r-bloggers, February 2018. Run as published — the buffer-and-intersect
baseline and the tiled points_in_distance() function — at the sizes the post ran.
Its 12-million figures are the post’s own projections, quoted as printed and marked as such.
Correctness. Every rgSpatial count is compared point-for-point with sf’s exact
st_is_within_distance at 25,000 and 100,000 points, and with the post’s
methods at every size the post ran. A wrong answer would disqualify the timing.
Machine. One desktop: AMD Ryzen 9 7950X3D, 16 cores, 64 GB, Windows 11. R 4.5.2,
sf 1.1.2, data.table 1.18.2.1, Python 3.13.5, rg 0.5.1. Timings are compute only; file
reads are excluded for every method. Measured 2026-09-07.
We will run your workload on your data, in R or Python, and report the result alongside a reference computation.
Request Access