How do I find the overlap between two shapefiles?

Question

I have two shapefiles (sf), one with polygons and one with points. As output I want a df showing which points fall within which polygons, something like this:

polygon   overlap  geometry
polygon1  point34  c(3478,234872)
polygon1  point56  c(23423,234982)
polygon2  point23  c(23498,2334)
polygon3  point45  c(872348,23847)
polygon3  point87  c(234982,1237)
polygon3  point88  c(234873,2873)

I assume I'll have to do something with st_intersection() but up to now I did not manage to get the desired output.

Answer 1

After fiddling around I came up with this solution, but I'm pretty sure it is not the most elegant. x and y are shapefiles, x with points and y with polygons.

count_overlap <- function(x, y){
  f1 <- function(z){
    r <- st_intersects(x,y[z,])
    return(r)
    }
  l1 <- c(1:nrow(y))
  l2 <- lapply(l1, f1)
  l3 <- lapply(l2, unlist)
  r <- sapply(l3, sum)
  y$overlap <- r
  return(y)
}

The result is the original y sf/dataframe with an added column called 'overlap' that shows the counts of points from x that fall within the polygon. Not exactly what I asked for in the question but a good outcome for me personally.