计算多个多边形之间共享边界的长度

Question

我有一个shapefile，我想知道每个多边形还有哪些其他多边形接触它。为此我有这样的代码：

require("rgdal") 
require("rgeos") 

download.file("https://www.dropbox.com/s/vbxx9dic34qwz63/Polygons.zip?dl=1", "Polygons.zip") 

Shapefile <- readOGR(".","Polygons") 

Touching_List <- gTouches(Shapefile, byid = TRUE, returnDense=FALSE) 

Touching_DF <- setNames(stack(lapply(Touching_List, as.character)), c("TOUCHING", "ORIGIN")) 

我现在要更进一步，了解在何种程度上每个多边形触摸其他多边形。我在Touching_DF的每一行之后的每个ORIGIN多边形的总长度/周长以及每个多边形都接触原始多边形的总长度。这将允许计算共享边界的百分比。我可以想象这个输出将是Touching_DF中的3个新列（例如，对于第一行，它可以是起始参数1000m，触摸长度500m，共享边界50％）。谢谢。

编辑1

我已经申请@ StatnMap的回答让我真正的数据集。看起来gTouches返回的结果是多边形共享一个边和一个点。这些问题由于没有长度而引起问题。我已经修改了StatnMap的部分代码来处理它，但是当涉及到在最后创建数据框时，gTouches返回多少个共享边/顶点以及多少个边具有不同长度。

下面是一些代码用我的实际数据集的样本来说明这个问题：

library(rgdal) 
library(rgeos) 
library(sp) 
library(raster) 

download.file("https://www.dropbox.com/s/hsnrdfthut6klqn/Sample.zip?dl=1", "Sample.zip") 

unzip("Sample.zip") 
Shapefile <- readOGR(".","Sample") 

Touching_List <- gTouches(Shapefile, byid = TRUE, returnDense=FALSE) 

# ---- Calculate perimeters of all polygons ---- 
perimeters <- sp::SpatialLinesLengths(as(Shapefile, "SpatialLines")) 

# ---- All in a lapply loop ---- 
all.length.list <- lapply(1:length(Touching_List), function(from) { 
    lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
    if(class(lines) != "SpatialLines"){lines <- lines@lineobj} 
    l_lines <- sp::SpatialLinesLengths(lines, longlat=FALSE) 
    results <- data.frame(origin = from, 
        perimeter = perimeters[from], 
        touching = Touching_List[[from]], 
        t.length = l_lines, 
        t.pc = 100*l_lines/perimeters[from]) 
    results 
}) 

这尤其显示了多边形的一个问题：

from <- 4 
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
if(class(lines) != "SpatialLines"){lines <- lines@lineobj} 
l_lines <- sp::SpatialLinesLengths(lines, longlat=FALSE) 

plot(Shapefile[c(from, Touching_List[[from]]),]) 
plot(lines, add = TRUE, col = 1 + 1:length(Touching_List[[from]]), lwd = 2) 

两个可能的解决方案我请参阅1.获取gTouches仅返回长度大于零的共享边，或2.返回长度为零（而不是错误），当遇到点而不是边时。到目前为止，我找不到任何可以做这些事情的事情。

EDIT 2

@ StatnMap修订后的解决方案的伟大工程。但是，如果多边形不与它相邻的多边形共享抓拍边界（即它关系到一个点，然后创建一个岛屿滑行多边形），则与此错误后lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE)

Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false, : 
     Geometry collections may not contain other geometry collections 

我一直在找上来对于能够识别具有严重绘制边界的多边形并且不执行任何计算并在res中返回“NA”（因此它们仍然可以在稍后识别）的解决方案。但是，我一直无法找到一个区分这些有问题的多边形与“正常”多边形的命令。

运行@ StatnMap与这些多边形8修订后的解决方案演示了这个问题：

download.file("https://www.dropbox.com/s/ttg2mi2nq1gbbrq/Bad_Polygon.zip?dl=1", "Bad_Polygon.zip") 
unzip("Bad_Polygon.zip") 
Shapefile <- readOGR(".","Bad_Polygon") 

Answer 1

两个多边形的交集仅触摸自己是一条线。用R中的空间库函数计算线长很容易。
当您开始使用库sp的示例时，您会发现该库有一个命题。不过，我也给你一个新图书馆sf的建议。

共享边界计算多边形与库长度sp

require("rgdal") 
require("rgeos") 
library(sp) 
library(raster) 

download.file("https://www.dropbox.com/s/vbxx9dic34qwz63/Polygons.zip?dl=1", "Polygons.zip") 

unzip("Polygons.zip") 
Shapefile <- readOGR(".","Polygons") 

Touching_List <- gTouches(Shapefile, byid = TRUE, returnDense=FALSE) 

# Touching_DF <- setNames(utils::stack(lapply(Touching_List, as.character)), c("TOUCHING", "ORIGIN")) 

# ---- Calculate perimeters of all polygons ---- 
perimeters <- sp::SpatialLinesLengths(as(Shapefile, "SpatialLines")) 

# ---- Example with the first object of the list and first neighbor ---- 
from <- 1 
to <- 1 
line <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]][to],]) 
l_line <- sp::SpatialLinesLengths(line) 

plot(Shapefile[c(from, Touching_List[[from]][to]),]) 
plot(line, add = TRUE, col = "red", lwd = 2) 

# ---- Example with the first object of the list and all neighbors ---- 
from <- 1 
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
l_lines <- sp::SpatialLinesLengths(lines) 

plot(Shapefile[c(from, Touching_List[[from]]),]) 
plot(lines, add = TRUE, col = 1 + 1:length(Touching_List[[from]]), lwd = 2) 

# ---- All in a lapply loop ---- 
all.length.list <- lapply(1:length(Touching_List), function(from) { 
    lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
    l_lines <- sp::SpatialLinesLengths(lines) 
    res <- data.frame(origin = from, 
        perimeter = perimeters[from], 
        touching = Touching_List[[from]], 
        t.length = l_lines, 
        t.pc = 100*l_lines/perimeters[from]) 
    res 
}) 

# ---- Retrieve as a dataframe ---- 
all.length.df <- do.call("rbind", all.length.list) 

在上面的表中，t.length是触摸长度及t.pc是与问候触摸百分比到原点多边形的周长。

编辑：一些共享边界是点（与sp）

作为评论的，一些边界可以是唯一点而不是线。为了说明这种情况，我建议将点的坐标加倍以创建一条长度为0的线。当出现这种情况时，这需要逐个计算与其他多边形的交点。
对于一个多边形，我们可以测试这个：

# Example with the first object of the list and all neighbours 
from <- 4 
lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
# If lines and points, need to do it one by one to find the point 
if (class(lines) == "SpatialCollections") { 
    list.Lines <- lapply(1:length(Touching_List[[from]]), function(to) { 
    line.single <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]][to],]) 
    if (class(line.single) == "SpatialPoints") { 
     # Double the point to create a line 
     L1 <- rbind(line.single@coords, line.single@coords) 
     rownames(L1) <- letters[1:2] 
     Sl1 <- Line(L1) 
     Lines.single <- Lines(list(Sl1), ID = as.character(to)) 
    } else if (class(line.single) == "SpatialLines") { 
     Lines.single <- line.single@lines[[1]] 
     Lines.single@ID <- as.character(to) 
    } 
    Lines.single 
    }) 
    lines <- SpatialLines(list.Lines) 
} 

l_lines <- sp::SpatialLinesLengths(lines) 

plot(Shapefile[c(from, Touching_List[[from]]),]) 
plot(lines, add = TRUE, col = 1 + 1:length(Touching_List[[from]]), lwd = 2) 

对于所有在lapply循环：

# Corrected for point outputs: All in a lapply loop 
all.length.list <- lapply(1:length(Touching_List), function(from) { 
    lines <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]],], byid = TRUE) 
    if (class(lines) == "SpatialCollections") { 
    list.Lines <- lapply(1:length(Touching_List[[from]]), function(to) { 
     line.single <- rgeos::gIntersection(Shapefile[from,], Shapefile[Touching_List[[from]][to],]) 
     if (class(line.single) == "SpatialPoints") { 
     # Double the point to create a line 
     L1 <- rbind(line.single@coords, line.single@coords) 
     rownames(L1) <- letters[1:2] 
     Sl1 <- Line(L1) 
     Lines.single <- Lines(list(Sl1), ID = as.character(to)) 
     } else if (class(line.single) == "SpatialLines") { 
     Lines.single <- line.single@lines[[1]] 
     Lines.single@ID <- as.character(to) 
     } 
     Lines.single 
    }) 
    lines <- SpatialLines(list.Lines) 
    } 
    l_lines <- sp::SpatialLinesLengths(lines) 
    res <- data.frame(origin = from, 
        perimeter = perimeters[from], 
        touching = Touching_List[[from]], 
        t.length = l_lines, 
        t.pc = 100*l_lines/perimeters[from]) 
    res 
}) 

all.length.df <- do.call("rbind", all.length.list) 

这也与库sf被应用，但你显然选择了与合作sp，我不会更新这部分的代码。也许以后...

----编辑结束----

共享边界计算多边形与库长度sf

图和输出是相同的。

library(sf) 
Shapefile.sf <- st_read(".","Polygons") 

# ---- Touching list ---- 
Touching_List <- st_touches(Shapefile.sf) 
# ---- Polygons perimeters ---- 
perimeters <- st_length(Shapefile.sf) 

# ---- Example with the first object of the list and first neighbour ---- 
from <- 1 
to <- 1 
line <- st_intersection(Shapefile.sf[from,], Shapefile.sf[Touching_List[[from]][to],]) 
l_line <- st_length(line) 

plot(Shapefile.sf[c(from, Touching_List[[from]][to]),]) 
plot(line, add = TRUE, col = "red", lwd = 2) 

# ---- Example with the first object of the list and all neighbours ---- 
from <- 1 
lines <- st_intersection(Shapefile.sf[from,], Shapefile.sf[Touching_List[[from]],]) 
lines <- st_cast(lines) # In case of multiple geometries (ex. from=71) 
l_lines <- st_length(lines) 

plot(Shapefile.sf[c(from, Touching_List[[from]]),]) 
plot(lines, add = TRUE, col = 1:length(Touching_List[[from]]), lwd = 2) 

# ---- All in a lapply loop ---- 
all.length.list <- lapply(1:length(Touching_List), function(from) { 
    lines <- st_intersection(Shapefile.sf[from,], Shapefile.sf[Touching_List[[from]],]) 
    lines <- st_cast(lines) # In case of multiple geometries 
    l_lines <- st_length(lines) 
    res <- data.frame(origin = from, 
        perimeter = as.vector(perimeters[from]), 
        touching = Touching_List[[from]], 
        t.length = as.vector(l_lines), 
        t.pc = as.vector(100*l_lines/perimeters[from])) 
    res 
}) 

# ---- Retrieve as dataframe ---- 
all.length.df <- do.call("rbind", all.length.list)