# Using OGR and Shapely more efficiently? [closed]

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 osgeo import ogr
import shapefile

pointFile = "C:\\NSF_Stuff\\NLTK_Scripts\\Gazetteer_New\\NationalFile_20110404.txt"
shapeFolder = "C:\NSF_Stuff\NLTK_Scripts\Gazetteer_New"
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")

#opens the write file
openWriteFile = open(writeFile, "w")

#uses Python Shapefile Library to read attributes from .dbf
records = sf.records()
print "Starting loop..."

#This will loop through the points in pointFile
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:
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"
``````

## closed as too broad by PolyGeo♦Jan 8 '18 at 21:11

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer. Avoid asking multiple distinct questions at once. See the How to Ask page for help clarifying this question. If this question can be reworded to fit the rules in the help center, please edit the question.

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.

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.

• 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
• 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