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

the sf package provides great approach working geographic features, can't figure out simple equivalent poly.counts function gistools package desires sp objects.

poly.counts computes number of points spatialpointsdataframe fall within polygons of spatialpolygonsdataframe , can used 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 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 can use poly.counts

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

add sp object:

sp_us_states$points.in.state <- points_in_states 

now i've finished i'd convert sf objects , visualise 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 perform operation without tedious conversion between sf , sp objects?

try 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 list indexes of points covered each state; lengths returns vector of lenghts of these vectors, or point count. warnings you'll see indicate although have geographic coordinates, underlying software assumes cartesian (which, case, not problematic; move projected coordinates if want rid of proper way)