R SF Intersects – Using ST_Intersects to Classify Points with R and SF

Question

I have two objects, a shapefile of location points, points, and a shapefile of polygons, world_buffer (both +proj=utm +datum=WGS84). The world_buffer is the world's coastlines, found here, with a 5 km buffer around it, computing using gBuffer().

I am trying to classify the "points" as either within the world polygons (land), within the buffer (coast), or outside of the buffer (ocean).

I have tried using st_intersects() to get the points within the buffer:

library(sf)
intersect <- st_intersects(points, world_buffer, sparse = TRUE, prepared = FALSE)

But this returns a large sgbp that has two columns, a column with an intersect number and then a column full of integer [0].

Answer 1

Let's use some sample data from sf to do this.

library(sf)
# North Carolina - lets' pretend its an island...
nc = st_read(system.file("shape/nc.shp", package="sf"))

# Turn it into a rough metre-length system by projection
nc = st_transform(nc, 3857)

# and make a 5km buffer
ncbuff = st_buffer(nc, 5000)

# scatter 1000 points around.
pts = st_jitter(st_sample(nc, 1000), factor=0.2)

The intersection tests return lists of which polygons a point is in (because polygons can overlap, so a point could be in more than one polygon). If a point is in no polygons then you get a zero-length vector.

One way of testing if a point is in any polygon is to use lengths on the list, which returns a vector of how many elements is in each element. If this is zero then the point is not in any polygons. So we can do:

buff = lengths(st_intersects(pts, ncbuff)) > 0
land = lengths(st_intersects(pts, nc)) > 0

this gives us two TRUE/FALSE vectors telling us if the corresponding pts element is in each set of polygons. Logic then says:

coast = buff & !land
ocean = !buff

Now we can plot subsets of the points in different colours over the polygon geometry:

plot(nc$geom)
plot(pts[ocean], pch=19, col="blue", add=TRUE);
plot(pts[coast], pch=19, col="red", add=TRUE);
plot(pts[land], pch=19, col="green", add=TRUE);

To get: