Cómo arreglar los agujeros huérfanos en R

Question

Estoy tratando de realizar una unión en un campo común después de fusionar dos shapefiles adyacentes. Los shapefiles terminan con al menos una delgada franja de espacio entre ellos. Cuando intento una unión obtengo el siguiente error de agujero huérfano:

Error en createPolygonsComment(p) : rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 17

He subido un ejemplo reproducible a Dropbox en este enlace .

Aquí está el código para recrear el problema:

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

Devuelve:

Error en createPolygonsComment(p) : rgeos_PolyCreateComment: orphaned agujero, no puede encontrar el polígono que contiene el agujero en el índice 17

Probando el arreglo propuesto aquí y aquí :

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

Esto devuelve el mismo error que viene del intento de unión pero con diferente número de índice:

rgeos_PolyCreateComment: agujero huérfano, no puede encontrar el polígono que lo contiene para el agujero en el índice 30

Probando el arreglo propuesto en El útil tutorial de Roger Bivand

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

Devuelve el mismo error en el índice 30 que el anterior.

Otros han planteado este problema aquí y aquí Y aunque las soluciones planteadas anteriormente parecen funcionar para algunos casos, otros no se resuelven. Un usuario utilizó QGIS para solucionar el problema, y el otro tuvo 2 de 3 elementos solucionados, pero sin resolución para el último.

Parece que la gente sigue teniendo problemas a pesar de que este código funciona de vez en cuando. ¿Alguien ha encontrado una solución dentro de R?

He realizado la herramienta "reparar geometría" en ArcGIS, y ha corregido el problema, pero parece que debería haber una solución en R.

Answer 1

He analizado los problemas de geometría en los datos adjuntos, y parece que no SOLO tiene orphaned holes pero también geometry validity issues . Es cierto que un orphaned hole es de alguna manera un problema de validez de la geometría, pero rgeos no lo maneja de la misma manera, ya que para los agujeros huérfanos, se levanta un error, en lugar de una simple advertencia. Como indicas, son pistas para comprobar los agujeros de los polígonos, pero no siempre tiene éxito cuando se aplica para arreglar los agujeros huérfanos.

Así que, vamos:

1. limpiar los datos (lo cual es necesario si se desea hacer un geoprocesamiento como la unión)

2. utilizar los datos depurados con su proceso sindical

1. Geometría de la limpieza La fijación de las geometrías en R puede ser a veces un reto, por lo que he tratado de construir un paquete R experimental (ver https://github.com/eblondel/cleangeo ) que pretende facilitar la limpieza de sp objetos (por ahora limitado a formas poligonales). Puede instalar el paquete con:

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

Para empezar, es bueno que veas cuáles son los problemas de geometría de tus datos de origen. Para ello, puedes ejecutar lo siguiente (tus datos son grandes, así que puede llevar algo de tiempo):

#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,]

Con esto, verás que tus datos tienen 2 tipos de problemas: orphaned holes y geometry validity issues . Ambos (y no sólo los huecos huérfanos) son susceptibles de hacer el union proceso fallido, por lo que los datos deben limpiarse antes, de forma automatizada cuando sea posible. Para una reproducción rápida, el primer código de ejemplo que aparece a continuación sólo toma el subconjunto de características que están etiquetadas como sospechosas (excepto la última, con índice = 9002 en los datos originales - véase mi nota más abajo sobre esto)

#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_Clean hace bien el trabajo, deberías tener todas las geometrías válidas ahora. Puedes aplicar esto al conjunto de datos completo (excepto el índice de característica = 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. Proceso de unión Ahora, veamos si el union funciona en este conjunto de datos:

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

Nota: como he dicho antes, he eliminado una característica (índice = 9002): plot(sp[9002,]) Verás que esta función es muy (muy) compleja. La he excluido de la muestra sólo porque comprobar los agujeros me llevaba demasiado tiempo. Veamos ahora si se produce el mismo problema utilizando readShapePoly (de maptools ) para leer los datos...

3. Cambiar a readShapePoly vs. readOGR para la lectura de datos (UPDATE)

readOGR no es la única función disponible para leer archivos shape. También se puede utilizar readShapePoly de maptools que, por lo general, es más eficaz que el primero:

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

Aparte de correr más rápido:

• si utiliza el código anterior basado en clgeo_CollectionReport No existe el problema de los agujeros huérfanos, pero sí el de la geometría.

• Limpieza de la geometría con clgeo_Clean también funciona bien, y ahora no se atasca con el índice de características 9002

• Y... el proceso sindical funciona.

Vea a continuación el resultado de la trama:

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

Conclusión: : prefiero maptools para leer sus datos shapefile, y considere el uso de cleangeo para limpiar sus datos antes de cualquier geoprocesamiento.

Answer 2

Una solución conveniente que me sigue funcionando en R es aplicar un búfer de ancho cero :

#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 tarda un poco con este conjunto de datos, pero parece funcionar bien.