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

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 多边形与原点多边形接触的总长度。这将允许计算共享边界的百分比。我可以想象它的输出将是 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. 当遇到点而不是边时返回长度为零(而不是错误)。到目前为止，我找不到任何可以做这些事情的东西。

编辑 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”(以便以后仍可以识别它们)。但是，我一直无法找到将这些有问题的多边形与“正常”多边形区分开来的命令。

使用这 8 个多边形运行 @StatnMap 的修订解决方案演示了这个问题:

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

只是为了补充 Sébastien Rochette 的答案，我认为 st_length 包中的函数 sf 不适用于多边形(请参阅此 post )。相反，我建议在 st_perimeter 包中使用函数 lwgeom。

(我想评论答案，但我没有足够的声誉)