Can't retrieve geometry (Python, OGR, Rtree)

Question

ISSUE FIXED SEE: '.Simplify' causes 'null geometry' errors (Python, OGR)

Im trying to create a spatial index on one of my layers for the purpose of performing an intersect analysis. However, I get the following error:

AttributeError: 'NoneType' object has no attribute 'GetEnvelope'

It seems that when I iterate over the features and try to use '.GetGeometryRef()' no geometry is returned. When I test '.GetGeometryRef()' outside of the loop it doesn't throw me an error. What could be the issue? Here is my code:

import sys, os, rtree
from osgeo import ogr

inFile = r'C:\Data\intermediates\class.shp' 
driver = ogr.GetDriverByName('ESRI Shapefile')
inData = driver.Open(inFile, 0)
inLayer = inData.GetLayer()

feature = inLayer.GetFeature(0)
geom = feature.GetGeometryRef()
if geom is None:
    sys.exit(1)
feature = None
geom = None

index = rtree.index.Index(interleaved = False)
for fid1 in range(0, inLayer.GetFeatureCount()):
    feature = inLayer.GetFeature(fid1)
    geom = feature.GetGeometryRef()
    xmin, xmax, ymin, ymax = geom.GetEnvelope()
    index.insert(fid1, (xmin, xmax, ymin, ymax))

UPDATE

Apparently the 'class.shp' contains geometry errors according to the Check Geometry tool in ArcMap, (71 issues): 'null geometry' and ' could not find spatial index'. I created class.shp from an original shapefile only by simplifying the features. I created the layer with this code, did I do something wrong here?:

if.os.path.exists('class.shp'):
    driver.DeleteDataSource('class.shp')
outFile = driver.CreateDataSource('class.shp')
outLayer = outFile.CreateLayer('class.shp', srs, geom_type = inLayer.GetLayerDefn().GetGeomType())

selectLayerDef = selectData.GetLayerDefn()
for i in range(0, selectLayerDef.GetFieldCount()):
outLayer.CreateField(selectLayerDef.GetFieldDefn(i))

outLayerDef = outLayer.GetLayerDefn()
for i in range(0, selectData.GetFeatureCount()):
    inFeature = selectData.GetNextFeature()
    outFeature = ogr.Feature(outLayerDef)
    for i in range(0, outLayerDef.GetFieldCount()):
        outFeature.SetField(outLayerDef.GetFieldDefn(i).GetNameRef(),
                            inFeature.GetField(i)
    geom = inFeature.GetGeometryRef()
    geomSimple = geom.Simplify(10)
    outFeature.SetGeometry(geomSimple.Clone())
    outFeature.SetField('FEAT_TYPE', 'LAV')
    outLayer.CreateFeature(outFeature)
    inFeature.Destroy()
    outFeature.Destroy()

Answer 1

The issue you're running into suggests that your shapefile has some features with no geometry attached.

You can use exception handling within your script to record a feature's FID if the GetGeometryRef() method fails and then investigate those features further.

Another option, since you have access to ArcMap, is the Check Geometry tool. You can use this to scan your shapefile and report on any geometry issues found.

Assuming the tool does indeed flag some null geometries, you will then have to examine your workflow to see how they were introduced.

Answer 2

I recommend using gdal.UseExceptions() to enable errors for try-except. If an operation on a geometry throws an exception you coul catch it with a continue. This discussion may help you: GDAL and Python - Don't print GDAL Error messages or in the python cookbook https://pcjericks.github.io/py-gdalogr-cookbook/gdal_general.html

For your code, you may use it like this:

gdal.UseExceptions()

# ...

for fid1 in range(0, inLayer.GetFeatureCount()):
    try:
        feature = inLayer.GetFeature(fid1)
        geom = feature.GetGeometryRef()
        xmin, xmax, ymin, ymax = geom.GetEnvelope()
    except RuntimeError as err:
        print "Error in reading feature %d - %s"%(fid1, err)
        continue # skip this feature

    index.insert(fid1, (xmin, xmax, ymin, ymax))