R - incorrect results with raster::intersect for spatial lines with polygons

Question

I need to test (i.e., TRUE/FALSE), one pair at a time, whether SpatialLinesDataFrame elements intersect with SpatialPolygonsDataFrame elements. If TRUE, I then erase the portion of each polygon that intersects each line (one of each at a time) and save to a new list.

The polygons are stored in a multifeature SpatialPolygonsDataFrame and the lines are stored in a list of single SpatialLinesDataFrames. I have created a nested loop to iterate thru each of the polygons and line elements. I am using the raster::intersect, raster::extent and raster::erase functions.

library(sp)
library(raster)

#create  multi-feature SpatialPolygonDataFrame
p <- 
SpatialPolygons(list(Polygons(list(Polygon(cbind(c(2,4,3,2),c(2,2,4,2)))), 
"1"),

Polygons(list(Polygon(cbind(c(12,14,13,12),c(0,0,2,0)))), "2"),

Polygons(list(Polygon(cbind(c(25,24,22,25),c(22,23,22,22)))), "3"),

Polygons(list(Polygon(cbind(c(0,2,1,0),c(12,12,14,12)))), "4"),

Polygons(list(Polygon(cbind(c(5,7,6,5),c(5,5,7,5)))), "5"))) 
# Create a dataframe and display default rownames
p.df <- data.frame( ID=1:length(p)) 
rownames(p.df)
# Extract polygon ID's
pid <- sapply(slot(p, "polygons"), function(x) slot(x, "ID")) 
# Create dataframe with correct rownames
p.df <- data.frame( ID=1:length(p), row.names = pid)   
# coersion 
p <- SpatialPolygonsDataFrame(p, p.df)

#create list of single-feature SpatialLineDataFrame
l1 <- cbind(c(0,3), c(0,3))
l2 <- cbind(c(0, 13), c(0, 1))
l3 <- cbind(c(0, 24), c(0,22.5))
l4 <- cbind(c(0, 1), c(0,13))
l5 <- cbind(c(0, 6), c(0,6))

Sl1 <- Line(l1)
Sl2 <- Line(l2)
Sl3 <- Line(l3)
Sl4 <- Line(l4)
Sl5 <- Line(l5)

Sl1 <- Lines(list(Sl1), ID = "1")
Sl2 <- Lines(list(Sl2), ID = "2")
Sl3 <- Lines(list(Sl3), ID = "3")
Sl4 <- Lines(list(Sl4), ID = "4")
Sl5 <- Lines(list(Sl5), ID = "5")

Sl <- SpatialLines(list(Sl1, Sl2, Sl3, Sl4, Sl5))

# Create a dataframe and display default rownames
Sl.df <- data.frame( ID=1:length(Sl)) 
rownames(Sl.df)
# Extract polygon ID's
pidl <- sapply(slot(Sl, "lines"), function(x) slot(x, "ID")) 
# Create dataframe with correct rownames
Sl.df <- data.frame( ID=1:length(Sl), row.names = pidl)   
# coersion 
Sldf <- SpatialLinesDataFrame(Sl, Sl.df)

#convert multipart SpatialLineDataFrame feature to individual features in 
list
linel <- list()
linel[[1]] <- Sldf[1,]
linel[[2]] <- Sldf[2,]
linel[[3]] <- Sldf[3,]
linel[[4]] <- Sldf[4,]
linel[[5]] <- Sldf[5,]

#if a line intersects with a polygon, erase part of line where polygon 
#overlaps, save to new list, and add original line ID + ID of intersected 
#polygon
list1 <- list()
index=0
for (i in 1:length(linel)) {  
  for (j in 1:length(p)) {
    if (!is.null(raster::intersect(extent(p[j,]), linel[[i]]))) {
      index=index+1
      list1[[index]] <- erase(linel[[i]],p[j,])
      list1[[index]]@data["id.parent"] <- linel[[i]]@data$ID
      #add intersected polygon ID to the line attribute table 
      list1[[index]]@data["id.intersect"] <- p[j,]@data$ID
    }}}
list1

#check results by plotting
plot(p)
plot(list1[[1]], add=T) #parent 1-1
plot(list1[[2]], add=T) #parent 2-2
plot(list1[[3]], add=T) #parent-intersect 3-1
plot(list1[[4]], add=T) #??? p-i 3-2 BUT DOES NOT INTERSECT WITH p[2,]!
plot(list1[[5]], add=T) #parent 3-3
plot(list1[[6]], add=T) #???? p-i 3-4 BUT DOES NOT INTERSECT WITH p[4,]!
plot(list1[[7]], add=T) #p-i 3-5
plot(list1[[8]], add=T) #parent 4-4
plot(list1[[9]], add=T) #p-i 5-1
plot(list1[[10]], add=T) #parent 5-5

#test raster::intersect function on linel[[3]] based on plotted errors 
#above:
!is.null(raster::intersect(extent(p[1,]), linel[[3]])) 
#output = 'TRUE' (correct)
!is.null(raster::intersect(extent(p[2,]), linel[[3]])) 
#output = 'TRUE' (INCORRECT)
!is.null(raster::intersect(extent(p[3,]), linel[[3]])) 
#output = 'TRUE' (correct)
!is.null(raster::intersect(extent(p[4,]), linel[[3]])) 
#output = 'TRUE' (INCORRECT)
!is.null(raster::intersect(extent(p[5,]), linel[[3]])) 
#output = 'TRUE' (correct)

Most of the results are as expected, but several of the raster::intersect boolean results are incorrect; i.e., they say a polygon intersects a line when it does not. Why is this happening?

Robert Hijmans Robert Hijmans · Accepted Answer · 2019-04-05T08:14:06

Here is a simplified way to create these objects

library(raster)
p1 <- cbind(c(2,4,3,2),c(2,2,4,2))
p2 <- cbind(c(12,14,13,12),c(0,0,2,0))
p3 <- cbind(c(25,24,22,25),c(22,23,22,22))
p4 <- cbind(c(0,2,1,0),c(12,12,14,12))
p5 <- cbind(c(5,7,6,5),c(5,5,7,5)) 
p <- spPolygons(p1, p2, p3, p4, p5, attr=data.frame(polID=1:5))

#create list of single-feature SpatialLineDataFrame
l1 <- cbind(c(0,3), c(0,3))
l2 <- cbind(c(0, 13), c(0, 1))
l3 <- cbind(c(0, 24), c(0,22.5))
l4 <- cbind(c(0, 1), c(0,13))
l5 <- cbind(c(0, 6), c(0,6))
Sldf <- spLines(l1, l2, l3, l4, l5, attr=data.frame(lineID=1:5))

#linel <- lapply(1:5, function(i) Sldf[i,])

Why not do this in one step?

x <- intersect(Sldf, p)
plot(x)

data.frame(x)
#  lineID polID
#1      1     1
#2      2     2
#3      3     1
#4      3     3
#5      3     5
#6      4     4
#7      5     1
#8      5     5

You can get the same result like this:

result <- list()
index <- 1
for (i in 1:length(Sldf)) {  
   for (j in 1:length(p)) {
     if (is.null(raster::intersect(extent(p[j,]), extent(Sldf[i, ])))) next
     #if (is.null(raster::intersect(p[j,], Sldf[i, ]))) next
     if (!rgeos::gIntersects(p[j,], Sldf[i, ])) next
     e <- erase(Sldf[i,], p[j,])
     e$polID <-  p[j,]$polID
     result[[index]] <- e 
     index <- index + 1
}}

t(sapply(result, data.frame))
#     lineID polID
#[1,] 1      1    
#[2,] 2      2    
#[3,] 3      1    
#[4,] 3      3    
#[5,] 3      5    
#[6,] 4      4    
#[7,] 5      1    
#[8,] 5      5

Note that the difference with your code is that you do not check whether the polygon intersects the the line. You only check if it intersects with the extent of the line. That can be a quick check, but it is not sufficient for your purposes.

R - incorrect results with raster::intersect for spatial lines with polygons

1 Answers