Эффективное нахождение соседей 1-го порядка из 200k полигонов

14

Для каждой из 208 781 группы блоков переписи я хотел бы получить идентификаторы FIPS всех соседей 1-го порядка. Я загрузил все границы TIGER и объединил их в один шейп-файл объемом 1 ГБ.

Я попробовал скрипт ArcPython, который использует SelectLayerByLocation для BOUNDARY_TOUCHES в своей основе, но он занимает более 1 секунды для каждой группы блоков, что медленнее, чем хотелось бы. Это даже после того, как я ограничиваю поиск SelectLayerByLocation, чтобы блокировать группы в том же состоянии. Я нашел этот скрипт , но он также использует SelectLayerByLocation внутри, так что он не быстрее.

Решение не должно быть на основе Arc - я открыт для других пакетов, хотя мне удобнее программировать на Python.

dmahr
источник
2
Начиная с версии 9.3, в наборе инструментов Пространственная статистика были инструменты для этого. Начиная с 10.0, они очень эффективны. Я помню, как выполнял аналогичную операцию для шейп-файла сопоставимого размера (все блоки в одном состоянии), и он завершился за 30 минут, 15 из которых только для дискового ввода-вывода - и это было два года назад на гораздо более медленной машине. Исходный код Python также доступен.
whuber
Какой инструмент геообработки в пространственной статистике вы использовали?
dmahr
1
Я забыл его имя; это специально для создания таблицы отношений соседей полигонов. Справочная система рекомендует вам создать эту таблицу перед запуском любого из инструментов пространственной статистики на основе соседей, чтобы инструментам не приходилось пересчитывать эту информацию на лету при каждом запуске. Существенным ограничением, по крайней мере в версии 9.x, было то, что вывод был в формате .dbf. Для большого входного шейп-файла это не сработает, и в этом случае вам нужно либо разбить операцию на части, либо взломать код Python для вывода в лучшем формате.
whuber
Да это оно. Код Python полностью использует внутренние возможности ArcGIS (которые используют пространственные индексы), что делает алгоритм довольно быстрым.
whuber

Ответы:

3

Если у вас есть доступ к ArcGIS 10.2 for Desktop или, возможно, ранее, я думаю, что инструмент Polygon Neighbors (Analysis) , который:

Создает таблицу со статистикой на основе смежности многоугольников (перекрытия, совпадающие ребра или узлы).

Соседи многоугольника

может сделать эту задачу намного проще сейчас.

PolyGeo
источник
Спасибо, PolyGeo. Я обновил принятый ответ, чтобы инструмент «Соседи многоугольника» получил немного больше информации. Это определенно более надежно, чем мой ручной метод на основе Python, хотя я не уверен, как сравнивается масштабируемость с большими наборами данных.
dmahr
В настоящее время я использую 10.3, и он не работает в моем шейп-файле с ~ 300K блоками переписи.
Soandos
@soandos Звучит так, что, возможно, стоит исследовать / задавать новый вопрос.
PolyGeo
8

Для решения, избегающего ArcGIS, используйте pysal . Вы можете получить веса непосредственно из шейп-файлов, используя:

w = pysal.rook_from_shapefile("../pysal/examples/columbus.shp")

или

w = pysal.queen_from_shapefile("../pysal/examples/columbus.shp")

Отправляйтесь в документы за дополнительной информацией.

Радек
источник
3

Просто обновление. Следуя совету Уубер, я обнаружил, что Матрица генерирования пространственных весов просто использует петли и словари Python для определения соседей. Я воспроизвел процесс ниже.

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

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

# Create dictionary of vertex coordinate : [...,IDs,...]
BlockGroupVertexDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupsDescription.ShapeFieldName
#For every block group...
for BlockGroupItem in BlockGroupCursor :
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex...
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt empty interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                #If coordinate is already in dictionary, append this BG's ID
                if PointText in BlockGroupVertexDictionary:
                    BlockGroupVertexDictionary[PointText].append(BlockGroupID)
                #If coordinate is not already in dictionary, create new list with this BG's ID
                else:
                    BlockGroupVertexDictionary[PointText] = [BlockGroupID]
del BlockGroupItem
del BlockGroupCursor


#Create dictionary of ID : [...,neighbors,...]
BlockGroupNeighborDictionary = {}
BlockGroupCursor = arcpy.SearchCursor(BlockGroups.shp)
BlockGroupDescription = arcpy.Describe(BlockGroups.shp)
BlockGroupShapeFieldName = BlockGroupDescription.ShapeFieldName
#For every block group
for BlockGroupItem in BlockGroupCursor:
    ListOfBlockGroupNeighbors = []
    BlockGroupID = BlockGroupItem.getValue("BKGPIDFP00")
    BlockGroupFeature = BlockGroupItem.getValue(BlockGroupShapeFieldName)
    for BlockGroupPart in BlockGroupFeature:
        #For every vertex
        for BlockGroupPoint in BlockGroupPart:
            #If it exists (and isnt interior hole signifier)...
            if BlockGroupPoint:
                #Create string version of coordinate
                PointText = str(BlockGroupPoint.X)+str(BlockGroupPoint.Y)
                if PointText in BlockGroupVertexDictionary:
                    #Get list of block groups that have this point as a vertex
                    NeighborIDList = BlockGroupVertexDictionary[PointText]
                    for NeighborID in NeighborIDList:
                        #Don't add if this BG already in list of neighbors
                        if NeighborID in ListOfBGNeighbors:
                            pass
                        #Add to list of neighbors (as long as its not itself)
                        elif NeighborID != BlockGroupID:
                            ListOfBGNeighbors.append(NeighborID)
    #Store list of neighbors in blockgroup object in dictionary
    BlockGroupNeighborDictionary[BlockGroupID] = ListOfBGNeighbors

del BlockGroupItem
del BlockGroupCursor
del BlockGroupVertexDictionary

Оглядываясь назад, я понимаю, что мог бы использовать другой метод для второй части, который не требовал бы повторения шейп-файла снова. Но это то, что я использовал, и он работает довольно хорошо даже для тысяч групп блоков одновременно. Я не пробовал делать это со всеми США, но он может быть выполнен для всего штата.

dmahr
источник
2

Альтернативой может быть использование PostgreSQL и PostGIS . Я задал несколько вопросов о том, как выполнить аналогичные вычисления на этом сайте:

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

djq
источник
1

Просто некоторые комментарии ... метод esri / ArcGIS в настоящее время использует словари для хранения информации, но основные вычисления выполняются в C ++ с использованием инструмента Polygon Neighbour Tool. Этот инструмент генерирует таблицу, которая содержит информацию о смежности, а также необязательные атрибуты, такие как длина общей границы. Вы можете использовать инструмент Generate Spatial Weights Matrix Tool, если хотите сохранить, а затем повторно использовать информацию снова и снова. Вы также можете использовать эту функцию в WeightsUtilities для создания словаря [произвольного доступа] с информацией о смежности:

contDict = polygonNeighborDict(inputFC, masterField, contiguityType = "ROOK")

где inputFC = любой класс объектов полигонов любого типа, masterField - это поле "уникальный идентификатор" целых чисел и contiguityType в {"ROOK", "QUEEN"}.

В esri предпринимаются попытки пропустить табличный аспект для пользователей Python и перейти непосредственно к итератору, который значительно ускорил бы многие случаи использования. PySAL и пакет spdep в R - фантастические альтернативы [см . Ответ radek] . Я думаю, что вы должны использовать шейп-файлы в качестве формата данных в этих пакетах, который настроен с этим форматом ввода потоков. Не уверен, как они справляются с перекрывающимися полигонами, а также с полигонами внутри полигонов. Генерация SWM, а также функция, которую я описал, будут считать эти пространственные отношения соседями «ROOK» И «QUEEN».

Марк Яникас
источник