Usprawnienie kodu Pythona dla dużych zbiorów danych

Mam kod w języku Python, który został zaprojektowany do pobierania punktowych plików kształtów poprzez następujący przepływ pracy:

1. Scal punkty
2. Zintegruj punkty, tak aby wszystkie punkty w odległości 1 m od siebie stały się jednym punktem
3. Utwórz warstwę obiektów, w której zaznaczone są punkty z <10
4. Punkty buforowe
5. Rozdzielczość wielokąta do rastra 1m
6. Przeklasyfikuj, gdzie 1–9 = 1; NoData = 0

Każdy plik kształtu ma około 250 000 do 350 000 punktów obejmujących ~ 5x7 km. Dane punktowe użyte jako dane wejściowe reprezentują lokalizacje drzew. Każdy punkt (tj. Drzewo) ma przypisaną wartość „z”, która reprezentuje promień korony i jest używana w procesie buforowania. Moim zamiarem jest wykorzystanie ostatecznego wyjścia binarnego w osobnym procesie do stworzenia rastra opisującego osłonę czaszy.

Przeprowadziłem test z czterema plikami kształtu, który wygenerował raster 700 MB i zajął 35 minut (procesor i5 i 8 GB pamięci RAM). Widząc, że będę musiał uruchomić ten proces na plikach kształtowych 3500, byłbym wdzięczny za wszelkie porady w celu usprawnienia tego procesu (patrz załączony kod). Ogólnie rzecz biorąc, jaki jest najlepszy sposób radzenia sobie z geoprzetwarzaniem dużych zbiorów danych? Mówiąc dokładniej, czy są jakieś poprawki w kodzie lub przepływie pracy, które mogłyby pomóc w zwiększeniu wydajności?

Edytuj :

Czas (% całości) na zadania geoprzetwarzania:

• Scalanie = 7,6%
• Zintegruj = 7,1%
• Cecha Lyr = 0
• Bufor = 8,8%
• Poly na Raster = 74,8%
• Przeklasyfikowanie = 1,6%

# Import arcpy module
import arcpy

# Check out any necessary licenses
arcpy.CheckOutExtension("spatial")

# Script arguments
temp4 = arcpy.GetParameterAsText(0)
if temp4 == '#' or not temp4:
    temp4 = "C:\\gdrive\\temp\\temp4" # provide a default value if unspecified

Reclassification = arcpy.GetParameterAsText(1)
if Reclassification == '#' or not Reclassification:
    Reclassification = "1 9 1;NODATA 0" # provide a default value if unspecified

Multiple_Value = arcpy.GetParameterAsText(2)
if Multiple_Value == '#' or not Multiple_Value:
    Multiple_Value = "C:\\t1.shp;C:\\t2.shp;C:\\t3.shp;C:\\t4.shp" # provide a default value if unspecified

# Local variables:
temp_shp = Multiple_Value
Output_Features = temp_shp
temp2_Layer = Output_Features
temp_Buffer = temp2_Layer
temp3 = temp_Buffer

# Process: Merge
arcpy.Merge_management(Multiple_Value, temp_shp, "x \"x\" true true false 19 Double 0 0 ,First,#,C:\\#########omitted to save space

# Process: Integrate
arcpy.Integrate_management("C:\\gdrive\\temp\\temp.shp #", "1 Meters")

# Process: Make Feature Layer
arcpy.MakeFeatureLayer_management(temp_shp, temp2_Layer, "z <10", "", "x x VISIBLE NONE;y y VISIBLE NONE;z z VISIBLE NONE;Buffer Buffer VISIBLE NONE")

# Process: Buffer
arcpy.Buffer_analysis(temp2_Layer, temp_Buffer, "z", "FULL", "ROUND", "NONE", "")

# Process: Polygon to Raster
arcpy.PolygonToRaster_conversion(temp_Buffer, "BUFF_DIST", temp3, "CELL_CENTER", "NONE", "1")

# Process: Reclassify
arcpy.gp.Reclassify_sa(temp3, "Value", Reclassification, temp4, "DATA")

Niektóre zmiany algorytmu, które powinny ci pomóc.

Dokonaj wyboru przed scaleniem lub scaleniem. Pozwoli to znacznie ograniczyć późniejsze funkcje, które są najdroższe.

Scalanie i integracja są kosztowne dla pamięci, więc chcesz nadal eliminować funkcje, wprowadzając klasy elementów, i próbuj wykonywać scalenia w drzewie binarnym, aby zachować rozmiar scalania i integracji. np. dla czterech plików kształtów scalasz dwa pliki kształtowe i integrujesz; połącz dwa kolejne pliki kształtów i zintegruj; Scal dwie wynikowe klasy obiektów i zintegruj.

Twoja kolejka zadań zaczyna się jako kolejka odniesień do pliku kształtu. Masz również kolejkę wyników, w której możesz umieścić wyniki. Metoda run () dla procesora przetwarzania równoległego wykona następujące operacje: Zdejmij dwa elementy z kolejki. Jeśli nie zostanie pobrany żaden element (kolejka jest pusta), zakończ proces roboczy. Jeśli jeden element zostanie zabrany, umieść go bezpośrednio w kolejce wyników.

Jeśli zostaną pobrane dwa elementy, dla każdego elementu: jeśli jest to plik kształtu, wybierz dla z <10 i utwórz klasę obiektów in_memory; w przeciwnym razie jest to już klasa funkcji in_memory i pomija krok wyboru. Scal dwie klasy funkcji in_memory, aby utworzyć nową klasę funkcji in_memory. Usuń oryginalne dwie klasy obiektów. Wykonaj integrację w nowej klasie funkcji. Umieść tę klasę obiektów w kolejce wyników.

Następnie uruchom zewnętrzną pętlę while. Pętla rozpoczyna się od kolejki shapefile i sprawdza długość większą niż 1. Następnie uruchamia kolejkę przez proces roboczy. Jeśli kolejka wynikowa ma długość większą niż 1, pętla while wykonuje inne równoległe przetwarzanie przebiegające przez procesy robocze, aż kolejka wynikowa będzie miała klasę 1 in_memory.

np. jeśli zaczynasz z 3500 plików kształtów, twoja pierwsza kolejka będzie miała 3500 zadań. Drugi będzie miał 1750 miejsc pracy. 875, 438, 219, 110, 55, 28, 14, 7, 4, 2, 1. Twoim dużym wąskim gardłem będzie pamięć. Jeśli nie masz wystarczającej ilości pamięci (w takim przypadku zabraknie pamięci podczas tworzenia pierwszej kolejki wyników), zmodyfikuj algorytm, aby scalić więcej niż 2 klasy elementów na raz, a następnie zintegruj, co zmniejszy rozmiar pierwszej kolejki wyników w zamian za dłuższy czas przetwarzania. Opcjonalnie możesz pisać pliki wyjściowe i pomijać za pomocą klas funkcji in_memory. Spowolni to znacznie, ale ominie wąskie gardło pamięci.

Dopiero po przeprowadzeniu scalania i integracji na wszystkich plikach kształtów, kończących się na jednej klasie elementów, wykonujesz bufor, od poly do raster i ponownie klasyfikujesz. W ten sposób te trzy operacje są wykonywane tylko raz, a twoja geometria jest prosta.

Pierwszą rzeczą, którą bym zrobił, to monitorowanie wykorzystania zasobów twojego systemu za pomocą Monitora zasobów w Windows 7 lub perfmon w Vista / XP, aby przekonać się, czy jesteś związany z procesorem , pamięcią lub IO .

Jeśli jesteś związany z pamięcią lub IO, prawdopodobnie niewiele możesz zrobić, ale zaktualizuj sprzęt, zmniejsz rozmiar problemu lub całkowicie zmień podejście.

Jeśli stwierdzisz, że jesteś związany z procesorem, eksperymentowałbym z multiprocessingmodułem lub jednym z wielu dostępnych pakietów przetwarzania równoległego opartych na języku Python , aby sprawdzić, czy możesz użyć większej liczby rdzeni procesora, aby przyspieszyć swoje operacje.

Sztuką wieloprocesowości i równoległości jest ogólnie znalezienie dobrego schematu partycjonowania, który:

1. Pozwala podzielić dane wejściowe na mniejsze zestawy robocze, a następnie ponownie połączyć wyniki w sposób, który ma sens,
2. Dodaje najmniejszy narzut (niektóre są nieuniknione w porównaniu do przetwarzania szeregowego), oraz
3. Pozwala dostosować rozmiar zestawu roboczego, aby jak najlepiej wykorzystać zasoby systemu w celu uzyskania optymalnej wydajności.

Możesz użyć skryptu, który utworzyłem w tej odpowiedzi, jako punktu wyjścia: Kod Porting Avenue do tworzenia Building Shadows do ArcPy / Python dla ArcGIS Desktop?

Zobacz także ten post na blogu ESRI Geoprocessing na ten temat: Python Multiprocessing - Podejścia i uwagi

Myślę, że twoja sprawa będzie jeszcze trudniejsza ze względu na bardziej „czarną skrzynkę” narzędzi, których używasz, a nie z bardziej precyzyjnych układów geometrii, z którymi pracowałem. Być może przydaje się praca z tablicami NumPy .

Natknąłem się również na interesujący materiał do czytania, jeśli chcesz spojrzeć poza arkadę:

• Algorytmy przetwarzania równoległego dla GIS . CRC Press, 1997. Richard Healey, Steve Dowers, Bruce Gittings, Mike J Mineter
• Przyspieszenie przetwarzania rastrowego z drobną i zgrubną równoległością ziarna w GRASS . Materiały z konferencji użytkowników FOSS / GRASS 2004 - Bangkok, Tajlandia, 12-14 września 2004 r. Onil Nazra Persada, Thierry Goubier.
• Przyspieszenie przetwarzania wsadowego przestrzennej analizy rastrowej za pomocą GPU . Computers & Geosciences, tom 45, sierpień 2012, strony 212–220. Mathias Steinbach, Reinhard Hemmerling