Covers vs intersects

# Packages
library(sf)
#> Linking to GEOS 3.13.0, GDAL 3.10.1, PROJ 9.5.1; sf_use_s2() is TRUE

# Let's work on the unit square
unit_square <- st_as_sfc(st_bbox(c(xmin = 0, xmax = 1, ymin = 0, ymax = 1)))

# Let's divide the unit square into 1600 (regular) cells
cells <- st_as_sf(st_make_grid(unit_square, cellsize = c(0.025, 0.025)))

# The value of each cell is equal to x_centroid ^ 2 + y_centroid ^ 2
cells$vals <- st_centroid(cells) |> 
  st_coordinates() |>
  apply(1, function(x) x[1]^2 + x[2]^2)

# Check the output
plot(cells)

# Let's sample 2000 points into the unit square
set.seed(1)
points <- st_as_sf(st_sample(unit_square, 2000))

# 1. Determine which cell covers which point
match_cells_points <- st_covers(cells, points)

# 2. Determine how many points fall in each cell
times_pts <- lengths(match_cells_points)

# 3. Replicate each cell value for each point that falls in the cell
points$vals <- NA_real_
points[unlist(match_cells_points), "vals"] <- rep(cells$vals, times = times_pts)

# Check the output
plot(points, pch = 16)

# This is slightly faster than a st_intersects approach
library(bench)

# Let's consider a wild example with 1'000'000 cells
cells <- st_make_grid(unit_square, cellsize = c(0.001, 0.001))

# and 100'000 points
points <- st_sample(unit_square, 100000)

mark(
  "intersects" = {
    st_intersects(points, cells)
  }, 
  "covers" = {
    st_covers(cells, points) |> lengths()
  }, 
  check = FALSE
)
#> Warning: Some expressions had a GC in every iteration; so filtering is
#> disabled.
#> # A tibble: 2 × 6
#>   expression      min   median `itr/sec` mem_alloc `gc/sec`
#>   <bch:expr> <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
#> 1 intersects   10.74s   10.74s    0.0931      50MB    0.745
#> 2 covers        2.93s    2.93s    0.342     19.8MB    2.05

^{Created on 2025-04-07 with reprex v2.1.1.9000}

agila5/st_covers.md