Retain spatial points outside of multiple polygons

Question

We can use st_intesection to retain points that intersect with one or more polygons.

tmp1 = st_intersection(st_as_sf(track_sp), st_as_sf(cols.buf_sp))
pts1 = as(tmp1, "Spatial")
plot(pts1)

What is the function that does the opposite? I.e. retain only points that DO NOT intersect with the polygons?

I tried st_disjoint, but this does not return an sf dataframe.

tmp2 = st_disjoint(st_as_sf(track_sp), st_as_sf(cols.buf_sp), sparse = F)

I feel like this should be something so straightforward, yet I cannot find a solution.

I have seabird GPS spatial points (> 250,000 points) and want to remove points that lie within a certain distance from the colony (3 buffer polygons).

I have seen this answer here, but do not quite follow: Opposite of ST_intersection

Answer 1

st_intersection does a sort of database join. The "point" you get back isn't the "point" you put in. Its the intersection of the polygon geometry and the point geometry. If you do st_intersection of two overlapping polygon geometries you only get the geometry of the overlap back. Same thing here.

st_intersection is like a database join where non-joining rows aren't included - a left join. You can do a full spatial join with st_join, which by default does a full join. You get back NA values for rows that had no intersection with polygons, for example:

   Z  FIPS FIPSNO CRESS_ID BIR74                   geometry
1  A  <NA>     NA       NA    NA POINT (-82.26209 34.38398)
2  B  <NA>     NA       NA    NA POINT (-81.56811 34.75251)
3  C 37059  37059       30  1207  POINT (-80.45774 36.0282)
4  D 37033  37033       17  1035 POINT (-79.41677 36.34003)
5  E 37001  37001        1  4672 POINT (-79.45147 35.94315)
6  F 37099  37099       50  1143 POINT (-83.07752 35.22026)
7  G  <NA>     NA       NA    NA POINT (-82.52233 36.73691)
8  H  <NA>     NA       NA    NA POINT (-80.14545 37.37475)
9  I 37117  37117       59  1549 POINT (-77.16134 35.84393)
10 J  <NA>     NA       NA    NA POINT (-76.77965 34.17137)

So you could filter out the ins and outs from that lot. BUT...

What you probably really want is the predicate functions for geometry - these are "verbs" rather than "nouns". So st_intersects returns if one set of geometries intersects another set, as opposed to st_intersection which is a noun.

st_intersects(pts, pols) returns a list because each point could be intersecting with multiple polygons (if the polygons overlap, for example). List elements of zero length are what you want:

st_intersects(pts, nc) Sparse geometry binary predicate list of length 10, where the predicate was `intersects'

 1: (empty)
 2: (empty)
 3: 40
 4: 11
 5: 27
 6: 66
 7: (empty)
 8: (empty)
 9: 36
 10: (empty)

which you can get with:

> lengths(st_intersects(pts, nc)) == 0
 [1]  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE FALSE  TRUE

So the points not in any polygons are:

> pts[lengths(st_intersects(pts, nc)) == 0,]
Simple feature collection with 5 features and 1 field
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -82.52233 ymin: 34.17137 xmax: -76.77965 ymax: 37.37475
Geodetic CRS:  NAD27
   Z                   geometry
1  A POINT (-82.26209 34.38398)
2  B POINT (-81.56811 34.75251)
7  G POINT (-82.52233 36.73691)
8  H POINT (-80.14545 37.37475)
10 J POINT (-76.77965 34.17137)