Counting points in polygons with sf package of R?

I'm quite new to R and especially GIS with R. So recently I had to proceed to a quite simple analysis: counting health centers (points) per administrative unit (polygons). Quite common process in arcgis or qgis but I didn't found any similar tool in R, at least with the "sf" package that I'm trying to work with exclusively for now (for an easier learning of the way this package works).

So I googled it and found some people with the same problem but I didn't understand most of the methods they were suggesting, it looked oddly complicated for such a simple operation to me. So I created a function to solve my problem, which works now. But I would like to know in order to get better at this:

Is there an easier way that I missed?

Here's the code below for the function which works simply by: selecting polygons one by one, intersecting it with the point layer, counting the size of the extracted points, adding that count to a vector, repeating with the next polygon... And in the end, binding that "count" vector to my initial polygon dataframe.

CountPointsInPolygons <- function(pts, polygons){
  countPts = c()
  for (i in 1:nrow(polygons)) {
    polySelect <- polygons[i,]
    pts2 <- st_intersection(pts, polySelect)
    countPts[i] = nrow(pts2)


  }

  return(cbind(polygons,countPts))
}

Answer

It's faster and easier to use st_intersects, no need to loop. The output is a bit obscure, essentially a classed list of feature IDs intersected, so we get the lengths() which is the number of points inside each feature.

First polygon has six points in it, in this example.

  library(sf)
#> Linking to GEOS 3.7.0, GDAL 2.4.0, PROJ 5.2.0

poly <- read_sf(system.file("shape/nc.shp", package = "sf"))
set.seed(77)
pts <- st_sample(poly, size = 750)

## add point count to each polygon
(poly$pt_count <- lengths(st_intersects(poly, pts)))
#>   [1]  6  8 12  3  9  2  3  4  9  9 10  5  6  5  4  9  4  7  4  1  4  4  4
#>  [24]  5  4 10 10 11  9  1  9  3  4  3  4  5 11  5  6  4  3  7 10  8  3  5
#>  [47]  7  9  6 14  9  9  6 11  6  9 12  7  5  7  9 11 10  6  6 12 12  6  4
#>  [70]  7  9  7  4  3  8  5  3  6 19  5  7 12  5 12  9 10  8 23  6  5  5  7

#>  [93] 12  9 12 16 18 17  2 14

plot(poly[1, 1], reset = FALSE, col = "grey")
plot(pts, add = TRUE)

^{Created on 2019-05-23 by the reprex package (v0.3.0)}