Utiliser OGR et Shapely plus efficacement? [fermé]

Je cherche des suggestions sur la façon de rendre mon code python plus efficace. Normalement, l'efficacité n'a pas d'importance pour moi, mais je travaille maintenant avec un fichier texte de sites américains avec plus de 1,5 million de points. Avec la configuration donnée, il faut environ 5 secondes pour exécuter des opérations sur un point; J'ai besoin de descendre ce chiffre.

J'utilise trois packages SIG python différents pour effectuer quelques opérations différentes sur les points et générer un nouveau fichier texte délimité.

1. J'utilise OGR pour lire un fichier de formes de limites de comté et accéder à la géométrie des limites.
2. Vérifie en forme pour voir si un point se trouve dans l'un de ces comtés.
3. Si elle en fait partie, j'utilise la bibliothèque de fichiers de formes Python pour extraire les informations d'attribut de la frontière .dbf.
4. J'écris ensuite des informations des deux sources dans un fichier texte.

Je soupçonne que l'inefficacité réside dans le fait d'avoir une boucle à deux ou trois niveaux ... je ne sais pas trop quoi faire à ce sujet. Je recherche en particulier de l'aide avec une personne expérimentée dans l'utilisation de l'un de ces 3 packages, car c'est la première fois que j'utilise l'un d'entre eux.

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"

La première étape serait de déplacer le fichier de formes ouvert en dehors de la boucle de lignes, vous ouvrez et fermez le fichier de formes 1,5 million de fois.

Pour être honnête, je fourrerais tout dans PostGIS et le ferais en utilisant SQL sur des tables indexées.

Un rapide coup d'œil à votre code vous rappelle quelques optimisations:

• Vérifiez d'abord chaque point par rapport à la boîte / enveloppe englobante des polygones, pour éliminer les valeurs aberrantes évidentes. Vous pouvez aller plus loin et compter le nombre de bboxes dans lesquelles se trouve un point, si c'est exactement un, alors il n'a pas besoin d'être testé par rapport à la géométrie plus complexe (enfin, ce le sera en fait s'il se trouve dans plus plus d'un, il faudra le tester plus loin. Vous pouvez faire deux passes pour éliminer les cas simples des cas complexes).

• Au lieu de parcourir chaque point et de le tester contre des polygones, parcourez les polygones et testez chaque point. Le chargement / conversion de la géométrie est lent, vous devez donc le faire le moins possible. De plus, créez initialement une liste de points à partir du CSV, encore une fois pour éviter d'avoir à le faire plusieurs fois par point, puis à ignorer les résultats à la fin de cette itération.

• Indexez spatialement vos points, ce qui implique de les convertir en un fichier de formes, un fichier SpatialLite ou quelque chose comme une base de données PostGIS / PostgreSQL . Cela a l'avantage que des outils comme OGR pourront faire la plupart du travail pour vous.

• N'écrivez la sortie qu'à la fin: print () est une fonction coûteuse dans le meilleur des cas. Au lieu de cela, stockez les données sous forme de liste et écrivez-les à la toute fin à l'aide des fonctions de décapage Python standard ou des fonctions de vidage de liste.