在R中相交空间线

Question

我在 R 中有下一个空间对象。

library(sp)
library(rgeos)
poly1 <- structure(c(-3.25753225, -3.33532866, -3.33503723, -3.35083008, 
                      -3.35420388, -3.407372, -3.391667, -3.254167, -3.248129, -3.25753225, 
                      47.78513433, 47.73738617, 47.73793803, 47.74440261, 47.74004583, 
                      47.803846, 47.866667, 47.866667, 47.806292, 47.78513433),
                    .Dim = c(10L, 2L), .Dimnames = list(NULL, c("x", "y")))
poly2 <- structure(c(-3.101871, -3.097764, -3.20532, -3.260711, -3.248129, 
                      -3.101871, 47.777041, 47.735975, 47.709087, 47.777982, 47.806292, 47.777041),
                    .Dim = c(6L, 2L), .Dimnames = list(NULL, c("x", "y")))
sobj <- SpatialPolygons(
    list(
        Polygons(list(Polygon(poly1)), ID = '1'),
        Polygons(list(Polygon(poly2)), ID = '2')),
    proj4string = CRS('+proj=merc'))
plot(sobj)

我想获取一个空间对象，其中包含两个多边形共有的边界线，即下一张图像中的绿色线。

lines <- matrix(c(-3.248129, -3.25753225, 47.806292, 47.78513433), 2, 2)
lobj <- SpatialLines(
    list(
        Lines(list(Line(lines)), ID = '1')),
    proj4string = CRS('+proj=merc'))

plot(lobj, col = 'green', add = TRUE)

lines <- matrix(c(-3.248129, -3.25753225, 47.806292, 47.78513433), 2, 2)
lobj <- SpatialLines(
    list(
        Lines(list(Line(lines)), ID = '1')),
    proj4string = CRS('+proj=merc'))

plot(lobj, col = 'green', add = TRUE)

到目前为止，我已尝试使用 rgeos 包中的 gIntersection 函数，但它并没有达到我的要求。我怎么会得到这个？

Answer 1

如果您的线条完全重叠，我认为rgeos::gIntersection 将是首选方法。考虑以下简单示例：

l1 <- SpatialLines(list(Lines(list(Line(rbind(c(1, 1), c(5, 1)))), 1)))
l2 <- SpatialLines(list(Lines(list(Line(rbind(c(3, 1), c(10, 1)))), 1)))

plot(0, 0, ylim = c(0, 2), xlim = c(0, 10), type = "n")
lines(l1, lwd = 2, lty = 2)
lines(l2, lwd = 2, lty = 3)
lines(gIntersection(l1, l2), col = "red", lwd = 2)

您的问题的一个解决方案，虽然不是完美的，也许其他人有更好的解决方案，就是添加一个小缓冲区。

xx <- as(sobj, "SpatialLines")
xx <- gBuffer(xx, width = 1e-5, byid = TRUE)
xx <- gIntersection(xx[1, ], xx[2, ])

plot(sobj)
plot(xx, border = "red", add = TRUE, lwd = 2)