1

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

1 Answer 1

4

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)
2
  • the st_join function doesnt seem to work for my data - it just returns all my points with no distinction made anywhere in the dataframe between points that intersect the polygons and points that dont
    – user303287
    Commented Aug 8, 2022 at 17:36
  • If st_join returns a dataframe with no NA values in the columns from the polygons then none of the points are outside the polygons. Anyway, the thing to use for testing is st_intersects if you then want to select for in/out. As in the answer.
    – Spacedman
    Commented Aug 8, 2022 at 20:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.