Trouver efficacement les voisins de 1er ordre de 200k polygones

Pour chacun des 208 781 groupes de blocs du recensement, j'aimerais récupérer les ID FIPS de tous ses voisins de 1er ordre. J'ai toutes les limites de TIGER téléchargées et fusionnées en un seul fichier de formes de 1 Go.

J'ai essayé un script ArcPython qui utilise SelectLayerByLocation pour BOUNDARY_TOUCHES à sa base, mais cela prend plus de 1 seconde pour chaque groupe de blocs, ce qui est plus lent que je ne le souhaiterais. C'est même après avoir limité la recherche SelectLayerByLocation pour bloquer des groupes dans le même état. J'ai trouvé ce script , mais il utilise également SelectLayerByLocation en interne, donc ce n'est pas plus rapide.

La solution n'a pas besoin d'être basée sur Arc - je suis ouvert à d'autres packages, bien que je sois plus à l'aise avec le codage avec Python.

Si vous avez accès à ArcGIS 10.2 for Desktop, ou peut-être plus tôt, je pense que l' outil Polygon Neighbors (Analysis) qui:

Crée une table avec des statistiques basées sur la contiguïté des polygones (chevauchements, arêtes coïncidentes ou nœuds).

peut rendre cette tâche beaucoup plus facile maintenant.

Pour une solution évitant ArcGIS, utilisez pysal . Vous pouvez obtenir les poids directement à partir des fichiers de formes en utilisant:

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

ou

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

Dirigez-vous vers les documents pour plus d'informations.

Juste une mise à jour. Après avoir suivi les conseils de Whuber, j'ai constaté que la matrice de génération de poids spatiaux utilise simplement des boucles et des dictionnaires Python pour déterminer les voisins. J'ai reproduit le processus ci-dessous.

La première partie parcourt chaque sommet de chaque groupe de blocs. Il crée un dictionnaire avec les coordonnées de sommet comme clés et une liste d'ID de groupe de blocs qui ont un sommet à cette coordonnée comme valeur. Notez que cela nécessite un ensemble de données topologiquement soigné, car seul un chevauchement sommet / sommet parfait sera enregistré en tant que relation de voisinage. Heureusement, les fichiers de formes de groupe de blocs TIGER du Census Bureau sont OK à cet égard.

La deuxième partie parcourt à nouveau tous les sommets de chaque groupe de blocs. Il crée un dictionnaire avec les ID de groupe de blocs comme clés et les ID de voisin de ce groupe de blocs comme valeurs.

# 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

Avec le recul, je me rends compte que j'aurais pu utiliser une méthode différente pour la deuxième partie qui ne nécessitait pas de boucler à nouveau dans le fichier de formes. Mais c'est ce que j'ai utilisé, et cela fonctionne assez bien même pour des milliers de groupes de blocs à la fois. Je n'ai pas essayé de le faire avec l'ensemble des États-Unis, mais il peut s'exécuter pour un État entier.

Une alternative pourrait être d'utiliser PostgreSQL et PostGIS . J'ai posé quelques questions sur la façon d'effectuer des calculs similaires sur ce site:

• Calcul du voisin le plus proche dans PostGIS?
• Identification des relations topologiques à l'aide de PostGIS

J'ai trouvé qu'il y avait une courbe d'apprentissage abrupte pour comprendre comment les différentes pièces du logiciel s'emboîtaient, mais je l'ai trouvé merveilleux pour faire des calculs sur de grandes couches vectorielles. J'ai effectué des calculs de voisin le plus proche sur des millions de polygones et cela a été rapide par rapport à ArcGIS.

Juste quelques commentaires ... la méthode esri / ArcGIS utilise actuellement des dictionnaires pour contenir les informations mais les calculs de base sont effectués en C ++ à l'aide de l'outil Polygon Neighbors. Cet outil génère une table qui contient les informations de contiguïté ainsi que des attr facultatifs comme la longueur de la limite partagée. Vous pouvez utiliser l'outil Générer une matrice de pondérations spatiales si vous souhaitez stocker puis réutiliser les informations à plusieurs reprises. Vous pouvez également utiliser cette fonction dans WeightsUtilities pour générer un dictionnaire [accès aléatoire] avec les informations de contiguïté:

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

où inputFC = tout type de classe d'entités surfaciques, masterField est le champ "ID unique" d'entiers et de contiguityType dans {"ROOK", "QUEEN"}.

Des efforts sont déployés à esri pour ignorer l'aspect tabulaire pour les utilisateurs de Python et passer directement à un itérateur qui rendrait de nombreux cas d'utilisation beaucoup plus rapides. PySAL et le package spdep dans R sont des alternatives fantastiques [voir la réponse de radek ] . Je pense que vous devez utiliser des fichiers de formes comme format de données dans ces packages qui sont en accord avec le format d'entrée de ces threads. Je ne sais pas comment ils traitent les polygones qui se chevauchent ainsi que les polygones à l'intérieur des polygones. Générer SWM ainsi que la fonction que j'ai décrite compteront ces relations spatiales comme voisines "ROOK" ET "QUEEN".