Problem statement
I'm using python 2.7 and ogr 1.11.0 on Windows 7.
I have two shapefiles: one for streets and one for junctions. The streets shapefile has fields listing the IDs of the start and end junctions, as well as additional information about the streets.
I need to loop over each junction, and collect information about all of the streets that are bounded by it. I started off with a spatial filter, like this:
for junction in junctions:
streets.SetSpatialFilter(junction.GetGeometryRef())
connected_streets = [street for street in streets]
But this was slow (about 1 second per junction, and I've got 191,000 junctions). So I tried to take advantage of the SRT_JUNC and END_JUNC fields in the roads shapefile, to see if this would work faster.
for junction in junctions:
junctionID = junction.GetField("ID")
connected_streets = [street for street in streets if street.GetField("SRT_JUNC") == junctionID or street.GetField("END_JUNC") == junctionID]
The list comprehension works okay in the first iteration of the loop, but for all of the following iterations it returns nothing (or a list with length zero?).
The only way I could get it to work is to Open the driver and GetLayer inside the loop. But I've got a large shapefile and the performance of this method is also too slow.
I show some simplified examples below.
Can I make multiple list comprehensions work for an ogr layer, without reloading the entire layer each time? If not, why not? (I would also accept suggestions for an entirely different approach to the main problem, but I still am looking for an explanation regarding list comprehension.)
Example 1: List comprehension without ogr layer
The following basic script
my_list = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
for x in range(3):
print len([a for a in my_list])
produces the expected result:
10
10
10
Example 2: List comprehension with ogr layer
The following script is parallel to Example 1, except it uses an ogr layer instead of a simple list.
from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(file, 0)
streets = dataSource.GetLayer()
for x in range(3):
print len([street for street in streets])
It produces an unexpected result:
191017
0
0
Example 3: move GetLayer inside the loop
I moved the GetLayer call into the loop:
from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(file, 0)
for x in range(3):
streets = dataSource.GetLayer()
print len([street for street in streets])
But still it produces the unexpected result:
191017
0
0
Example 4: move driver.Open inside the loop
I also moved the driver Open call into the loop:
from osgeo import ogr
driver = ogr.GetDriverByName("ESRI Shapefile")
for x in range(3):
dataSource = driver.Open(file, 0)
streets = dataSource.GetLayer()
print len([street for street in streets])
And it produces the expected result:
191017
191017
191017
