我正在寻找有关如何提高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"


评论

您可以考虑发布此@ Code Review:codereview.stackexchange.com

#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