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

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的修订方案非常好用。然而，如果一个多边形与其相邻的多边形没有共享的边缘（即它到达一个点，然后创建一个岛屿slither多边形），则在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")

- Chris

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)

- Sébastien Rochette

---

谢谢您。根据我的样本数据，它完全符合我的要求。然而，我已经将您的解决方案应用到了我的“真实”数据集中，并遇到了一些问题。我已经编辑了我的问题以展示我的意思。 - Chris

谢谢。你的代码完美地运行了。我现在唯一的问题是我的一些输入多边形被糟糕地绘制了。我无法解决这个问题，所以我需要找到一个解决方案，要么跳过有问题的多边形（即在res中返回NA），要么能够处理滑动岛屿多边形。对我来说，第一种选择似乎是更好的解决方案，我只是找不到一个if语句，在“rgeos :: gIntersection（Shapefile [n，]，Shapefile [Touching_List [[n]]，]，byid = TRUE）”之前，能够识别问题多边形。我已经添加了一组8个多边形到我的问题中，以演示我的意思。 - Chris

使用try或tryCatch，如果lines返回错误，则使用if语句的逐个多边形方法，并在其中再次使用try。 - Sébastien Rochette

Answer 2

只是补充一下Sébastien Rochette的答案，我认为sf包中的st_length函数不适用于多边形（参见这个post）。相反，我建议使用lwgeom包中的st_perimeter函数。

（我想评论这个答案，但我没有足够的声望）

- FAmorim