Исправление осиротевших отверстий в R

18

Я пытаюсь выполнить объединение на общем поле после объединения двух смежных шейп-файлов. Шейп-файлы заканчиваются как минимум одним узким промежутком между ними. Когда я пытаюсь объединиться, я получаю следующую ошибку:

Ошибка в createPolygonsComment (p): rgeos_PolyCreateComment: потерянная дыра, не удается найти содержащий полигон для дыры с индексом 17

Я загрузил воспроизводимый пример в Dropbox по этой ссылке .

Вот код для воссоздания проблемы:

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

Возвращает:

Ошибка в createPolygonsComment (p): rgeos_PolyCreateComment: потерянная дыра, не удается найти содержащий полигон для дыры с индексом 17

Попытка исправить предложенное здесь и здесь :

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

Это возвращает ту же ошибку, которая приходит от попытки объединения, но с другим порядковым номером:

rgeos_PolyCreateComment: потерянная дыра, не может найти содержащий полигон для дыры с индексом 30

Попытка исправить предложенное в полезном уроке Роджера Биванда

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

Возвращает ту же ошибку с индексом 30, как указано выше.

Другие поднимали эту проблему здесь и здесь , и, хотя приведенные выше решения, по-видимому, работают для некоторых случаев, другие случаи не решаются. Один пользователь использовал QGIS для решения проблемы, а другому было исправлено 2 из 3 пунктов, но окончательного решения не было.

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

Я выполнил инструмент «восстановить геометрию» в ArcGIS, и он исправил проблему, но, похоже, в R. должно быть исправление

Люк Маколей
источник
Без ваших данных трудно сказать, где проблема.
@Pascal, я только что загрузил ссылку dropbox с уменьшенным шейп-файлом в 10 Мб в сжатом виде и 16 Мб в разархивированном виде, чтобы воспроизвести проблему. Я не был уверен, как предоставить данные, так как оригинал был 1,5 ГБ, но мне удалось использовать ArcGIS, чтобы сузить проблему до меньшего файла. Существует ли хороший протокол для создания и обмена воспроизводимыми примерами контролируемого размера?
Люк Маколей
Попытка различных подходов с R не сработала. И Qgis зависает при проверке геометрии.

Ответы:

25

Я проанализировал проблемы геометрии в приложенных данных, и кажется, что это не только, orphaned holesно и geometry validity issues. Это правда, что в orphaned holeнекотором роде проблема достоверности геометрии, но rgeos не обрабатывает ее так же, как и для осиротевших отверстий, вместо простого предупреждения возникает ошибка. Как вы указываете, они являются подсказками для проверки многоугольных отверстий, но это не всегда успешно, когда применяется для исправления осиротевших отверстий.

Итак, начнем:

  1. очистите ваши данные (что необходимо, если вы хотите выполнять геообработку, например, объединение)

  2. используйте очищенные данные в процессе объединения

1. Очистка геометрии Исправление геометрии в R может быть иногда сложным, поэтому я попытался создать экспериментальный пакет R (см. Https://github.com/eblondel/cleangeo ), который предназначен для облегчения очистки spобъектов (в настоящее время ограничен многоугольники). Вы можете установить пакет с помощью:

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

Для начала, хорошо, что вы видите какие проблемы с геометрией у ваших исходных данных. Для этого вы можете выполнить следующее (ваши данные большие, поэтому это может занять некоторое время):

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

При этом вы увидите, что ваши данные имеют 2 вида проблем: orphaned holesи geometry validity issues. Обе (и не только осиротевшие дыры) могут привести к unionсбою процесса, поэтому данные должны быть очищены до этого, при возможности, в автоматизированном режиме. Для быстрого воспроизведения первый пример кода ниже принимает только подмножество функций, которые помечены как подозрительные (за исключением последней, с исходными данными с индексом = 9002 - см. Мое примечание ниже об этом)

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

Если вы clgeo_Cleanхорошо справляетесь с работой, вы должны получить все геометрические фигуры. Вы можете применить это к полному набору данных (кроме индекса объекта = 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. Процесс объединения Теперь, давайте посмотрим, unionработает ли этот набор данных:

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

Примечание: как было сказано ранее, я удалил одну функцию (индекс = 9002). Построив ее plot(sp[9002,]), вы увидите, что эта функция очень (очень) сложна. Я исключил его из образца только потому, что проверка отверстий заняла слишком много времени. Посмотрим теперь, возникает ли такая же проблема при использовании readShapePoly(из maptools) для чтения данных ...

3. Переключитесь на readShapePoly против readOGR для чтения данных (ОБНОВЛЕНИЕ)

readOGRэто не единственная функция, доступная для чтения шейп-файлов. Вы также можете использовать readShapePolyиз maptoolsпакета, как правило, более производительный, чем первый:

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

Помимо бега быстрее:

  • Если вы используете приведенный выше код, основанный на clgeo_CollectionReportпроблемах, проблем с сиротами нет, но проблемы с геометрией остаются.

  • Очистка геометрии с помощью clgeo_Cleanтакже работает хорошо, и теперь она не застревает с индексом функции 9002

  • И ... процесс объединения работает.

Смотрите ниже результат графика:

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

Союзный результат

Вывод : предпочитайте maptools читать ваши данные шейп- файла и подумайте об использовании cleangeo для очистки ваших данных перед любой геообработкой.

eblondel
источник
Спасибо, еблондель! Я попробую это. Спасибо за разработку пакета!
Люк Маколей
Привет, Эблондель, это сработало хорошо, но я хотел сообщить, что при корректировке геометрии он часто создает очень большой многоугольник при работе со сложными и сложными объектами. Например, сеть дорог была скорректирована до большого многоугольника, который в основном был протяженностью сети. Я не уверен, насколько легко это исправить, но хотел, чтобы вы знали.
Люк Маколей
Вау. Очень впечатляющий пакет. Должно быть, это было много работы.
nograpes
3
Спасибо @nograpes за ваш отзыв. Когда эта проблема была опубликована, я создал этот пакет с нуля, а также потому, что очистка геометрии не всегда является легкой задачей. Если вы находитесь на Github, я бы приветствовал вашу «звезду» :-), я хотел бы еще улучшить пакет в будущем и, возможно, выпустить его на CRAN.
eblondel
7
Просто, чтобы вы знали, что cleangeo был опубликован в CRAN ( cran.r-project.org/package=cleangeo ), чтобы все люди, которые его используют, могли свободно сообщать о запросах на улучшение или ошибках в Github.
Eblondel
1

Удобное решение, которое продолжает работать для меня в R, это применить буфер нулевой ширины :

#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 занимает некоторое время с этим набором данных, но, кажется, работает просто отлично.

FGG
источник