複数のポリゴン間の共有境界の長さを計算する

Question

シェイプファイルがあり、他のポリゴンにどのポリゴンが接触するかを各ポリゴンについて知りたい。その目的のために私はこのコードを持っています：

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％）。ありがとう。私は私の本当のデータセットにStatnMapの答え@適用されている1

EDIT。ポリゴンがエッジとポイントの両方を共有する場合、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) 

2つの解決策I 1. gTouchesに0または2より大きい長さの共有エッジのみを返すようにする。エッジではなくポイントに遭遇したときにエラーの代わりにゼロの長さを返す。これまでのところ、私はこれらのことのいずれかを行うことはできません。 StatnMapの改訂ソリューション@

EDIT 2

は素晴らしい作品。しかし、ポリゴンはその隣接ポリゴンとスナップ国境を共有していない場合（すなわち、それがポイントになり、その後、島のスリザーリンクのポリゴンを作成します）、それは私が探している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") 

Answer 1

2つのポリゴンの交差が自分自身だけに触れることはラインです。 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有する）点

あり、いくつかの国境ではなく線の一意点であってもよいです。このケースを説明するために、ポイントの座標を倍にしてlength = 0の行を作成することをお勧めします。この場合、他のポリゴンとの交点を1つずつ計算する必要があります。これはまた、ライブラリsfに適用されるかもしれないが、あなたは明らかで動作することを選んだよう

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

：lapplyループ内のすべてのために

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

：単一の多角形の場合は、我々はこれをテストすることができます
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)