我有一个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")
try
或tryCatch
,如果lines
返回错误,则使用if
语句的逐个多边形方法,并在其中再次使用try
。 - Sébastien Rochette