How to keep NA values when using st_intersection in R

Question

I have a polygon layer representing different site treatments and a point layer representing plots taken within those treatments and several control plots outside of the treatments. I need to extract the treatment information from the polygons and add it to the point layer, but there are several control plots that are outside of the polygons that I need to keep. I have been using the sf package in R. Is there a way to use st_intersection and keep NA values for the points that don't intersect the polygon layer, or do I need a different strategy?

Answer 1

You may need to use st_intersects to test intersection first. Example, let's use a couple of polygons from the nc data set:

> library(sf)
> nc = st_read(system.file("shape/nc.shp", package="sf"))[c(23,42),]
Reading layer `nc' from data source `/nobackup/rowlings/RLibrary/R/x86_64-pc-linux-gnu-library/3.6/sf/shape/nc.shp' using driver `ESRI Shapefile'
Simple feature collection with 100 features and 14 fields
geometry type:  MULTIPOLYGON
dimension:      XY
bbox:           xmin: -84.32385 ymin: 33.88199 xmax: -75.45698 ymax: 36.58965
epsg (SRID):    4267
proj4string:    +proj=longlat +datum=NAD27 +no_defs

and some random points in the bounding box, given the same CRS:

> set.seed(123)
> pts = st_as_sf(data.frame(x=runif(50,-80.8,-80.0),y=runif(50,35.5, 36.37),N=paste0("A",1:50)),coords=1:2)
> st_crs(pts)=st_crs(nc)

Then use st_intersects:

> overs = st_intersects(pts, nc)
although coordinates are longitude/latitude, st_intersects assumes that they are planar
although coordinates are longitude/latitude, st_intersects assumes that they are planar
> overs[1:5]
[[1]]
integer(0)

[[2]]
[1] 2

[[3]]
[1] 1

[[4]]
[1] 2

[[5]]
[1] 2

This is a list of length "number of points". Since its possible for a point to overlap more than one polygon (if the polygons overlap) each element could be more than one item, but in this case because the polygons don't overlap each element is length one when it does overlap, and length zero when they don't.

You can then do things like lengths(overs)>0 (note, not length but lengths to get the length of each element in a list) to test which points overlapped:

> lengths(overs) > 0
 [1] FALSE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
[13]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
[25]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE
[37] FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE  TRUE FALSE
[49] FALSE  TRUE

and then work on the values in overs to lookup rows in your polygons for each point conditional on lengths(overs)>0.