计算多个多边形之间共享边界的长度
Calculate the length of shared boundaries between multiple polygons
我有一个 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 多边形的总 length/perimeter 以及每个 TOUCHING 多边形接触原点多边形的总长度。这将允许计算共享边界的百分比。我可以想象它的输出将是 Touching_DF 中的 3 个新列（例如，对于第一行，它可能类似于原始参数 1000m，接触长度 500m，共享边界 50%）。谢谢。
编辑 1
我已将@StatnMap 的答案应用到我的真实数据集。如果多边形共享边和点，则 gTouches 似乎是 returning 结果。这些点引起问题是因为它们没有长度。我修改了 StatnMap 的代码部分来处理它，但是在最后创建数据框时，共享 edges/vertices gTouches returns 的数量与长度的边数量不匹配.
这里有一些代码使用我的实际数据集的样本来演示问题：
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 设置为 return 仅共享长度大于零或 2. return 长度为零（而不是错误）时遇到一个点而不是一条边。到目前为止，我找不到任何可以做这些事情的东西。
编辑 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
我一直在寻找一种解决方案，该解决方案能够识别边界绘制不当的多边形并且不执行任何计算并且 return 'NA' 在 res 中（因此它们仍然可以稍后识别）。但是，我一直无法找到将这些有问题的多边形与 'normal' 多边形区分开来的命令。
运行 @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")
两个只相互接触的多边形的交点是一条线。使用 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)
在上面的table中，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 的回答，我认为 sf 包中的函数 st_length 不适用于多边形（参见 post）。相反，我建议在 lwgeom 包中使用函数 st_perimeter。
（我想评论答案，但我没有足够的声誉）
我有一个 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 多边形的总 length/perimeter 以及每个 TOUCHING 多边形接触原点多边形的总长度。这将允许计算共享边界的百分比。我可以想象它的输出将是 Touching_DF 中的 3 个新列（例如，对于第一行，它可能类似于原始参数 1000m，接触长度 500m，共享边界 50%）。谢谢。
编辑 1
我已将@StatnMap 的答案应用到我的真实数据集。如果多边形共享边和点，则 gTouches 似乎是 returning 结果。这些点引起问题是因为它们没有长度。我修改了 StatnMap 的代码部分来处理它，但是在最后创建数据框时，共享 edges/vertices gTouches returns 的数量与长度的边数量不匹配.
这里有一些代码使用我的实际数据集的样本来演示问题：
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 设置为 return 仅共享长度大于零或 2. return 长度为零（而不是错误）时遇到一个点而不是一条边。到目前为止，我找不到任何可以做这些事情的东西。
编辑 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
我一直在寻找一种解决方案，该解决方案能够识别边界绘制不当的多边形并且不执行任何计算并且 return 'NA' 在 res 中（因此它们仍然可以稍后识别）。但是，我一直无法找到将这些有问题的多边形与 'normal' 多边形区分开来的命令。
运行 @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")
两个只相互接触的多边形的交点是一条线。使用 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)
在上面的table中，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 的回答，我认为 sf 包中的函数 st_length 不适用于多边形（参见 post）。相反，我建议在 lwgeom 包中使用函数 st_perimeter。
（我想评论答案，但我没有足够的声誉）