我正在寻找有关如何提高python代码效率的一些建议。通常,效率对我而言并不重要,但是我现在正在处理一个美国点的文本文件,该文件的得分超过150万。使用给定的设置,在一点上运行操作大约需要5秒钟。我需要弄清楚这个数字。
我正在使用三个不同的python GIS软件包对这些点进行一些不同的操作,并输出一个新的带分隔符的文本文件。
我使用OGR读取县边界shapefile并访问边界几何。
检查形状是否在这些县中的任何一个县之内。
如果在其中,我使用Python Shapefile库从边界.dbf中提取属性信息。
然后我将两个来源的一些信息都写到文本文件中。
我怀疑效率低下在于2-3层循环...不太确定该怎么做。我特别希望与有经验的人一起使用这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 楼
第一步是将shapefile打开到行循环之外,您要打开和关闭shapefile 150万次。老实说,尽管我把很多东西都塞进了PostGIS并在索引表上使用SQL来完成。
#2 楼
快速浏览一下代码,您会想到一些优化方法:首先根据多边形的边界框/包络检查每个点,以消除明显的异常值。您可以更进一步,计算一个点所在的bbox的数量,如果恰好是一个,则不需要针对更复杂的几何体进行测试(嗯,实际上它应该位于更多的几何体中)。而不是一个,将需要进一步测试。您可以进行两次遍历以从复杂的情况中消除简单的情况。
而不是遍历每个点并针对多边形进行测试,而是遍历多边形并对其进行测试点。几何体的加载/转换很慢,因此您要尽可能少地做。另外,请先从CSV创建一个点列表,再次避免每次点都要做多次,然后在该迭代结束时丢弃结果。
空间索引您的点,这涉及将其转换为shapefile,SpatialLite文件或类似PostGIS / PostgreSQL数据库的文件。这样做的好处是,像OGR这样的工具将能够为您完成大部分工作。
直到最后都不要编写输出:print()在最佳情况下是一个昂贵的函数。而是将数据存储为列表,并在最后使用标准Python酸洗功能或列表转储功能将其写出。
评论
前两个将获得丰厚回报。您还可以通过对所有内容使用ogr而不是Shapely和Shapefile来加快速度。
–傻瓜
2011年5月20日在18:22
对于任何与“ Python”和“空间索引”相关的东西,与Rtree相似的是它,因为它可以很快地找到其他形状附近的形状
– Mike T
2011年5月27日下午3:29
评论
您可以考虑发布此@ Code Review:codereview.stackexchange.com