Equivalent of `poly.counts` to count lat/long pairs falling inside of polygons with the sf package

Question

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?

Answer 1

Try the following:

sf_obs = st_as_sf(observations_tibble, coords = c("long", "lat"), 
    crs = st_crs(sf_us_states))
lengths(st_covers(sf_us_states, sf_obs))
# check:
summary(points_in_states - lengths(st_covers(sf_us_states, sf_obs)))

st_covers returns a list with the indexes of points covered by each state; lengths returns the vector of the lenghts of these vectors, or the point count. The warnings you'll see indicate that although you have geographic coordinates, the underlying software assumes they are cartesian (which, for this case, will be most likely not problematic; move to projected coordinates if you want to get rid of it the proper way)