The sf package provides a great approach to working with geographic features, but I can't figure out a simple equivalent to the poly.counts function from GISTools package which desires sp objects.

poly.counts computes the number of points from a SpatialPointsDataFrame fall within the polygons of a SpatialPolygonsDataFrame and can be used as follows:

Data

## Libraries
library("GISTools")
library("tidyverse")
library("sf")
library("sp")
library("rgdal")
## Obtain shapefiles
download.file(url = "https://www2.census.gov/geo/tiger/TIGER2016/STATE/tl_2016_us_state.zip", destfile = "data-raw/states.zip")
unzip(zipfile = "data-raw/states.zip", exdir = "data-raw/states")
sf_us_states <- read_sf("data-raw/states")

## Our observations:
observations_tibble <- tribble(
  ~lat, ~long,
31.968599,  -99.901813,
35.263266,  -80.854385,
35.149534,  -90.04898,
41.897547,  -84.037166,
34.596759,  -86.965563,
42.652579,  -73.756232,
43.670406,  -93.575858
)

Calculate points per polygon

I generate both my sp objects:

sp_us_states <- as(sf_us_states, "Spatial")

observations_spdf <- observations_tibble %>%
  select(long, lat) %>% # SPDF want long, lat pairs
  SpatialPointsDataFrame(coords = .,
                         data = .,
                         proj4string = sp_us_states@proj4string)

Now I can use poly.counts

points_in_states <-
  poly.counts(pts = observations_spdf, polys = sp_us_states)

Add this into the sp object:

sp_us_states$points.in.state <- points_in_states

Now I've finished I'd convert back to sf objects and could visualise as follows:

library("leaflet")
updated_sf <- st_as_sf(sp_us_states)
updated_sf %>%
  filter(points.in.state > 0) %>%
  leaflet() %>%
  addPolygons() %>%
  addCircleMarkers(
    data = observations_tibble
  )

Question

Can I perform this operation without tedious conversion between sf and sp objects?