Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

I'm looking for some suggestions about how to make my python code more efficient. Normally efficiency doesn't matter for me but I am now working with a text file of US locations with over 1.5 million points. With the given setup it is taking about 5 seconds to run operations on one point; I need to get this figure way down.

I'm using three different python GIS packages to do a few different operations on the points and output a new delimited text file.

  1. I use OGR to read a county boundary shapefile and get access to the boundary geometry.
  2. Shapely checks to see if a point is within any of these counties.
  3. If it is within one, I use the Python Shapefile Library to pull attribute information from the boundary .dbf.
  4. I then write some information from both sources to a text file.

I suspect that the inefficiency lies in having a 2-3 tiered loop... not quite sure what to do about it. I'm particularly looking for help with someone experienced in using any of these 3 packages, as it is my first time to use any of them.

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"
share|improve this question
1  
You might consider posting this @ Code Review: codereview.stackexchange.com –  RyanDalton May 20 '11 at 17:35
add comment

2 Answers

First step would be to move the shapefile open outside the rows loop, you are opening and closing the shapefile 1.5 million times.

To be honest though I'd stuff the whole lot into PostGIS and do it using SQL on indexed tables.

share|improve this answer
add comment

A quick look at your code brings a few optimisations to mind:

  • Check each point against the bounding box/envelope of the polygons first, to eliminate obvious outliers. You could go a step further and count the number of bboxes a point lies in, if it is exactly one, then it doesn't need to be tested against the more complex geometry (well, it'll actually be if it lies in more than one, it will need to be tested further. You could do two passes to eliminate the simple cases from the complex cases).

  • Instead of looping through each point and testing it against polygons, loop through the polygons and test each point. Loading/converting of geometry is slow, so you want to do it as little as possible. Also, create a list of Points from the CSV initially, again to avoid having to do it multiple times per point then discarding the results at the end of that iteration.

  • Spatially index your points, which involves converting it to either a shapefile, SpatialLite file, or something like a PostGIS/PostgreSQL database. This has the advantage that tools like OGR will be able to do most of the work for you.

  • Don't write the output until the end: print() is an expensive function at the best of times. Instead store the data as a list, and write it out at the very end using the standard Python pickling functions or list-dumping functions.

share|improve this answer
4  
The first two will pay off big. You could also speed things up a little bit by using ogr for everything instead of Shapely and Shapefile. –  sgillies May 20 '11 at 18:22
1  
For anything "Python" and "spatial index" related, look no further than Rtree as it is very fast at finding shapes near other shapes –  Mike T May 27 '11 at 3:29
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.