1

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
  • Are you trying to return the number of features in the shapefile? If so, you should use layer.GetFeatureCount() function. That would be much more efficient. – Matej Jul 29 '14 at 20:57
  • @Matej, good point, and thanks for the tip. To be more clear: I'm not trying to return the number of features. I need to collect information (name, class, etc) about the roads connected to each junction to do some comparisons per junction. The counting in my examples above is just to show that the list comprehensions are not returning results after the first run. – andytilia Jul 29 '14 at 21:28
1

I think that when you iterate over features in a gdal layer you're using OGRLayer::GetNextFeature() so after you iterated over all features, next feature will always be empty. That explains why you get your expected result when you call driver.Open(file, 0) inside the loop.

You can use OGRLayer::ResetReading() to start from the top

On the other hand, iterating over features is always gonna be slow. Doing this in postgis is not an option? It would be faster

Hope this helps

*edit: added reference to ResetReading

  • Thanks - I'll try out ResetReading. In the mean time, for the postgis option: do you mean loading the shapefiles into postgis and running the same loop in python with postgis as the datasource? – andytilia Jul 30 '14 at 0:34
  • 1
    Well, that should be faster than shapfiles I think. But what I meant was directly on postgres, with sql. – plablo09 Jul 30 '14 at 2:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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