1

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()
  • Works for me with a shapefile of mine. Have you tried with different input data? Generally I highly recommend using Fiona instead of the awkward GDAL bindings. – bugmenot123 Mar 13 '17 at 9:13
  • IT seems to be the dataset. When I try a different one it does work. However, the dataset was created from another shapefile using ogr. – Tins Mar 13 '17 at 9:23
  • Is it synced, closed and not used by any other process? Maybe things get weird that way. – bugmenot123 Mar 13 '17 at 9:31
  • Are you able to recreate the class.shp file from the original shapefile? It may be that something was corrupted during processing. – Ali Mar 13 '17 at 11:07
  • 'class.shp' is created from an original dataset. When I use the original it works fine but not on the new one. 'class.shp' loads in ArcMap fine though. – Tins Mar 13 '17 at 11:31
2

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.

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))
  • When I do this 'index' turns out empty. – Tins Mar 13 '17 at 11:34
  • Better DO NOT do a catch-all except and DO output the exception message! – bugmenot123 Mar 13 '17 at 13:23
  • True, didn't find error messages from gdal, if it is all RuntimeError, then that would propably better. Changed it, but not sure... – Andreas Müller Mar 13 '17 at 14:58

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.