Correction des trous orphelins dans R

J'essaie d'effectuer une union sur un champ commun après avoir fusionné deux fichiers de formes adjacents. Les fichiers de formes se retrouvent avec au moins un mince ruban d'espace entre eux. Lorsque je tente une union, j'obtiens l'erreur de trou orphelin suivante:

Erreur dans createPolygonsComment (p): rgeos_PolyCreateComment: trou orphelin, impossible de trouver le polygone contenant le trou pour l'index 17

J'ai téléchargé un exemple reproductible sur Dropbox à ce lien .

Voici le code pour recréer le problème:

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

Retour:

Erreur dans createPolygonsComment (p): rgeos_PolyCreateComment: trou orphelin, impossible de trouver le polygone contenant le trou pour l'index 17

Essayer le correctif proposé ici et ici :

slot(example, "polygons") <- lapply(slot(example, "polygons"), checkPolygonsHoles)

Cela renvoie la même erreur qui provient de la tentative d'union mais avec un numéro d'index différent:

rgeos_PolyCreateComment: trou orphelin, impossible de trouver le polygone contenant le trou à l'index 30

Essayer le correctif proposé dans le tutoriel utile de Roger Bivand

fix <- slot(example, "polygons")
fixa <- lapply(fix, checkPolygonsHoles)

Renvoie la même erreur à l'index 30 que ci-dessus.

D'autres ont soulevé ce problème ici et ici , et bien que les solutions proposées ci-dessus semblent fonctionner pour certains cas, d'autres cas ne sont pas résolus. Un utilisateur a utilisé QGIS pour résoudre le problème, et l'autre a eu 2 des 3 éléments corrigés, mais aucune résolution pour le dernier.

Il semble que les gens continuent d'avoir des problèmes malgré le fait que ce code fonctionne de temps en temps. Quelqu'un a-t-il trouvé une solution au sein de R?

J'ai exécuté l'outil "réparer la géométrie" dans ArcGIS, et cela a corrigé le problème, mais il semble qu'il devrait y avoir un correctif dans R.

J'ai analysé les problèmes de géométrie dans les données jointes, et il semble que non seulement orphaned holesmais aussi geometry validity issues. Il est vrai que an orphaned holeest en quelque sorte un problème de validité de la géométrie, mais rgeos ne le gère pas de la même manière, comme pour les trous orphelins, une erreur est déclenchée, au lieu d'un simple avertissement. Comme vous l'avez indiqué, ce sont des astuces pour vérifier les trous de polygone, mais ce n'est pas toujours réussi lorsqu'ils sont appliqués afin de réparer les trous orphelins.

Alors, disons:

1. nettoyer vos données (ce qui est nécessaire si vous souhaitez faire du géotraitement comme union)

2. utiliser les données nettoyées avec votre processus d'union

1. Nettoyage de la géométrie La correction des géométries dans R peut parfois être difficile, j'ai donc essayé de construire un package R expérimental (voir https://github.com/eblondel/cleangeo ) qui vise à faciliter le nettoyage des spobjets (actuellement limité sur formes polygonales). Vous pouvez installer le package avec:

require(devtools)
install_github("eblondel/cleangeo")
require(cleangeo)

Pour commencer, il est bon de voir quels sont les problèmes de géométrie avec vos données source. Pour cela, vous pouvez exécuter ce qui suit (vos données sont volumineuses, cela peut donc prendre un certain temps):

#get a report of geometry validity & issues for a sp spatial object
report <- clgeo_CollectionReport(sp)
summary <- clgeo_SummaryReport(report)
issues <- report[report$valid == FALSE,]

Avec cela, vous verrez que vos données ont 2 types de problèmes: orphaned holeset geometry validity issues. Les deux (et pas seulement les trous orphelins) sont susceptibles de faire unionéchouer le processus, donc les données doivent être nettoyées avant, de manière automatisée lorsque cela est possible. Pour une reproduction rapide, le premier exemple de code ci-dessous ne prend que le sous-ensemble de fonctionnalités marquées comme suspectes (sauf la dernière, avec index = 9002 dans les données d'origine - voir ma note ci-dessous à ce sujet)

#get suspicious features (indexes)
nv <- clgeo_SuspiciousFeatures(report)
mysp <- sp[nv[-14],]

#try to clean data
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

Si clgeo_Cleancela fonctionne bien, vous devriez obtenir toutes les géométries valides maintenant. Vous pouvez l'appliquer à l'ensemble de données complet (sauf l'index d'entité = 9002)

#try to clean data
mysp <- sp[-9002,]
mysp.clean <- clgeo_Clean(mysp, print.log = TRUE)

#check if they are still errors
report.clean <- clgeo_CollectionReport(mysp.clean)
summary.clean <- clgeo_SummaryReport(report.clean)

2. Processus d'union Maintenant, voyons si les uniontravaux sur cet ensemble de données:

#Attempting a UnionSpatialPolygons based on the COUNTY field
mysp.df <- as(mysp, "data.frame")
countycol <- mysp.df$COUNTY
mysp.diss <- unionSpatialPolygons(mysp.clean, countycol)

Remarque: comme indiqué précédemment, j'ai supprimé une fonctionnalité (index = 9002). En la traçant:, plot(sp[9002,])vous verrez que cette fonctionnalité est très (très) complexe. Je l'ai exclu de l'échantillon uniquement parce que la vérification des trous prenait trop de temps. Voyons maintenant si le même problème se produit en utilisant readShapePoly(à partir de maptools) pour lire les données ...

3. Passez à readShapePoly vs readOGR pour lire les données (UPDATE)

readOGRn'est pas la seule fonction disponible pour lire les fichiers de formes. Vous pouvez également utiliser à readShapePolypartir du maptoolspackage, généralement plus performant que le premier:

require(maptools)
mysp <- readShapePoly("ReproducibleExample.shp")

En plus de courir plus vite:

• si vous utilisez le code ci-dessus en fonction clgeo_CollectionReport, il n'y a pas de problème de trous orphelins, mais toujours des problèmes de géométrie.

• Le nettoyage de la géométrie avec clgeo_Cleanfonctionne également bien, et maintenant il ne se bloque pas avec l'indice de fonctionnalité 9002

• Et ... le processus syndical fonctionne.

Voir ci-dessous le résultat du tracé:

#plot the result
plot(mysp, border= "lightgray")
plot(mysp.diss, border="red", add = TRUE)

Conclusion : préférez maptools pour lire vos données de fichiers de formes et envisagez d'utiliser cleangeo pour nettoyer vos données avant tout géotraitement.

Une solution pratique qui continue de fonctionner pour moi dans R consiste à appliquer un tampon de largeur nulle :

#loading required packages
require(sp)    
require(rgdal)
require(maptools)
require(rgeos)

#load example data, set "dsn=" to your working directory or specify the path
example <- readOGR(dsn=".",layer="ReproducibleExample")

#project your data (I'm using California Albers here) and apply a zero-width buffer
example <- spTransform(example, CRS("+init=epsg:3310"))
example <- gBuffer(example, byid = T, width = 0)

#Attempting a UnionSpatialPolygons based on the COUNTY field
example.df <- as(example, "data.frame")
countycol <- example.df$COUNTY
example.diss <- unionSpatialPolygons(example, countycol)

unionSpatialPolygons prend un certain temps avec cet ensemble de données, mais semble fonctionner très bien.