Using st_intersects to classify points in a multi polygon

I have a list of lat/lons and want to determine which county they belong to, by intersecting them with multi polygon dataset using sf::st_intersect. I can successfully perform the intersection and get a binary list, but am struggling with the final step of converting that to a data frame where each point is associated with the polygon information.

Here is a reproducible example, largely borrowed from an answer to a related question:

``````library(sf)
nc <- st_read(system.file("shape/nc.shp", package="sf")) #sample data
nc <- st_transform(nc, 3857)
pts <- st_jitter(st_sample(nc, 1000), factor=0.2) #sample points

# intersect data
intersect_list <- st_intersects(pts, nc)
``````

intersect_list is a list; the second column is a binary integer (presumably indicating whether or not the point intersected with a polygon), and the third column is a list containing integers (my guess that this is the id of the specific polygon that was intersected). At this point, I'd like to get the CNTY_ID and FIPS associated with each point, but I'm not sure how to complete this step. In sp, I would have done something like:

``````df <- pts #make a data frame with the points
intersect_list <- over(pts, nc)
df\$countyID <- intersect_list\$CNTY_ID
df\$fips <- intersect_list\$FIPS
``````

The result `intersect_list` is a prettied-up R list where each element represents a point in `pts` and the contents of the element are the row numbers of `nc` that the point intersects:

``````> intersect_list
Sparse geometry binary predicate list of length 1000, where the
predicate was `intersects'
first 10 elements:
1: 32
2: 50
3: 51
4: 57
5: (empty)
6: (empty)
7: (empty)
``````

So the first point intersects with `nc[30,]` only, the fifth point doesn't intersect any of `nc` polygons (so returns an empty element). Since polygons can in general overlap, its possible for a single point to intersect multiple polygons, and then you'd get a list like this one I've made up:

``````> intersect_list
Sparse geometry binary predicate list of length 1000, where the
predicate was `intersects'
first 10 elements:
1: 32 33
2: 50 51 52
3: (empty)
4: (empty)
5: (empty)
``````

showing that point 1 is in polygons 32 and 33, and point 2 is in polygons 50, 51 and 52.

At this point in your analysis you have to decide what to do with non-matched points and multiply matched points. You can test for these by getting the length of each element in the list using the `lengths` (not `length`) function. This returns a vector of the number of elements. For example the first 7 elements gives me:

``````> lengths(intersect_list)[1:7]
[1] 1 1 1 1 0 0 0
``````

showing the first four points are only in one polygon and then there's some in no polygons (hence "empty"). You can do some tests and in general do anything.

Another approach is to make sure you have two spatial data frames and use `st_join`, which will do all the intersecting and matching and return a data frame joined on a geometry operation:

``````> pts = st_as_sf(pts) # make it dataframey
> st_join(pts, nc)
Simple feature collection with 1000 features and 14 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -9512554 ymin: 3815132 xmax: -8259388 ymax: 4562025
Projected CRS: WGS 84 / Pseudo-Mercator
First 10 features:
AREA PERIMETER CNTY_ CNTY_ID     NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
1  0.059     1.319  1927    1927 Mitchell 37121  37121       61   671     0
2  0.134     1.590  1980    1980    Rowan 37159  37159       80  4606     3
3  0.168     1.791  1984    1984     Pitt 37147  37147       74  5094    14
4  0.203     3.197  2004    2004 Beaufort 37013  37013        7  2692     7
5     NA        NA    NA      NA     <NA>  <NA>     NA       NA    NA    NA
6     NA        NA    NA      NA     <NA>  <NA>     NA       NA    NA    NA
7     NA        NA    NA      NA     <NA>  <NA>     NA       NA    NA    NA
8  0.241     2.214  2083    2083  Sampson 37163  37163       82  3025     4
9     NA        NA    NA      NA     <NA>  <NA>     NA       NA    NA    NA
10 0.091     1.470  2068    2068   Gaston 37071  37071       36  9014    11
NWBIR74 BIR79 SID79 NWBIR79                 geometry
1        1   919     2       4 POINT (-9159955 4316755)
2     1057  6427     8    1504 POINT (-8983465 4251687)
3     2620  6635    11    3059 POINT (-8601253 4224719)
4     1131  2909     4    1163 POINT (-8535427 4240341)
5       NA    NA    NA      NA POINT (-8838087 3905710)
6       NA    NA    NA      NA POINT (-8576943 4020740)
7       NA    NA    NA      NA POINT (-8456819 4233280)
8     1396  3447     4    1524 POINT (-8702113 4178122)
9       NA    NA    NA      NA POINT (-8816557 4083780)
10    1523 11455    26    2194 POINT (-9057212 4219754)
``````

This is a table with POINT geometry with attributes taken from the intersecting POLYGON geometry. You'll get NA when a point falls in no polygons, and (I think) multiple matches when a point is in multiple polygons - check this with an example or it should be in the docs.