Encontrar com eficiência os vizinhos de primeira ordem de 200 mil polígonos

Para cada um dos 208.781 grupos de blocos do Censo, gostaria de recuperar os IDs do FIPS de todos os seus vizinhos de primeira ordem. Eu tenho todos os limites do TIGER baixados e mesclados em um único shapefile de 1 GB.

Eu tentei um script ArcPython que usa SelectLayerByLocation para BOUNDARY_TOUCHES em seu núcleo, mas leva mais de 1 segundo para cada grupo de blocos mais lento do que eu gostaria. Isso ocorre mesmo depois de eu limitar a pesquisa SelectLayerByLocation para bloquear grupos no mesmo estado. Encontrei esse script , mas ele também usa SelectLayerByLocation internamente, para que não seja mais rápido.

A solução não precisa ser baseada em Arc - estou aberto a outros pacotes, apesar de me sentir mais confortável em codificar com Python.

Se você tem acesso ao ArcGIS 10.2 for Desktop, ou possivelmente mais cedo, acho que a ferramenta Polygon Neighbors (Analysis) que:

Cria uma tabela com estatísticas baseadas na contiguidade de polígono (sobreposições, arestas coincidentes ou nós).

pode tornar essa tarefa muito mais fácil agora.

Para uma solução que evita o ArcGIS, use pysal . Você pode obter os pesos diretamente dos shapefiles usando:

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

ou

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

Vá para os documentos para obter mais informações.

Apenas uma atualização. Depois de seguir o conselho de Whuber, descobri que a matriz Generate Spatial Weights simplesmente usa loops e dicionários Python para determinar vizinhos. Reproduzi o processo abaixo.

A primeira parte percorre todos os vértices de cada grupo de blocos. Ele cria um dicionário com coordenadas de vértice como as chaves e uma lista de IDs de grupos de blocos que possuem um vértice nessa coordenada como valor. Observe que isso requer um conjunto de dados topologicamente limpo, pois apenas a sobreposição perfeita de vértices / vértices será registrada como um relacionamento vizinho. Felizmente, os shapefiles do grupo de blocos TIGER do Census Bureau estão bem nesse aspecto.

A segunda parte percorre todos os vértices de cada grupo de blocos novamente. Ele cria um dicionário com os IDs do grupo de blocos como as chaves e os IDs vizinhos do grupo de blocos como os valores.

# 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

Em retrospectiva, percebo que poderia ter usado um método diferente para a segunda parte que não exigisse um loop pelo shapefile novamente. Mas é isso que eu usei e funciona muito bem, mesmo para milhares de grupos de blocos de cada vez. Eu não tentei fazê-lo com os EUA inteiros, mas pode ser executado para um estado inteiro.

Uma alternativa pode ser usar o PostgreSQL e o PostGIS . Fiz algumas perguntas sobre como executar cálculos semelhantes neste site:

• Cálculo de vizinho mais próximo no PostGIS?
• Identificando relacionamentos topológicos usando o PostGIS

Descobri que havia uma curva de aprendizado acentuada para descobrir como as várias partes do software se encaixavam, mas achei maravilhoso fazer cálculos em grandes camadas vetoriais. Fiz alguns cálculos de vizinhos mais próximos em milhões de polígonos e foi rápido em comparação com o ArcGIS.

Apenas alguns comentários ... o método esri / ArcGIS atualmente usa dicionários para armazenar as informações, mas os cálculos principais são feitos em C ++ usando a ferramenta Polygon Neighbors. Essa ferramenta gera uma tabela que contém as informações de contiguidade, bem como atributos opcionais, como o comprimento do limite compartilhado. Você pode usar a Ferramenta de matriz Gerar pesos espaciais se desejar armazenar e subsequentemente reutilizar as informações repetidamente. Você também pode usar esta função no WeightsUtilities para gerar um dicionário [acesso aleatório] com as informações de contiguidade:

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

onde inputFC = qualquer tipo de classe de recurso de polígono, masterField é o campo "ID exclusivo" de números inteiros e contiguityType em {"ROOK", "QUEEN"}.

Há esforços na esri para ignorar o aspecto tabular dos usuários do Python e ir diretamente para um iterador que tornaria muitos casos de uso muito mais rápidos. PySAL e o pacote spdep em R são alternativas fantásticas [veja a resposta de radek] . Eu acho que você é obrigado a usar shapefiles como o formato de dados nesses pacotes que está em sintonia com este formato de entrada de threads. Não sabe ao certo como eles lidam com polígonos sobrepostos e polígonos dentro de polígonos. Gere SWM, bem como a função que descrevi contará essas relações espaciais como vizinhos "ROOK" E "QUEEN".