OGR und Shapely effizienter nutzen? [geschlossen]

Ich suche einige Vorschläge, wie ich meinen Python-Code effizienter gestalten kann. Normalerweise spielt Effizienz für mich keine Rolle, aber ich arbeite jetzt mit einer Textdatei von US-Standorten mit über 1,5 Millionen Punkten. Mit der gegebenen Konfiguration dauert es ungefähr 5 Sekunden, um Operationen an einem Punkt auszuführen. Ich muss diese Zahl deutlich senken.

Ich verwende drei verschiedene Python-GIS-Pakete, um einige verschiedene Operationen an den Punkten durchzuführen und eine neue Textdatei mit Trennzeichen auszugeben.

1. Ich verwende OGR, um ein County Boundary Shapefile zu lesen und Zugriff auf die Grenzgeometrie zu erhalten.
2. Formschön überprüft, ob sich ein Punkt in einem dieser Landkreise befindet.
3. Wenn es innerhalb einer liegt, verwende ich die Python-Shapefile-Bibliothek, um Attributinformationen von der Grenze .dbf abzurufen.
4. Ich schreibe dann einige Informationen aus beiden Quellen in eine Textdatei.

Ich vermute, dass die Ineffizienz darin liegt, eine 2-3-stufige Schleife zu haben ... nicht ganz sicher, was man dagegen tun soll. Ich bin besonders auf der Suche nach Hilfe für jemanden, der Erfahrung mit diesen 3 Paketen hat, da ich zum ersten Mal eines davon benutze.

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"

Der erste Schritt wäre, das Shapefile außerhalb der Zeilenschleife zu öffnen. Sie öffnen und schließen das Shapefile 1,5 Millionen Mal.

Um ehrlich zu sein, würde ich die ganze Zeit in PostGIS stecken und SQL für indizierte Tabellen verwenden.

Ein kurzer Blick auf Ihren Code bringt einige Optimierungen in den Sinn:

• Überprüfen Sie jeden Punkt zuerst mit dem Begrenzungsrahmen / Umschlag der Polygone, um offensichtliche Ausreißer zu beseitigen. Sie könnten noch einen Schritt weiter gehen und die Anzahl der B-Boxen zählen, in denen ein Punkt liegt. Wenn es genau einer ist, muss er nicht gegen die komplexere Geometrie getestet werden (na ja, wenn er in mehr liegt, ist es tatsächlich so.) als eins, muss es weiter getestet werden. Sie können zwei Durchgänge durchführen, um die einfachen Fälle aus den komplexen Fällen zu entfernen.)

• Anstatt jeden Punkt zu durchlaufen und gegen Polygone zu testen, durchlaufen Sie die Polygone und testen Sie jeden Punkt. Das Laden / Konvertieren von Geometrie ist langsam, daher möchten Sie so wenig wie möglich tun. Erstellen Sie außerdem zunächst eine Liste mit Punkten aus der CSV, um zu vermeiden, dass Sie dies mehrmals pro Punkt tun und die Ergebnisse am Ende dieser Iteration verwerfen müssen.

• Indizieren Sie Ihre Punkte räumlich, indem Sie sie entweder in ein Shapefile, eine SpatialLite-Datei oder in eine PostGIS / PostgreSQL- Datenbank konvertieren . Dies hat den Vorteil, dass Tools wie OGR den größten Teil der Arbeit für Sie erledigen können.

• Schreiben Sie die Ausgabe erst am Ende: print () ist im besten Fall eine teure Funktion. Speichern Sie die Daten stattdessen als Liste und schreiben Sie sie ganz am Ende mit den standardmäßigen Python-Beiz- oder List-Dumping-Funktionen aus.