Arreglando agujeros huérfanos en R

18

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

  

Error en createPolygonsComment (p): rgeos_PolyCreateComment: agujero huérfano, no se puede encontrar el polígono que contiene el agujero en el índice 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)

Devoluciones:

  

Error en createPolygonsComment (p): rgeos_PolyCreateComment: huérfano   agujero, no se puede encontrar el polígono que contiene para el agujero en el índice 17

Intentando la solución propuesta here y here :

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

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

  

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

Intente la solución propuesta en útil tutorial de Roger Bivand

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

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

Otros han planteado este problema here y here , y mientras Las soluciones planteadas anteriormente parecen funcionar para algunos casos, otros casos no se resuelven. Un usuario usó QGIS para solucionar el problema, y el otro tenía 2 de 3 elementos fijos, pero ninguna solución para el final.

Parece que las personas continúan 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 corrigió el problema, pero parece que debería haber una solución en R.

    
pregunta Luke Macaulay 16.09.2014 - 11:04

1 respuesta

25

He analizado los problemas de geometría en los datos adjuntos, y parece que SÓLO no tiene orphaned holes sino también geometry validity issues . Es cierto que un orphaned hole es de alguna manera un problema de validez de geometría, pero rgeos no lo maneja de la misma manera, como para los huérfanos huecos, se genera un error, en lugar de una simple advertencia. Como indicas, son consejos para verificar los orificios de polígonos, pero no siempre es exitoso cuando se aplican para reparar orificios huérfanos.

Entonces, vamos a:

  1. limpie sus datos (que se requieren si desea hacer geoprocesamiento como unión)

  2. use los datos limpiados con su proceso de unión

1. Geometría de limpieza La reparación de geometrías en R puede ser a veces un desafío, por lo que he intentado construir un paquete experimental R (consulte enlace ) que pretende facilitar limpieza de objetos sp (ahora limitados en formas poligonales). Puedes instalar el paquete con:

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

Para comenzar, es bueno que vea cuáles son los problemas de geometría con sus datos de origen. Para esto, puede ejecutar lo siguiente (sus datos son grandes, por lo que puede tomar algún 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á que sus datos tienen 2 tipos de problemas: orphaned holes y geometry validity issues . Es probable que ambos (y no solo los orificios huérfanos) hagan que el proceso union falle, por lo que los datos deben limpiarse antes, de manera automatizada cuando sea posible. Para una reproducción rápida, el primer código de muestra a continuación solo toma el subconjunto de características que están etiquetadas como sospechosas (excepto la última, con un índice = 9002 en los datos originales; consulte mi nota a continuación)

#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ía obtener todas las geometrías válidas ahora. Puede aplicar esto al conjunto de datos completo (excepto el índice de características = 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 sindical Ahora, veamos si 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 se dijo antes, eliminé una función (índice = 9002). Al trazarla: plot(sp[9002,]) , verás que esta función es muy (muy) compleja. Lo he excluido de la muestra solo porque la comprobación de los agujeros llevaba demasiado tiempo. Veamos ahora si ocurre el mismo problema utilizando readShapePoly (de maptools ) para leer los datos ...

3. Cambie a readShapePoly en lugar de readOGR para leer datos (ACTUALIZAR)

readOGR no es la única función disponible para leer shapefiles. También puede usar readShapePoly del paquete maptools , generalmente más eficiente que el primero:

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

Aparte de correr más rápido:

  • si usa el código anterior basado en clgeo_CollectionReport , no hay problema de orificios huérfanos, pero sí problemas de geometría.

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

  • Y ... el proceso de unión 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:prefiera maptools para leer su shapefile datos, y considere usar cleangeo para limpiar sus datos antes de cualquier geoprocesamiento.

    
respondido por el eblondel 25.09.2014 - 21:07

Lea otras preguntas en las etiquetas