我有一个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_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 <- [email protected]}
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 <- [email protected]}
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.返回长度为零(而不是错误),当遇到点而不是边时。到目前为止,我找不到任何可以做这些事情的事情。
EDIT 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”(因此它们仍然可以在稍后识别)的解决方案。但是,我一直无法找到一个区分这些有问题的多边形与“正常”多边形的命令。
运行@ 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")
感谢您的支持。根据我的样本数据,它完全符合我的要求。但是,我已将您的解决方案应用于“真实”数据集并遇到了一些问题。我编辑了我的问题来展示我的意思。 – Chris
我编辑了我的答案来解释这种情况。 –
谢谢。你的代码完美地工作。我现在唯一的问题是一些输入多边形被严重绘制。我无法对此做任何事情,所以我需要找到一个解决方案,它可以跳过不规则的多边形(即,在资源库中返回NA),或者可以处理滑动的岛多边形。第一个选项似乎是一个更好的解决方案,我只是在'rgeos :: gIntersection(Shapefile [n,],Shapefile [Touching_List [[n]],],byid = TRUE)'之前找不到if语句'能够识别问题多边形。我在我的问题中添加了一小组8个多边形来展示我的意思。 – Chris