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

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多边形接触原点多边形的总长度.然后,这将允许计算共享边界的百分比.我可以想象这将是3个新列的输出Touching_DF(例如,对于第一行,它可能是原点参数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

我一直在寻找能够识别具有严重绘制边界的多边形并且不执行任何计算并返回"NA"的解决方案res(因此它们仍可以在以后识别).但是,我无法找到区分这些有问题的多边形与"普通"多边形的命令.

使用这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

仅触及自身的两个多边形的交点是一条线.使用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)