Python GDAL/OGR: retrieved feature is NoneType

Question

I'm doing some scripts in Python to deal with shapefiles, and I have an issue I don't understand. I'm looping over a shapefile, and I want to retrieve the geometry of each feature, but very often I get an error saying:

    pointGeo = pointFeat.GetGeometryRef()
AttributeError: 'NoneType' object has no attribute 'GetGeometryRef'

Here's the way I did it:

        for j in range(pointsLayer.GetFeatureCount()):
                pointFeat = pointsLayer.GetNextFeature()
                pointGeo = pointFeat.GetGeometryRef()
                if pointGeo.Within(whatever):
                    ...

I don't think my shapefile is corrupted or damaged, as I tried with several... and I often get this error, but for some features I don't get it.
Any idea how to solve this?

Thanks

Answer 1

It is a pure Python problem and not an GDAL/OGR problem. But if you want to program in Python,you need to understand the difference between an iterator and an iterable (many, many references on the Web, but one of the best is Loop Like A Native by Ned Batchelder). These are fundamental concepts in Python.

"The important operation on an iterable is iter(), which will return an iterator. And the only operation available on an iterator is next(), which either returns the next value, or raises StopIteration, a special exception that means the iteration is finished." (from Ned Batchelder)

An iterable produces an iterator and an iterator produces a stream of values.

So what's the difference between the for loop and pointsLayer.GetNextFeature() ?

the for loops = iterator and pointsLayer= iterable

for feature in iterable:
   statements

so with your example:

from osgeo import ogr
source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
for feature in pointsLayer:
    geom =feature.GetGeometryRef()
    xy = geom.GetPoint()
    print xy

(272022.68669955182, 155404.12013724342, 0.0)
(272904.99338241993, 152881.6706538822, 0.0)
.....
(272796.14718989708, 152075.00336062044, 0.0)

But we can create another type of iterator with iter() and next():

source = ogr.Open('yourpointshapefile.shp')
pointsLayer = source.GetLayer()
iterator = iter(pointsLayer)
# first feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272022.68669955182, 155404.12013724342, 0.0)
# second feature in pointsLayer
feature = iterator.next()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
print xy
(272904.99338241993, 152881.6706538822, 0.0)
....
# end raises StopIteration error to signal that iteration is complete
feature = iterator.next()
Traceback (most recent call last):
...
StopIteration

Unlike the for loop that handles the StopIteration error, you need here another control structure (while loop,if,...)

Thus, what is pointsLayer.GetNextFeature() ?

It is an iterator and you can replace iterator = iter(pointsLayer) and iterator.next() by

feature = pointsLayers.GetNextFeature()
geom = feature.GetGeometryRef()
xy = geom.GetPoint()
.....
feature = pointsLayers.GetNextFeature()
....

I hope that I have been able to explain why you should not mix the for loops (iterator) with .GetNextFeature() (another iterator).

Answer 2

In general GDAL/OGR don't throw errors in Python when they can't retrieve a feature object. Instead they return None objects. This means that you have to manually test for a feature existing when you retrieve it.

Secondly, if you're using GetNextFeature, OGR reads only in one direction, like if you were reading a text file, so to go back to the beginning of the layer you need to use lyr.ResetReading() as otherwise you will be attempting to GetNextFeature from the end of file.

Happily you can do this a bit more easily as in newer versions of GDAL/OGR, OGR datasets and layes are indexable (so you can loop over them like lists). In your case you could try:

dataset = ogr.Open('some.shp')

for points_layer in dataset:
    for point_feat in points_layer:
        if point_feat: #check the existence of the point manually
            point_geo = point_feat.GetGeometryRef()
            if point_geo.Within(whatever):
                ...

Using the list interface method should avoid needing to check for features, but can be helpful if data might be a bit dodgey.