Calculer la longueur des limites partagées entre plusieurs polygones

Question

J'ai un fichier de formes et je veux savoir pour chaque polygone ce que d'autres polygones le touchent. À cette fin, j'ai ce code:

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

Je veux maintenant aller plus loin et comprendre la mesure dans laquelle chaque touche de polygone d'autres polygones. Ce que je suis après pour chaque ligne dans Touching_DF est une longueur totale/périmètre pour chaque polygone ORIGIN et la longueur totale que chaque polygone TOUCHING est en contact avec le polygone d'origine. Cela permettra ensuite de calculer le pourcentage de la limite partagée. Je peux imaginer la sortie de ce serait 3 nouvelles colonnes dans Touching_DF (par exemple pour la première rangée, il pourrait être quelque chose comme le paramètre d'origine 1000m, toucher la longueur 500m, frontière partagée 50%). Merci.

EDIT 1

J'ai appliqué @ réponse de StatnMap à mon vrai jeu de données. Il semble que gTouches renvoie des résultats si un polygone partage à la fois une arête et un point. Ces points causent des problèmes parce qu'ils n'ont pas de longueur. J'ai modifié la portion de code de StatnMap pour y faire face, mais quand il s'agit de créer la trame de données à la fin, il y a un décalage entre le nombre d'arêtes/sommets partagés renvoyés par gTouches et le nombre d'arêtes ayant des longueurs.

Voici un code pour démontrer le problème en utilisant un échantillon de mon jeu de données réelles:

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

Cela montre spécifiquement la question pour l'un des polygones:

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) 

Les deux solutions possibles I voir sont 1. obtenir gTouches pour retourner seulement les bords partagés avec une longueur supérieure à zéro ou 2. renvoyant une longueur de zéro (plutôt que d'erreur) quand un point plutôt qu'un bord est rencontré. Jusqu'à présent, je ne trouve rien qui puisse faire l'une ou l'autre de ces choses.

EDIT 2

@ solution révisée de StatnMap fonctionne très bien. Toutefois, si un polygone ne partage pas un boarder Snapped avec son polygone (il va à un point et crée alors un polygone de Slither de l'île), il arrive avec cette erreur après 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 

Je suis voisin à la recherche pour une solution qui est capable d'identifier les polygones avec des bordures mal dessinées et de ne pas effectuer de calculs et de retourner 'NA' dans res (afin qu'ils puissent être identifiés plus tard). Cependant, j'ai été incapable de trouver une commande qui distingue ces polygones problématiques des polygones 'normaux'.

Exécution @ solution révisée de StatnMap avec ces 8 polygones démontre la question:

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

L'intersection de deux polygones ne se toucher est une ligne. Le calcul d'une longueur de ligne est facile avec les fonctions des bibliothèques spatiales dans R.
Lorsque vous avez commencé votre exemple avec la bibliothèque sp, vous trouverez une proposition avec cette bibliothèque. Cependant, je vous donne également une proposition avec la nouvelle bibliothèque sf.

polygones Calculer les limites partagées longueurs avec bibliothèque 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) 

Dans le tableau ci-dessus, t.length est la longueur touchante et t.pc est le pourcentage touchant en ce qui concerne au périmètre du polygone d'origine.

Edit: Certaines limites partagées sont des points (avec sp)

Comme commenté, certaines frontières peuvent être un point unique, au lieu de lignes. Pour rendre compte de ce cas, je suggère de doubler les coordonnées du point pour créer une ligne de longueur = 0. Cela nécessite de calculer les intersections avec les autres polygones un par un, lorsque ce cas apparaît.
Pour un seul polygone, nous pouvons tester ceci:

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

Pour tous dans une boucle 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) 

Cela peut également être appliquée avec la bibliothèque sf, mais comme vous a apparemment choisi de travailler avec sp, je ne vais pas mettre à jour le code pour cette partie. Peut-être plus tard ...

---- Fin ---- Edit

polygones Calculer les limites partagées longueurs avec bibliothèque sf

Les chiffres et les sorties sont les mêmes.

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)