Effizientes Finden der Nachbarn 1. Ordnung von 200.000 Polygonen

14

Für jede der 208.781 Census-Blockgruppen möchte ich die FIPS-IDs aller Nachbarn erster Ordnung abrufen. Ich habe alle TIGER-Grenzen heruntergeladen und in einem einzigen 1-GB-Shapefile zusammengeführt.

Ich habe ein ArcPython-Skript ausprobiert, das im Kern SelectLayerByLocation für BOUNDARY_TOUCHES verwendet, aber für jede Blockgruppe, die langsamer ist als gewünscht, dauert es mehr als 1 Sekunde. Dies ist auch dann der Fall, wenn ich die Suche in SelectLayerByLocation auf das Blockieren von Gruppen im selben Status beschränkt habe. Ich habe dieses Skript gefunden , aber es verwendet SelectLayerByLocation auch intern, sodass es nicht schneller ist.

Die Lösung muss nicht auf Arc basieren - ich bin offen für andere Pakete, obwohl ich mit Python am besten programmieren kann.

dmahr
quelle
2
Seit Version 9.3 gibt es Tools in der Spatial Statistics-Toolbox, um dies zu tun. Ab 10.0 sind sie sehr effizient. Ich erinnere mich, dass eine ähnliche Operation mit einem Shapefile vergleichbarer Größe (alle Blöcke in einem Zustand) ausgeführt wurde und in 30 Minuten abgeschlossen war, 15 davon nur für Festplatten-E / A - und dies war vor zwei Jahren auf einer viel langsameren Maschine. Auf den Python-Quellcode kann ebenfalls zugegriffen werden.
whuber
Welches Geoverarbeitungswerkzeug in Spatial Statistics haben Sie verwendet?
dmahr
1
Ich vergesse seinen Namen; Es dient speziell zum Erstellen einer Tabelle mit Polygonnachbarnbeziehungen. Das Hilfesystem fordert Sie dazu auf, diese Tabelle zu erstellen, bevor Sie eines der auf Nachbarn basierenden räumlichen Statistik-Tools ausführen, damit die Tools diese Informationen nicht jedes Mal neu berechnen müssen, wenn sie ausgeführt werden. Zumindest in der 9.x-Version bestand eine wesentliche Einschränkung darin, dass die Ausgabe im DBF-Format erfolgte. Bei einem großen Eingabe-Shapefile funktioniert dies nicht. In diesem Fall müssen Sie entweder die Operation in Teile zerlegen oder den Python-Code hacken, um ein besseres Format auszugeben.
whuber
Ja das ist es. Der Python-Code nutzt die internen ArcGIS-Funktionen (die räumliche Indizes verwenden) vollständig aus, wodurch der Algorithmus recht schnell wird.
whuber

Antworten:

3

Wenn Sie Zugriff auf ArcGIS 10.2 for Desktop oder eine frühere Version haben, bietet das Tool " Polygon-Nachbarn (Analyse)" Folgendes:

Erstellt eine Tabelle mit Statistiken auf der Grundlage der Polygonkontiguität (Überlappungen, zusammenfallende Kanten oder Knoten).

Polygon Nachbarn

kann diese Aufgabe jetzt viel einfacher machen.

PolyGeo
quelle
Danke, PolyGeo. Ich habe die akzeptierte Antwort aktualisiert, damit das Werkzeug "Polygon-Nachbarn" mehr Aufmerksamkeit erhält. Es ist definitiv robuster als meine manuelle Python-basierte Methode, obwohl ich nicht sicher bin, wie die Skalierbarkeit mit großen Datenmengen verglichen wird.
dmahr
Ich verwende derzeit 10.3 und es schlägt in meinem Shapefile mit ca. 300.000 Zensusblöcken fehl.
Soandos
@soandos Das hört sich so an, als ob es sich lohnt, eine neue Frage zu stellen.
PolyGeo
8

Verwenden Sie pysal , um eine Lösung zu finden, die ArcGIS vermeidet . Sie können die Gewichte direkt aus Shapefiles abrufen, indem Sie Folgendes verwenden:

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

oder

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

Gehen Sie für die Dokumentation für weitere Informationen.

radek
quelle
3

Nur ein Update. Nachdem ich Whubers Rat gefolgt war, stellte ich fest, dass die Matrix zum Generieren von Raumgewichten einfach Python-Schleifen und -Wörterbücher verwendet, um Nachbarn zu bestimmen. Ich habe den folgenden Prozess reproduziert.

Der erste Teil durchläuft jeden Scheitelpunkt jeder Blockgruppe. Es erstellt ein Wörterbuch mit Scheitelpunktkoordinaten als Schlüsseln und einer Liste von Blockgruppen-IDs, die einen Scheitelpunkt an dieser Koordinate als Wert haben. Beachten Sie, dass hierfür ein topologisch sauberer Datensatz erforderlich ist, da nur eine perfekte Vertex / Vertex-Überlappung als Nachbarbeziehung registriert wird. Glücklicherweise sind die Shapefiles der TIGER-Blockgruppe des Census Bureau in dieser Hinsicht in Ordnung.

Der zweite Teil durchläuft erneut jeden Scheitelpunkt jeder Blockgruppe. Es erstellt ein Wörterbuch mit Blockgruppen-IDs als Schlüsseln und den Nachbar-IDs dieser Blockgruppe als Werten.

# 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

Im Nachhinein stelle ich fest, dass ich für den zweiten Teil eine andere Methode hätte verwenden können, bei der das Shapefile nicht erneut durchlaufen werden musste. Aber dies ist, was ich verwendet habe, und es funktioniert ziemlich gut, auch für Tausende von Blockgruppen gleichzeitig. Ich habe es nicht mit den ganzen USA versucht, aber es kann für einen ganzen Staat durchgeführt werden.

dmahr
quelle
2

Eine Alternative könnte die Verwendung von PostgreSQL und PostGIS sein . Ich habe ein paar Fragen zur Durchführung ähnlicher Berechnungen auf dieser Site gestellt:

Ich fand, dass es eine steile Lernkurve gibt, um herauszufinden, wie die verschiedenen Teile der Software zusammenpassen, aber ich fand es wunderbar, um Berechnungen auf großen Vektorebenen durchzuführen. Ich habe einige Berechnungen zum nächsten Nachbarn für Millionen von Polygonen durchgeführt und diese wurden schnell mit ArcGIS verglichen.

djq
quelle
1

Nur ein paar Kommentare ... Die esri / ArcGIS-Methode verwendet derzeit Wörterbücher, um die Informationen zu speichern, aber die Kernberechnungen werden in C ++ mit dem Polygon Neighbors Tool durchgeführt. Dieses Tool generiert eine Tabelle, die die Informationen zur Kontiguität sowie optionale Attribute wie die Länge der gemeinsam genutzten Grenze enthält. Sie können das Generate Spatial Weights Matrix Tool verwenden, wenn Sie die Informationen immer wieder speichern und anschließend wiederverwenden möchten. Sie können diese Funktion auch in WeightsUtilities verwenden, um ein Wörterbuch [Direktzugriff] mit den folgenden Informationen zu erstellen:

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

Wenn inputFC eine beliebige Art von Polygon-Feature-Class ist, ist masterField das Feld "eindeutige ID" für Ganzzahlen und contiguityType in {"ROOK", "QUEEN"}.

Es gibt Bestrebungen bei esri, den tabellarischen Aspekt für Python-Benutzer zu überspringen und direkt zu einem Iterator zu wechseln, der viele Anwendungsfälle weitaus schneller machen würde. PySAL und das spdep-Paket in R sind fantastische Alternativen [siehe radeks Antwort] . Ich denke, Sie müssen Shapefiles als Datenformat in diesen Paketen verwenden, das mit dem Eingabeformat dieses Threads übereinstimmt. Nicht sicher, wie sie mit überlappenden Polygonen sowie Polygonen innerhalb von Polygonen umgehen. Durch das Generieren von SWM sowie die von mir beschriebene Funktion werden diese räumlichen Beziehungen als "ROOK" - UND "QUEEN" -Nachbarn gezählt.

Mark Janikas
quelle