Использование OGR и Shapely более эффективно? [закрыто]

Я ищу несколько предложений о том, как сделать мой код Python более эффективным. Обычно эффективность не имеет значения для меня, но сейчас я работаю с текстовым файлом из США с более чем 1,5 миллионами точек. При заданной настройке выполнение операций в одной точке занимает около 5 секунд; Мне нужно уменьшить эту цифру.

Я использую три разных пакета Python GIS, чтобы выполнить несколько разных операций над точками и вывести новый текстовый файл с разделителями.

1. Я использую OGR, чтобы прочитать шейп-файл границы округа и получить доступ к геометрии границы.
2. Shapely проверяет, находится ли точка в каком-либо из этих округов.
3. Если он в пределах одного, я использую библиотеку шейп-файлов Python, чтобы извлечь информацию об атрибутах из границы .dbf.
4. Затем я записываю некоторую информацию из обоих источников в текстовый файл.

Я подозреваю, что неэффективность заключается в наличии 2-3-уровневой петли ... не совсем уверен, что с этим делать. Я особенно ищу помощи с кем-то, имеющим опыт использования любого из этих трех пакетов, так как я впервые использую любой из них.

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"

Первым шагом будет перемещение шейп-файла за пределы цикла строк, вы открываете и закрываете шейп-файл 1,5 миллиона раз.

Если честно, я бы все это вложил в PostGIS и сделал бы это с использованием SQL для индексированных таблиц.

Беглый взгляд на ваш код заставляет задуматься:

• Сначала проверьте каждую точку на ограничивающей рамке / конверте полигонов, чтобы устранить явные выбросы. Вы можете пойти еще дальше и посчитать количество bbox, в которых лежит точка, если она ровно одна, тогда ее не нужно проверять на более сложную геометрию (ну, на самом деле это будет, если она лежит в более кроме одного, его необходимо будет дополнительно протестировать. Вы можете сделать два прохода, чтобы исключить простые случаи из сложных случаев).

• Вместо того, чтобы проходить по каждой точке и проверять ее на полигоны, проходите по полигонам и проверяйте каждую точку. Загрузка / преобразование геометрии происходит медленно, поэтому вы хотите делать это как можно меньше. Кроме того, сначала создайте список точек из CSV, опять же, чтобы избежать необходимости делать это несколько раз на точку, а затем отбрасывать результаты в конце этой итерации.

• Пространственное индексирование ваших точек, которое включает преобразование его в шейп-файл, файл SpatialLite или что-то вроде базы данных PostGIS / PostgreSQL . Это имеет то преимущество, что такие инструменты, как OGR , смогут выполнить большую часть работы за вас.

• Не пишите вывод до конца: print () - дорогая функция в лучшие времена. Вместо этого сохраните данные в виде списка и запишите их в самом конце, используя стандартные функции выбора Python или функции вывода списка.