# R: Using st_intersects to classify points inside, outside and within a buffer

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]`.

• The easier you can make your question reproducible the better - there's lots of sample data included with R's spatial packages and show all your working up to the point you get stuck. The less we have to do to get to the point the more likely you are to get an answer! Commented Apr 26, 2021 at 16:44

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

``````library(sf)
# North Carolina - lets' pretend its an island...

# 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)