在R中求空间线的交点

4

我有一个在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函数,但它没有达到我所需的效果。如何实现?

输入图像描述

1个回答

4

如果您的线完全重叠,我认为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)

enter image description here

针对您的问题,虽然不是完美的解决方法,也许还有更好的解决方案,但是可以考虑添加一个小缓冲区。

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)

enter image description here


那真是聪明。我想我可以通过使非常接近的点成为唯一的点,将生成的多边形转换为线。 - Usobi
ж€–иЂ…дҢ еЏҮд»ӨдҢүз”ЁgIntersectionе‡Ңж•°дёҺеҺџе§‹зғүжқҰиү›иҰЊдғ¤й›†иүђз®—пәЊдң‹е¦‚gIntersection(as(sobj[1, ], "SpatialLines"), xx)пәЊиү™е°†д»…з»™е‡ғж‚ЁеЏҮиѓҢйњЂи¦Ѓзљ„иң№з•ЊйѓЁе€†гЂ‚ - johannes

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接