Keep spatial lines intersecting spatial polygons using R?

Question

I am trying to keep spatial Lines objects intersecting a spatial polygon in R. The intersect function from raster package is doing actually a crop. But what I need is to keep the lines intersecting a polygon (but not to crop them). In ArcMap there was a possibility for that (something like select by attributes using a polygon).

Some code I created:

library(sp)
library(raster)

# from the sp vignette:
l1 = cbind(c(-35,-32, -31, -18, -15, 0, 22, 24 ), c(38, 39, 41, 44, 45, 48, 55, 66))
rownames(l1) = letters[1:8]
l1a = cbind(l1[,1]+25,l1[,2]+25)
rownames(l1a) = letters[1:8]
l2 <- l1-66
rownames(l2) = letters[1:8]
l3 <- l2 - 28.2
rownames(l3) = letters[1:8]


Sl1 = Line(l1)
Sl2 = Line(l2)
Sl3 = Line(l3)
S1 = Lines(list(Sl1), ID="a")
S2 = Lines(list(Sl2), ID="b")
S3  = Lines(list(Sl3), ID  = 'c')
Sl = SpatialLines(list(S1,S2, S3))


poly <- extent(c(-45, 15, -15, 57))
poly <- as(poly, 'SpatialPolygons')  


plot(Sl, col = c("red", "blue", 'green'))
plot(poly, add  = T, lwd = 2) # here I should keep the red and blue lines

intersected <- raster::intersect(Sl, poly)

plot(intersected)
plot(poly, add=T) # I managed only to crop the lines which is not exactly what I wanted.

Answer 1

You need an intersection test that tells you if two geometries intersect. You can use gIntersects from rgeos:

> library(rgeos)
rgeos version: 0.3-28, (SVN revision 572)
 GEOS runtime version: 3.6.2-CAPI-1.10.2 4d2925d6 
 Linking to sp version: 1.3-1 
 Polygon checking: TRUE 

Use byid=TRUE to test if each feature in Sl intersects poly:

> gIntersects(Sl, poly, byid=TRUE)
     a    b     c
1 TRUE TRUE FALSE

You can then subset to get the intersecting features:

> Sli = Sl[gIntersects(Sl, poly, byid=TRUE),]
> plot(Sli, add=TRUE, col="black", lwd=2)

This shows the features are unchanged, just selected.

You can also do this with over:

> over(Sl, poly)
 a.1  b.1 c.NA 
   1    1   NA 

which tells you which element of poly each feature in Sl is over, or NA if the feature is not over anything in poly. To select here, use is.na on the output from over.

If you start using the sf package instead, you want st_intersects.