Obtención de excepción de topología: ¿La entrada geom 1 no es válida debido a la intersección automática en R?

24

El error de auto-intersección 'TopologyException: Input geom 1 is invalid' que surge de geometrías de polígonos inválidas ha sido ampliamente discutido. Sin embargo, no he encontrado una solución conveniente en la web que se base únicamente en la funcionalidad R.

Por ejemplo, he logrado crear un objeto 'SpatialPolygons' a partir del resultado de map("state", ...)seguir la buena respuesta de Josh O'Brien aquí .

library(maps)
library(maptools)

map_states = map("state", fill = TRUE, plot = FALSE)

IDs = sapply(strsplit(map_states$names, ":"), "[[", 1)
spydf_states = map2SpatialPolygons(map_states, IDs = IDs, proj4string = CRS("+init=epsg:4326"))

plot(spydf_states)

estados

El problema con este conjunto de datos ampliamente aplicado es ahora que la auto-intersección ocurre en el punto que se indica a continuación.

rgeos::gIsValid(spydf_states)
[1] FALSE
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point -122.22023214285259 38.060546477866055

Desafortunadamente, este problema impide el uso posterior de 'spydf_states', por ejemplo, al llamar rgeos::gIntersection. ¿Cómo puedo resolver este problema desde R?

fdetsch
fuente
1
Si hace zoom alrededor de ese punto: plot(spydf_states, xlim=c(-122.1,-122.3),ylim=c(38,38.1))verá que no hay "aparentemente" al respecto, hay una auto-intersección.
Spacedman

Respuestas:

39

El uso de un búfer de ancho cero limpia muchos problemas de topología en R.

spydf_states <- gBuffer(spydf_states, byid=TRUE, width=0)

Sin embargo, trabajar con coordenadas lat-long no proyectadas puede provocar rgeosadvertencias.

Aquí hay un ejemplo extendido que reproyecta primero a una proyección de Albers:

library(sp)
library(rgeos)

load("~/Dropbox/spydf_states.RData")

# many geos functions require projections and you're probably going to end
# up plotting this eventually so we convert it to albers before cleaning up
# the polygons since you should use that if you are plotting the US
spydf_states <- spTransform(spydf_states, 
                            CRS("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=37.5 +lon_0=-96"))

# simplify the polgons a tad (tweak 0.00001 to your liking)
spydf_states <- gSimplify(spydf_states, tol = 0.00001)

# this is a well known R / GEOS hack (usually combined with the above) to 
# deal with "bad" polygons
spydf_states <- gBuffer(spydf_states, byid=TRUE, width=0)

# any bad polys?
sum(gIsValid(spydf_states, byid=TRUE)==FALSE)

## [1] 0

plot(spydf_states)

ingrese la descripción de la imagen aquí

hrbrmstr
fuente
44
¿algún comentario / lectura adicional sobre por qué funciona el gBuffer"hack"?
MichaelChirico
¿Desea usar gSimplify ya que destruye data.frame y convierte el SPDF en un objeto de polígono espacial?
wnursal
55
Si está usando sftambién puede usarsf::st_buffer(x, dist = 0)
Phil
también funciona en algunos casos cuando se usaPostGIS
natsuapo