1

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
2

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.

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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