Skuteczne znajdowanie sąsiadów 1. rzędu 200 tys. Wielokątów

Dla każdej z 208,781 grup bloków spisu chciałbym odzyskać identyfikatory FIPS wszystkich sąsiadów pierwszego rzędu. Mam wszystkie granice TIGER pobrane i scalone w jeden plik kształtu 1GB.

Próbowałem skryptu ArcPython, który w swoim rdzeniu używa SelectLayerByLocation dla BOUNDARY_TOUCHES, ale zajmuje to ponad 1 sekundę dla każdej grupy bloków, co jest wolniejsze niż chciałbym. Dzieje się tak nawet po tym, jak ograniczę wyszukiwanie SelectLayerByLocation do blokowania grup w tym samym stanie. Znalazłem ten skrypt , ale korzysta on również wewnętrznie z SelectLayerByLocation, więc nie jest już szybszy.

Rozwiązanie nie musi być oparte na Arc - jestem otwarty na inne pakiety, chociaż najwygodniej koduję w Pythonie.

Jeśli masz dostęp do ArcGIS 10.2 dla komputerów stacjonarnych lub być może wcześniej, myślę, że narzędzie Polygon Neighbours (Analysis) , które:

Tworzy tabelę ze statystykami opartymi na ciągłości wielokątów (zakładki, zbieżne krawędzie lub węzły).

może znacznie ułatwić to zadanie.

Aby uzyskać rozwiązanie unikające ArcGIS, użyj pysal . Możesz uzyskać wagi bezpośrednio z plików kształtów, używając:

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

lub

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

Przejdź do dokumentacji, aby uzyskać więcej informacji.

Tylko aktualizacja. Po podążeniu za radą Whubera odkryłem, że macierz generowania wag przestrzennych po prostu wykorzystuje pętle i słowniki Pythona do określania sąsiadów. Proces odtworzyłem poniżej.

Pierwsza część przechodzi przez każdy wierzchołek każdej grupy bloków. Tworzy słownik ze współrzędnymi wierzchołków jako kluczami i listą identyfikatorów grup bloków, które mają wierzchołek na tej współrzędnej jako wartość. Zauważ, że wymaga to uporządkowanego topologicznie zestawu danych, ponieważ tylko idealne nachylenie wierzchołka / wierzchołka zostanie zarejestrowane jako relacja sąsiada. Na szczęście pliki kształtów grup bloków TIGER Biura Spisu Ludności są w tym względzie OK.

Druga część ponownie zapętla się w każdym wierzchołku każdej grupy bloków. Tworzy słownik z identyfikatorami grup bloków jako kluczami i identyfikatorami sąsiadów grupy bloków jako wartości.

# 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

Z perspektywy czasu zdaję sobie sprawę, że mogłem zastosować inną metodę do drugiej części, która nie wymagała ponownego zapętlania pliku kształtu. Ale tego właśnie użyłem i działa całkiem dobrze nawet dla 1000 grup bloków jednocześnie. Nie próbowałem tego robić z całymi USA, ale można to wykonać dla całego stanu.

Alternatywą może być użycie PostgreSQL i PostGIS . Zadałem kilka pytań na temat wykonywania podobnych obliczeń na tej stronie:

• Obliczenia najbliższego sąsiada w PostGIS?
• Identyfikacja zależności topologicznych za pomocą PostGIS

Odkryłem, że istnieje stroma krzywa uczenia się, aby dowiedzieć się, jak różne elementy oprogramowania pasują do siebie, ale uważam, że jest to wspaniałe do wykonywania obliczeń na dużych warstwach wektorowych. Przeprowadziłem obliczenia najbliższego sąsiada na milionach wielokątów i było to szybkie w porównaniu do ArcGIS.

Tylko kilka komentarzy ... metoda esri / ArcGIS używa obecnie słowników do przechowywania informacji, ale podstawowe obliczenia są wykonywane w C ++ za pomocą narzędzia Polygon Neighbours Tool. To narzędzie generuje tabelę, która zawiera informacje o ciągłości, a także opcjonalne atrybuty, takie jak długość wspólnej granicy. Możesz użyć narzędzia Wygeneruj macierz wag przestrzennych, jeśli chcesz przechowywać, a następnie ponownie wykorzystywać informacje w kółko. Tej funkcji można także użyć w WeightsUtilities do wygenerowania słownika [dostęp losowy] z informacjami o przyległości:

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

gdzie inputFC = dowolny typ klasy cech wielokąta, masterField jest polem liczb całkowitych „niepowtarzalny identyfikator” i contiguityType w {„ROOK”, „QUEEN”}.

W esri starano się pominąć aspekt tabelaryczny dla użytkowników Pythona i przejść bezpośrednio do iteratora, który znacznie przyspieszyłby wiele przypadków użycia. PySAL i pakiet spdep w R są fantastycznymi alternatywami [patrz odpowiedź radka ] . Myślę, że musisz użyć plików shapefile jako formatu danych w tych pakietach, który jest dostosowany do tego formatu wejściowego wątków. Nie jestem pewien, jak radzą sobie z nakładającymi się wielokątami, a także wielokątami wewnątrz wielokątów. Wygeneruj SWM, a także funkcja, którą opisałem, policzą te relacje przestrzenne jako „ROOK” I „QUEEN” Neighbours.