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

14

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.

dmahr
źródło
2
Od wersji 9.3 do tego celu w zestawie narzędzi Spatial Statistics. Począwszy od 10.0, są bardzo wydajne. Pamiętam, że uruchomiłem podobną operację na pliku kształtu o porównywalnym rozmiarze (wszystkie bloki w jednym stanie) i zakończyłem to w 30 minut, 15 z tego tylko dla dyskowych operacji we / wy - i to było dwa lata temu na znacznie wolniejszej maszynie. Kod źródłowy Pythona jest również dostępny.
whuber
Z jakiego narzędzia geoprzetwarzania w statystykach przestrzennych korzystałeś?
dmahr
1
Zapomniałem jego nazwy; służy do tworzenia tabeli relacji sąsiadów wielokąta. System pomocy zachęca do utworzenia tej tabeli przed uruchomieniem któregokolwiek z narzędzi statystyk przestrzennych sąsiada, aby narzędzia te nie musiały przeliczać tych informacji w locie za każdym razem, gdy są uruchamiane. Istotnym ograniczeniem, przynajmniej w wersji 9.x, było to, że dane wyjściowe były w formacie .dbf. W przypadku dużego pliku kształtu wejściowego, który nie będzie działał, w takim przypadku musisz podzielić operację na części lub zhakować kod Pythona, aby uzyskać lepszy format.
whuber
Tak, to jest to. Kod Python w pełni wykorzystuje wewnętrzne możliwości ArcGIS (które wykorzystują indeksy przestrzenne), dzięki czemu algorytm jest dość szybki.
whuber

Odpowiedzi:

3

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).

Wieloboki sąsiedzi

może znacznie ułatwić to zadanie.

PolyGeo
źródło
Dzięki, PolyGeo. Zaktualizowałem zaakceptowaną odpowiedź, więc narzędzie Polygon Neighbours zyskuje nieco więcej ekspozycji. Jest zdecydowanie bardziej solidny niż moja ręczna metoda oparta na języku Python, chociaż nie jestem pewien, jak wygląda skalowalność z dużymi zestawami danych.
dmahr
Obecnie używam 10.3, a mój plik kształtu nie działa z ~ 300 000 blokami spisu.
soandos
@soandos Wygląda na to, że warto zbadać / zadać nowe pytanie.
PolyGeo
8

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.

radek
źródło
3

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.

dmahr
źródło
2

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

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.

djq
źródło
1

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.

Mark Janikas
źródło