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))
}