Korzystasz z OGR i Shapely bardziej wydajnie? [Zamknięte]

Szukam sugestii, jak ulepszyć mój kod Python. Zwykle wydajność nie ma dla mnie znaczenia, ale teraz pracuję z plikiem tekstowym lokalizacji w USA z ponad 1,5 miliona punktów. Przy danej konfiguracji uruchomienie operacji w jednym punkcie zajmuje około 5 sekund; Muszę obniżyć tę liczbę.

Używam trzech różnych pakietów GIS Pythona, aby wykonać kilka różnych operacji na punktach i wygenerować nowy plik tekstowy z ogranicznikami.

1. Korzystam z OGR, aby odczytać plik kształtu granicy hrabstwa i uzyskać dostęp do geometrii granicy.
2. Zgrabnie sprawdza, czy punkt znajduje się w którymś z tych hrabstw.
3. Jeśli jest w jednym, używam Biblioteki Shapefile w Pythonie, aby pobrać informacje o atrybutach z granicy .dbf.
4. Następnie piszę informacje z obu źródeł do pliku tekstowego.

Podejrzewam, że nieefektywność polega na posiadaniu 2-3 warstwowej pętli ... nie jestem pewien, co z tym zrobić. Szczególnie szukam pomocy u kogoś, kto ma doświadczenie w korzystaniu z któregokolwiek z tych 3 pakietów, ponieważ po raz pierwszy korzystam z dowolnego z nich.

import os, csv
from shapely.geometry import Point
from shapely.geometry import Polygon
from shapely.wkb import loads
from osgeo import ogr
import shapefile

pointFile = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\NationalFile_20110404.txt"
shapeFolder = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New"
#historicBounds = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\US_Counties_1860s_NAD"
historicBounds = "US_Counties_1860s_NAD"
writeFile = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\NewNational_Gazet.txt"

#opens the point file, reads it as a delimited file, skips the first line
openPoints = open(pointFile, "r")
reader = csv.reader(openPoints, delimiter="|")
reader.next()

#opens the write file
openWriteFile = open(writeFile, "w")

#uses Python Shapefile Library to read attributes from .dbf
sf = shapefile.Reader("C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\US_Counties_1860s_NAD.dbf")
records = sf.records()
print "Starting loop..."

#This will loop through the points in pointFile    
for row in reader:
    print row
    shpIndex = 0
    pointX = row[10]
    pointY = row[9]
    thePoint = Point(float(pointX), float(pointY))
    #This section uses OGR to read the geometry of the shapefile
    openShape = ogr.Open((str(historicBounds) + ".shp"))
    layers = openShape.GetLayerByName(historicBounds)
    #This section loops through the geometries, determines if the point is in a polygon
    for element in layers:
        geom = loads(element.GetGeometryRef().ExportToWkb())
        if geom.geom_type == "Polygon":
            if thePoint.within(geom) == True:
                print "!!!!!!!!!!!!! Found a Point Within Historic !!!!!!!!!!!!"
                print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])
                print records[shpIndex]
                openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))
        if geom.geom_type == "MultiPolygon":
            for pol in geom:
                if thePoint.within(pol) == True:
                    print "!!!!!!!!!!!!!!!!! Found a Point Within MultiPolygon !!!!!!!!!!!!!!"
                    print str(row[1]) + ", " + str(row[2]) + ", " + str(row[5]) + " County, " + str(row[3])
                    print records[shpIndex]
                    openWriteFile.write((str(row[0]) + "|" + str(row[1]) + "|" + str(row[2]) + "|" + str(row[5]) + "|" + str(row[3]) + "|" + str(row[9]) + "|" + str(row[10]) + "|" + str(records[shpIndex][3]) + "|" + str(records[shpIndex][9]) + "|\n"))
        shpIndex = shpIndex + 1
    print "finished checking point"
    openShape = None
    layers = None


pointFile.close()
writeFile.close()
print "Done"

Pierwszym krokiem byłoby przesunięcie pliku kształtu poza pętlę wierszy, otwieranie i zamykanie pliku kształtu 1,5 miliona razy.

Szczerze mówiąc, włożyłbym wszystko do PostGIS i robiłem to przy użyciu SQL na indeksowanych tabelach.

Szybkie spojrzenie na kod przywołuje na myśl kilka optymalizacji:

• Najpierw sprawdź każdy punkt względem obwiedni / obwiedni wielokątów, aby wyeliminować oczywiste wartości odstające. Możesz pójść o krok dalej i policzyć liczbę pól, w których znajduje się punkt, jeśli jest dokładnie jeden, to nie trzeba go testować pod kątem bardziej złożonej geometrii (cóż, tak naprawdę będzie, jeśli leży w więcej niż jeden, będzie musiał zostać przetestowany dalej. Możesz zrobić dwa przejścia, aby wyeliminować proste przypadki ze złożonych przypadków).

• Zamiast zapętlać każdy punkt i testować go pod kątem wielokątów, zapętlaj wielokąty i testuj każdy punkt. Ładowanie / konwersja geometrii jest powolna, więc chcesz to zrobić tak mało, jak to możliwe. Utwórz też początkowo listę Punktów z CSV, ponownie, aby uniknąć konieczności wielokrotnego wykonywania tego zadania na punkt, a następnie odrzucając wyniki na końcu tej iteracji.

• Przestrzennie indeksuj swoje punkty, co wiąże się z konwersją do pliku shapefile, pliku SpatialLite lub czegoś w rodzaju bazy danych PostGIS / PostgreSQL . Ma to tę zaletę, że narzędzia takie jak OGR będą w stanie wykonać większość pracy za Ciebie.

• Nie zapisuj danych wyjściowych do końca: print () jest w najlepszym przypadku kosztowną funkcją. Zamiast tego przechowuj dane jako listę i wypisz je na samym końcu, używając standardowych funkcji wytrawiania w języku Python lub funkcji zrzutu listy.