# Keep spatial lines intersecting spatial polygons using R?

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.
``````
• You're missing `Sl1a` in your example. I guess its `Sl1a = Line(l1a)`? – Spacedman Jul 23 at 11:13
• Hi @Spacedman, thanks for the heads up. I corrected it now. – Andrei Niță Jul 23 at 11:16

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
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),]
``````

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

• I didn't know that. Thanks a lot @Spacedman ! Always helpful – Andrei Niță Jul 23 at 11:21