Couper le polygone et conserver les données?

13

J'ai ces deux polygones:

library(sp); library(rgeos); library(maptools)

coords1 <- matrix(c(-1.841960, -1.823464, -1.838623, -1.841960, 55.663696,
                    55.659178, 55.650841, 55.663696), ncol=2)
coords2 <- matrix(c(-1.822606, -1.816790, -1.832712, -1.822606, 55.657887,
                    55.646806, 55.650679, 55.657887), ncol=2)
p1 <- Polygon(coords1)
p2 <- Polygon(coords2)
p1 <- Polygons(list(p1), ID = "p1")
p2 <- Polygons(list(p2), ID = "p2")
myPolys <- SpatialPolygons(list(p1, p2))
spdf1 = SpatialPolygonsDataFrame(myPolys, data.frame(variable1 = c(232,
                                                                   242), variable2 = c(235, 464), row.names = c("p1", "p2")))
proj4string(spdf1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf1, col="red")

coords1a <- matrix(c(-1.830219, -1.833753, -1.821154, -1.830219, 55.647353,
                     55.656629, 55.652122, 55.647353), ncol=2)
p1a <- Polygon(coords1a)
p1a <- Polygons(list(p1a), ID = "p1a")
myPolys1 <- SpatialPolygons(list(p1a))
spdf2 = SpatialPolygonsDataFrame(myPolys1, data.frame(variable1 = c(2),
                                                      variable2 = c(3), row.names = c("p1a")))
proj4string(spdf2) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84
+towgs84=0,0,0")
plot(spdf2, col="yellow", add=T)

entrez la description de l'image ici

Je veux en couper des parties spdf1qui sont recoupées par spdf2. Cependant, je veux spdf1rester en tant que SpatialPolygonsDataFrame et conserver toute information contenue dans spdf1@data.

J'ai essayé gDifference comme suit, qui supprime certaines parties spdf1qui sont intersectées par spdf2, puis se convertit spdf1en SpatialPolygons (c'est-à-dire en supprimant les informations contenues dans spdf1@data).

gDifference(spdf1, spdf2, byid=T)

Comment puis-je découper spdf1avec spdf2mais conserver les données contenues dans spdf1@data? J'ai vérifié et essayé cette question similaire sans savoir comment superposer un polygone sur SpatialPointsDataFrame et conserver les données SPDF?

luciano
la source

Réponses:

9

L'approche la plus simple serait

  library(raster)
  x <- spdf1 - spdf2

  # or, more formally
  y <- erase(spdf1,  spdf2)

Voir? 'Raster-package' (section XIV) pour plus de fonctions traitant de la superposition de polygones. Ces fonctions utilisent les fonctions de base des rgeos sous le capot, dans les fonctions de "niveau utilisateur" (par opposition aux fonctions "de niveau développeur").

Robert Hijmans
la source
Qu'entendez-vous par "dans les fonctions" au niveau utilisateur "(par opposition aux fonctions" au niveau développeur ")"?
luciano
rgeosfournit des opérations géométriques mais ne traite pas des attributs des données. Ainsi, l'utilisation de ces fonctions nécessite beaucoup de travail pour tout garder ensemble. Les fonctions raster simplifient cela et se comportent comme des fonctions similaires dans le logiciel SIG,
Robert Hijmans
Oui, mais cela peut conduire à une erreur dans SpatialPolygonsDataFrame (part2, x @ data [match (row.names (part2),: row.names of data et Polygons IDs ne correspondent pas
jebyrnes
Ce serait un bug. Pouvez-vous en donner un exemple?
Robert Hijmans
4

Une solution de contournement serait de rajouter les attributs après avoir fait le clip, lors de la conversion de SpatialPolygonsà SpatialPolygonsDataFrame.

sp3 <- gDifference(spdf1, spdf2, byid = TRUE)
row.names(sp3) <- row.names(spdf1)

spdf3 <- SpatialPolygonsDataFrame(sp3, data = spdf1@data)

spdf3@data

   variable1 variable2
p1       232       235
p2       242       464

plot(spdf3, col="red")

entrez la description de l'image ici

Andre Silva
la source
Cela ressemble à un problème que j'ai, seul le clip de sortie dans mon instance particulière produit des noms de domaine à partir de spdf1 qui n'existent pas (comme un simple gsub pour se débarrasser du 2ème chiffre de la ligne. Les noms devraient prendre soin de la correspondance des noms de famille, non? )
jebyrnes
Erreur dans SpatialPolygonsDataFrame (clip, data = as.data.frame (all_spdfs_together @ data)): non-concordance de la longueur des objets: le clip a 1718 objets Polygones, mais as.data.frame (all_spdfs_together @ data) a 86 lignes
jebyrnes
Bien sûr - j'ai un tas de polygones de choses dans l'océan. Certains sont mal placés sur la terre ou se chevauchent avec la terre. Je veux couper cela, donc je n'ai que des zones qui se trouvent dans l'océan. J'ai un fichier de formes de côte à comparer. Voici le code - gist.github.com/jebyrnes/c2e8d2b6c82849dad3a813d952ab8bb0
jebyrnes
1
nevermind - la solution raster :: erase fonctionne (ce n'était pas le cas auparavant pour une raison étrange)
jebyrnes