Obtention de TopologyException: l'entrée geom 1 n'est pas valide, ce qui est dû à l'auto-intersection dans R?

24

L'erreur d'auto-intersection «TopologyException: Input geom 1 is invalid» qui provient de géométries de polygones invalides a été largement discutée. Cependant, je n'ai pas trouvé de solution pratique sur le Web qui repose uniquement sur la fonctionnalité R.

Par exemple, j'ai réussi à créer un objet «SpatialPolygons» à partir de la sortie de la map("state", ...)suite de la belle réponse de Josh O'Brien ici .

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)

États

Le problème avec cet ensemble de données largement appliqué est maintenant que l'auto-intersection se produit au point indiqué ci-dessous.

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

Malheureusement, ce problème empêche toute autre utilisation de «spydf_states», par exemple lors de l'appel rgeos::gIntersection. Comment puis-je résoudre ce problème depuis R?

fdetsch
la source
1
Si vous effectuez un zoom avant autour de ce point: plot(spydf_states, xlim=c(-122.1,-122.3),ylim=c(38,38.1))vous verrez qu'il n'y a pas "apparemment" à ce sujet - il y a une auto-intersection.
Spacedman

Réponses:

39

L'utilisation d'un tampon de largeur nulle nettoie de nombreux problèmes de topologie dans R.

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

Cependant, travailler avec des coordonnées lat non projetées peut provoquer des rgeosavertissements.

Voici un exemple étendu qui reprojete d'abord une projection d'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)

entrez la description de l'image ici

hrbrmstr
la source
4
des commentaires / lectures supplémentaires sur la raison pour laquelle le gBuffer"hack" fonctionne?
MichaelChirico
voulez-vous utiliser gSimplify car il détruit le data.frame et convertit le SPDF en objet polygonal spatial?
wnursal
5
Si vous utilisez, sfvous pouvez également utilisersf::st_buffer(x, dist = 0)
Phil
fonctionne également dans certains cas lors de l'utilisationPostGIS
natsuapo