Search code examples
rgisgeospatialrasterr-sf

Check if polygons intersect in R

Many users have asked how to check whether polygons intersect, however, current answers to those questions are not working for my case.

I have two objects. One is a "Large SpatialPolygons" named "farms". This object has several polygons (total of 2011), and each polygon indicates the limits of a different farm (see print screen).

My second object is a "Large SpatialPolygons DataFrame" named slope_RJ_100m. That object divides a large area into several square polygons with area of 10.000m^2 each (total of 310000 polygons).

For each of the squares (polygons) in "slope_RJ_100m", I would like to know whether they intersect any of the polygons in "farms". In other words, I want to know whether each particular square in "slope_RJ_100m" has a farm inside (even if just a piece of a farm). I was expecting the outcome to be something with 310000 rows and two variables, one indicating the polygon in slope_RJ_100m, and the other with TRUE or FALSE for whether that polygon has a farm.

I have tried:

inters = gIntersection(slope_RJ_100m, farms)

This one produces an output of about 1500 polygons. I am not sure how to use this to know which of my 310000 polygons has a farm in it.

inters = raster::intersect(slope_RJ_100m, farms)

The output has 29144 polygons. As in the previous case, not sure how I can use this to know whether the square has a farm.

and

inters = st_intersects(slope_RJ_100m, farms)

Error in UseMethod("st_intersects") : no applicable method for 'st_intersects' applied to an object of class "c('SpatialPolygonsDataFrame', 'SpatialPolygons', 'Spatial', 'SpatialPolygonsNULL', 'SpatialVector', 'SPNULL')"

I am sure my question is trivial and I apologize beforehand.

enter image description here

Solution

  • Here is an example with terra (to run this you need terra 1.1-17)

    library(terra)
# polygons
p1 <- vect("POLYGON ((0 0, 8 0, 8 9, 0 9, 0 0))")
p2 <- vect("POLYGON ((5 6, 15 6, 15 15, 5 15, 5 6))")
p3 <- vect("POLYGON ((8 2, 9 2, 9 3, 8 3, 8 2))")
p4 <- vect("POLYGON ((2 6, 3 6, 3 8, 2 8, 2 6))")
p5 <- vect("POLYGON ((2 12, 3 12, 3 13, 2 13, 2 12))")
p6 <- vect("POLYGON ((10 4, 12 4, 12 7, 11 7, 11 6, 10 6, 10 4))")

p <- rbind(p1, p2, p3, p4, p5, p6)
plot(p, col=rainbow(6, alpha=.5))
lines(p, lwd=2)
text(p)

    enter image description here

    relate(rbind(p1, p2), rbind(p3,p4,p5,p6), "intersects")
#      [,1]  [,2]  [,3]  [,4]
#[1,]  TRUE  TRUE FALSE FALSE
#[2,] FALSE FALSE FALSE  TRUE

    With your SpatialPolygons* you should be able to do the following:

    s <- vect(slope_RJ_100m)
f <- vect(farms)
x <- relate(s, f, "intersects")

    Likewise, to use st_intersects, you need to use sf objects, not Spatial* objects. Something like

    library(sf)
ss <- st_as_sf(slope_RJ_100m)
ff <- st_as_sf(farms)
inters <- st_intersects(ss, ff)

    With the example data from above

    s1 < st_as_sf(rbind(p1, p2))
s2 <- st_as_sf(rbind(p3,p4,p5,p6))
st_intersects(s1, s2)
#Sparse geometry binary predicate list of length 2, where the predicate was `intersects'
# 1: 1, 2
# 2: 4