.Length() method from ogr is crashing python

Question

I created a simple 1-feature polyline shapefile in QGIS 2.18. I would like to now read this shapefile in gdal using python and extract its length.

I have done the following

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
data_source = driver.Open('simple_line.shp')
layer = data_source.GetLayer()
line = layer.GetFeature(0).geometry()
line.Length()

The last line crashes python:

python.exe has stopped working

A problem caused the program to stop working correctly. Please close the program.

I guess it is to do with my shapefile since when I create a geometry using more "conventional" gdal means I get the following:

import ogr

line = ogr.Geometry(ogr.wkbLineString)
line.AddPoint(1116651.439379124, 637392.6969887456)
line.AddPoint(1188804.0108498496, 652655.7409537067)
line.AddPoint(1226730.3625203592, 634155.0816022386)
line.AddPoint(1281307.30760719, 636467.6640211721)
line.Length()

The output to the above is 170573.31786365344 which implies it works.

What could be the reason that my shapefile-derived geometry object is crashing python?

Answer 1

You have been bitten by a well known GDAL/OGR Python "gotcha":

Issue: Your feature reference has gone out of scope and has orphaned your geometry object which causes a crash.

Python crashes if you use an object after deleting an object it has a relationship with

This problem can manifest itself in subtle ways. For example, can occur if you try to instantiate a temporary dataset instance within a single line of code:

>>> print gdal.Open('C:\\RandomData.img').GetRasterBand(1).Checksum()
< Python crashes >

In this example, the dataset instance was no longer needed after the call to GetRasterBand() so Python deallocated it before calling Checksum().

This problem occurs because the GDAL and OGR objects are implemented in C++ and the relationships between them are maintained in C++ using pointers. When you delete the dataset instance in Python it causes the C++ object behind it to be deallocated. But the C++ object behind the band instance does not know that this happened, so it contains a pointer to the C++ dataset object that no longer exists. When the band tries to access the non-existing object, the process crashes.

The GDAL team knows that this design is not what Python programmers expect. Unfortunately the design is difficult to correct so it is likely to remain for some time.

The problem is not restricted to the GDAL band and dataset objects. It happens in other areas where objects have relationships with each other. Unfortunately there is no complete list, so you have to watch for it yourself. One other known place involves the OGR GetGeometryRef() function.

Solution: Assign the feature to a variable before accessing the geometry. Or use fiona, a much more pythonic interface to OGR.

Answer 2

I found out how to prevent it from crashing, but I have no idea why it was crashing in the first place, since what I have posted in my original question is valid python syntax.

It turns out that you cannot set the geometry all in one go, but need to extract the feature and assign it to a variable first. The working version of the code above is:

import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
data_source = driver.Open('simple_line.shp')
layer = data_source.GetLayer()
feature = layer.GetFeature(0)
line = feature.geometry()
line.Length()

If anyone can give a thorough explanation for why this is necessary (from a pythonic point of view), I would consider that a useful improvement of this answer.