Comparison of spatial index performance in {sf} and {terra}

spatial_indexes.R

library("sf")
library("terra")

n = c(10, 50, 100, 150, 200)
t = matrix(NA, length(n), 3)

for (i in seq_along(n)) {

  r = rast(xmin = 0, xmax = n[i], ymin = 0, ymax = n[i], ncol = n[i], nrow = n[i],
           crs = "epsg:3857") # plannar CRS
  v1 = as.polygons(r)
  v2 = shift(v1, 0.5, 0.5)
  g1 = st_as_sf(v1)
  g2 = st_as_sf(v2)
  t[i, 1] = system.time(st_intersects(g1, g2))[[3]]
  # arguments: (SpatVector v, std::string relation, bool prepared, bool index)
  t[i, 2] = system.time(v1@ptr$relate_between(v2@ptr, "intersects", TRUE, TRUE))[[3]]
  t[i, 3] = system.time(v1@ptr$relate_between(v2@ptr, "intersects", TRUE, FALSE))[[3]]

}

plot(NULL, xlim = c(0, max(n)^2), ylim = c(0, 1.1 * max(t)),
     ylab = "Time [s]", xlab = "Polygons", main = "Geometry intersects")
lines(n^2, t[, 1], type = "b", col = 1)
lines(n^2, t[, 2], type = "b", col = 2)
lines(n^2, t[, 3], type = "b", col = 3)
legend("topleft", horiz = TRUE, lty = 1, col = 1:3, bty = "n",
       legend = c("sf", "terra (with index)", "terra (old)"))