計算多個多邊形之間共享邊界的長度

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_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.返回長度爲零（而不是錯誤），當遇到點而不是邊時。到目前爲止，我找不到任何可以做這些事情的事情。

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") 

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)