Menggunakan OGR dan Shapely lebih efisien? [Tutup]

Saya mencari beberapa saran tentang cara membuat kode python saya lebih efisien. Biasanya efisiensi tidak masalah bagi saya, tetapi saya sekarang bekerja dengan file teks lokasi AS dengan lebih dari 1,5 juta poin. Dengan pengaturan yang diberikan, dibutuhkan sekitar 5 detik untuk menjalankan operasi pada satu titik; Saya perlu menurunkan angka ini.

Saya menggunakan tiga paket GIS python yang berbeda untuk melakukan beberapa operasi yang berbeda pada poin dan menghasilkan file teks baru yang dibatasi.

1. Saya menggunakan OGR untuk membaca shapefile batas county dan mendapatkan akses ke geometri batas.
2. Dengan cermat memeriksa untuk melihat apakah suatu titik ada di dalam salah satu dari kabupaten ini.
3. Jika ada di dalam satu, saya menggunakan Python Shapefile Library untuk menarik informasi atribut dari batas .dbf.
4. Saya kemudian menulis beberapa informasi dari kedua sumber ke file teks.

Saya menduga bahwa inefisiensi terletak pada memiliki loop 2-3 tingkat ... tidak yakin apa yang harus dilakukan tentang itu. Saya terutama mencari bantuan dengan seseorang yang berpengalaman dalam menggunakan salah satu dari 3 paket ini, karena ini adalah pertama kalinya saya menggunakannya.

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"

Langkah pertama adalah memindahkan shapefile terbuka di luar loop baris, Anda membuka dan menutup shapefile 1,5 juta kali.

Sejujurnya saya akan memasukkan semuanya ke PostGIS dan melakukannya menggunakan SQL pada tabel yang diindeks.

Melihat cepat pada kode Anda membawa beberapa optimasi ke pikiran:

• Periksa setiap titik terhadap kotak pembatas / amplop poligon terlebih dahulu, untuk menghilangkan outlier yang jelas. Anda bisa melangkah lebih jauh dan menghitung jumlah bbox di mana titik berada, jika tepat satu, maka itu tidak perlu diuji terhadap geometri yang lebih kompleks (baik, sebenarnya akan menjadi jika itu terletak di lebih dari satu, itu perlu diuji lebih lanjut. Anda bisa melakukan dua melewati untuk menghilangkan kasus-kasus sederhana dari kasus-kasus kompleks).

• Alih-alih mengulang melalui setiap titik dan mengujinya terhadap poligon, loop melalui poligon dan menguji setiap titik. Memuat / mengonversi geometri lambat, jadi Anda ingin melakukannya sesedikit mungkin. Juga, buat daftar Poin dari CSV pada awalnya, sekali lagi untuk menghindari keharusan melakukannya berkali-kali per poin kemudian membuang hasilnya di akhir iterasi itu.

• Secara spasial indeks poin Anda, yang melibatkan mengonversinya menjadi file shapefile, SpatialLite, atau sesuatu seperti database PostGIS / PostgreSQL . Ini memiliki keuntungan bahwa alat seperti OGR akan dapat melakukan sebagian besar pekerjaan untuk Anda.

• Jangan menulis hasilnya sampai selesai: print () adalah fungsi yang mahal pada saat terbaik. Alih-alih menyimpan data sebagai daftar, dan menuliskannya di bagian paling akhir menggunakan fungsi pengawetan Python standar atau fungsi daftar-dumping.