Мне нужно проверить (т.е. ИСТИНА / ЛОЖЬ) по одной паре за раз, пересекаются ли SpatialLinesDataFrame элементы с SpatialPolygonsDataFrame элементами. Если ИСТИНА, я стираю часть каждого многоугольника, пересекающего каждую линию (по одной каждой за раз), и сохраняю в новом списке.
Полигоны сохраняются в мультиэлементе 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 неверны; т.е. они говорят, что многоугольник пересекает линию, когда это не так. Почему это происходит?
