Encontrar eficientemente los vecinos de primer orden de polígonos de 200k

Para cada uno de los 208.781 grupos de bloque del Censo, me gustaría recuperar las identificaciones FIPS de todos sus vecinos de primer orden. Tengo todos los límites de TIGER descargados y fusionados en un solo archivo de forma de 1GB.

Probé un script de ArcPython que usa SelectLayerByLocation para BOUNDARY_TOUCHES en su núcleo, pero toma más de 1 segundo para cada grupo de bloques, que es más lento de lo que me gustaría. Esto es incluso después de limitar la búsqueda de SelectLayerByLocation para bloquear grupos en el mismo estado. Encontré este script , pero también usa SelectLayerByLocation internamente, por lo que no es más rápido.

La solución no tiene que estar basada en Arc: estoy abierto a otros paquetes, aunque me siento más cómodo codificando con Python.

Si tiene acceso a ArcGIS 10.2 for Desktop, o posiblemente antes, entonces creo que la herramienta Polygon Neighbours (Analysis) que:

Crea una tabla con estadísticas basadas en la contigüidad del polígono (superposiciones, bordes coincidentes o nodos).

puede hacer esta tarea mucho más fácil ahora.

Para una solución que evite ArcGIS, use pysal . Puede obtener los pesos directamente de los archivos de forma usando:

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

o

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

Dirígete a los documentos para obtener más información.

Solo una actualización. Después de seguir los consejos de Whuber, descubrí que la Matriz Generar Pesos Espaciales simplemente usa bucles y diccionarios de Python para determinar vecinos. Reproduje el proceso a continuación.

La primera parte recorre cada vértice de cada grupo de bloques. Crea un diccionario con coordenadas de vértice como las claves y una lista de ID de grupos de bloques que tienen un vértice en esa coordenada como valor. Tenga en cuenta que esto requiere un conjunto de datos topológicamente limpio, ya que solo la superposición perfecta de vértices / vértices se registrará como una relación vecina. Afortunadamente, los archivos de forma de grupo de bloques TIGER de la Oficina del Censo están bien en este sentido.

La segunda parte recorre cada vértice de cada grupo de bloques nuevamente. Crea un diccionario con ID de grupo de bloque como claves y las ID de vecino de ese grupo de bloque como 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

En retrospectiva, me doy cuenta de que podría haber usado un método diferente para la segunda parte que no requirió recorrer el archivo de forma nuevamente. Pero esto es lo que usé, y funciona bastante bien incluso para miles de grupos de bloques a la vez. No he intentado hacerlo con todo Estados Unidos, pero puede ejecutarse para todo un estado.

Una alternativa podría ser usar PostgreSQL y PostGIS . He hecho algunas preguntas sobre cómo realizar cálculos similares en este sitio:

• ¿Cálculo del vecino más cercano en PostGIS?
• Identificar relaciones topológicas usando PostGIS

Descubrí que hay una curva de aprendizaje empinada para descubrir cómo encajan las distintas partes del software, pero me pareció maravilloso para hacer cálculos en capas de vectores grandes. He realizado algunos cálculos de vecinos cercanos en millones de polígonos y ha sido rápido en comparación con ArcGIS.

Solo algunos comentarios ... el método esri / ArcGIS actualmente utiliza diccionarios para almacenar la información, pero los cálculos principales se realizan en C ++ con la herramienta Vecinos poligonales. Esta herramienta genera una tabla que contiene la información de contigüidad, así como atributos opcionales como la longitud del límite compartido. Puede usar la herramienta Generar matriz de pesos espaciales si desea almacenar y posteriormente reutilizar la información una y otra vez. También puede usar esta función en WeightsUtilities para generar un diccionario [acceso aleatorio] con la información de contigüidad:

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

donde inputFC = cualquier tipo de clase de entidad poligonal, masterField es el campo "ID único" de enteros y contiguityType en {"ROOK", "QUEEN"}.

Hay esfuerzos en esri para omitir el aspecto tabular para los usuarios de Python e ir directamente a un iterador que haría que muchos casos de uso sean mucho más rápidos. PySAL y el paquete spdep en R son alternativas fantásticas [ver la respuesta de radek] . Creo que debe usar archivos de forma como el formato de datos en estos paquetes que está en sintonía con este formato de entrada de hilos. No estoy seguro de cómo manejan los polígonos superpuestos, así como los polígonos dentro de los polígonos. Generar SWM así como la función que describí contará esas relaciones espaciales como vecinos "ROOK" Y "QUEEN".