Suite

Découper un SpatialPoygonsDataFrame à un autre SpatialPolygonsDataFrame

Découper un SpatialPoygonsDataFrame à un autre SpatialPolygonsDataFrame


Comment restreindre unSpatialPolygonsDataFramede sorte qu'aucun des polygones ne se trouve en dehors d'un "sous-jacent"SpatialPolygonsDataFrame?

L'application que j'ai en tête implique Voronoi Polygons à la Carson Farmer -- le résultat d'une telle décomposition est (généralement) rectangulaire par défaut, mais il est souvent plus attrayant visuellement de conformer le résultat à une géométrie sous-jacente (ville/comté/frontière du pays, etc.). Une autre application courante consiste à supprimer l'Alaska et Hawaï des géométries à l'échelle des États-Unis lorsque seuls les 48 inférieurs présentent un intérêt.

C'est similaire à cette question, mais plus compliqué, puisque les deux géométries d'entrée se composent de plusieurs polygones - un simplegIntersectionentraînera une dissolution inutile de la géométrie principale d'intérêt.

C'est aussi similaire à cette question mais je n'ai pas pu suivre exactement ce qui se passait là-bas car il n'y avait pas d'exemple reproductible. Je pense que l'exercice suivant subsume probablement cette question.

Voici un exemple simple avec lequel travailler :

bibliothèque(sp) bibliothèque(rgeos) bibliothèque(maptools) points<-list("A" = c(0, 0),"B" = c(1/2, 0),"C" = c(1, 0 ), "D" = c(0, 1/2), "E" = c(1/2, 1/2), "F" = c(1, 1/2), "G" = c(0 , 1), "H" = c(1/2, 1), "I" = c(1, 1), "J" = c(0, -1/2), "K" = c(3/ 2, -1/2), "L" = c(7/24, 3/8), "M" = c(5/8, 3/8), "N" = c(3/8, 5/ 8), "O" = c(5/8, 17/24), "P" = c(3/2, 1)) pts_poly <- function(x)Reduce(rbind, lapply(x, function(y) unlist(points[y]))) poly1_list<- list("poly11" = list(ID = "A", coords = pts_poly(c("A", "B", "E", "D" ))) , "poly12" = liste(ID = "B", coords = pts_poly(c("B", "C", "F", "E"))), "poly13" = liste(ID = "C", coords = pts_poly(c("E", "F", "I", "H"))), "poly14" = list(ID = "D", coords = pts_poly(c("D", "E" , "H", "G")))) poly2_list<- list("poly11" = list(ID = "a", coords = pts_poly(c("J", "M", "L"))), "poly12" = liste(ID = "b", coords = pts_poly(c("J", "K", "M"))), "poly13" = liste(ID = "c", coords = pts_poly(c ("M", "K", "P"))), "poly14" = list(ID = "d", coords = pts_poly(c("M", "P", "O"))), " poly15" = liste (ID = "e", coords = pts_poly(c("L", "M", "O", "N")))) poly1<-SpatialPolygons(lapply( poly1_list, function(pl){ with(pl, Polygones(list(Polygon(coords)), ID = ID))})) df1 <- data.frame(var1 = rnorm(length(poly1)), var2 = rnorm(length(poly1)), row.names = sapply ([email protected], fonction(x) [email protected])) polydf1 <- SpatialPolygonsDataFrame(poly1, df1) poly2<-SpatialPolygons(lapply( poly2_list, function(pl){ with(pl, Polygons(list(Polygon(coords)) , ID = ID))})) df2 <- data.frame(var1 = rnorm(length(poly2)), var2 = rnorm(length(poly2)), row.names = sapply([email protected], function(x) [email protected])) polydf2 <- SpatialPolygonsDataFrame(poly2, df2)

Beaucoup de code. Voici ce que nous avons visuellement :

(code pour les parcelles :)

par(mfrow = c(1, 3)) par(mar=c(0,0,0,0)) par(oma=c(0,0,0,0)) plot(polydf1, xlim = c(0 , 1.1), ylim = c(0, 1.1)) title("Polygones sous-jacents", line = -15, cex.main=1.8) text(Reduce(rbind, points[1:9]) + .03, labels = noms(points)[1:9], col = "rouge") text(matrice(c(.25, .25, .75, .25, .75, .75, .25, .75), ncol = 2 , byrow = TRUE), cex = 2, labels=sapply([email protected], slot, name = "ID")) plot(polydf2, xlim = c(0, 1.6), ylim = c(-.5, 1.1) ) title("Polygones à découper", line = -15, cex.main=1.8) text(Reduce(rbind, points[-(1:9)]) + .03, labels = names(points)[-( 1:9)], col = "red") text(coordinates(polydf2), cex = 2, labels=sapply([email protected], slot, name = "ID")) plot(polydf2) plot(polydf1, add = TRUE) title("Superposition de polygones", ligne = -15, cex.main=1.8)

Ce que je veux, c'est éliminer les parties des triangles E, F, G et H qui se trouvent à l'extérieur du premier ensemble de polygones, A, B, C et D.

Une simple intersection est incorrecte :

plot(gIntersection(polydf1, polydf2, byid = TRUE), main = "Simple Intersection : Overkill")

Voici à quoi ressemble ma sortie souhaitée (en éliminant essentiellement les bordures de A, B, C et D):


L'approche la plus simple pour ce faire est

bibliothèque(raster) x <- crop(polydf2, polydf1)

Cela fera l'intersection géométrique, mais assurera également que les attributs sont en ordre. Voir?'paquet-raster'(section XIV) pour plus de fonctions qui traitent de la superposition de polygones.


La réponse est assez simple une fois que vous la voyez, mais j'ai passé une demi-journée à me cogner la tête sur le bureau avant de réaliser le simple ajustement que nous devons faire ; Je ne pouvais le trouver nulle part ailleurs en ligne, alors je voulais laisser ceci ici pour la référence de toute autre personne en difficulté comme moi.

La clé est d'éliminer les frontières du premier ensemble de polygones avant de croiser les polygones importants :

clipped_polys <- gIntersection(polydf2, gUnaryUnion(polydf1), byid = TRUE, id = sapply([email protected], slot, name = "ID"))

Pour nous assurer que les polygones conservent leur nom d'origine, nous spécifions également l'emplacement d'identification nous-mêmes. Voici la sortie :

Et bien sûr, nous pouvons rajouter toutes les données associées àpoly2facilement comme indiqué dans cette réponse:

clipped_polys <- SpatialPolygonsDataFrame(clipped_polys, [email protected])

Voir la vidéo: ArcGIS zonal statistics: distribution of meanaverage spatial values to polygonu0026sub polygons