R - 带有多边形的空间线与raster :: intersect的结果不正确

Question

我需要测试（即，TRUE / FALSE），一次一对，SpatialLinesDataFrame元素是否与SpatialPolygonsDataFrame元素相交。如果为TRUE，则我擦除每个与每行相交的多边形部分（每次一个），并保存到新列表。

多边形存储在多功能SpatialPolygonsDataFrame中，并且这些线存储在单个SpatialLinesDataFrames列表中。我创建了一个嵌套循环来迭代每个多边形和线元素。我正在使用raster::intersect，raster::extent和raster::erase函数。

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)

大多数结果都是预期的，但是几个raster::intersect布尔结果是不正确的;也就是说，他们说多边形与一条线相交。为什么会这样？

Answer 1

这是创建这些对象的简化方法

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

为什么不一步到位呢？

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

您可以得到相同的结果：

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    

请注意，与代码的不同之处在于，您不检查多边形是否与该线相交。您只检查它是否与线条的范围相交。这可以是一个快速检查，但它不足以满足您的目的。