Получение TopologyException: входной геом 1 недействителен из-за самопересечения в R?

24

'TopologyException: Input geom 1 is invalid' ошибка самопересечения, которая возникает из-за неправильной геометрии многоугольника, широко обсуждалась. Однако в Интернете я не нашел удобного решения, которое бы опиралось исключительно на функциональность R.

Например, мне удалось создать объект «SpatialPolygons» из вывода map("state", ...)следующего хорошего ответа Джоша О'Брайена здесь .

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)

состояния

Проблема с этим широко применяемым набором данных заключается в том, что самопересечение происходит в точке, указанной ниже.

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

К сожалению, эта проблема предотвращает дальнейшее использование spydf_states, например, при вызове rgeos::gIntersection. Как я могу решить эту проблему изнутри R?

fdetsch
источник
1
Если вы приблизитесь к этой точке: plot(spydf_states, xlim=c(-122.1,-122.3),ylim=c(38,38.1))вы увидите, что в ней нет «на первый взгляд» - есть самопересечение.
Космонавт

Ответы:

39

Использование буфера нулевой ширины устраняет многие проблемы топологии в R.

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

Однако работа с непроецированными координатами в долготу может привести rgeosк выдаче предупреждений.

Вот расширенный пример, который сначала перепроецируется в проекцию Альберса:

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)

введите описание изображения здесь

hrbrmstr
источник
4
какие-либо дополнительные комментарии / чтения о том, почему gBuffer"взломать" работает?
MichaelChirico
Вы хотите использовать gSimplify, так как он разрушает data.frame и конвертирует SPDF в объект пространственного многоугольника?
wnursal
5
Если вы используете, sfвы также можете использоватьsf::st_buffer(x, dist = 0)
Фил
также работает в некоторых случаях при использованииPostGIS
natsuapo