Segmentation fault (core dumped) ogr in Python

[skolios]

I am creating a script which provides the Spatial Difference of 2 Polygon Layers (shapefiles, i.e. Test and Reference Layers), when 2 feature geometries match. When there are no overlaying geometries (no matching) between the two Layer's Geometries, the script adds to the output the Geometry of the Test Feature. When running in linux, and no Geometries are matched, I receive a "Segmentation fault (core dumped)", and even after some research I haven't managed to solve the problem.

The version of GDAL is 1.9.1

Please find the code below:

import sys
import os
import osgeo.ogr as ogr


def copyFields(fromLayer, toLayer):
    """Function that Copies field names and types from one layer to another"""

    featureDefn = fromLayer.GetLayerDefn()
    for i in range(featureDefn.GetFieldCount()):
        toLayer.CreateField(featureDefn.GetFieldDefn(i))


# Setting the File Format Driver
driver = ogr.GetDriverByName('ESRI Shapefile')

Path = os.getcwd()
print "Current Path: ", Path, "\n"

# Opening Test and Reference Data Sources
TestFile = "../Scratch/AND_2012_09_BuiltupAreas.shp"
RefFile = "../Scratch/AND_2012_06_BuiltupAreas.shp"
OutFileName = "../Scratch/BuiltupAreas_2012_09_Added.shp"


# Checking if the Result File exists and (re)creating the Output Shapefile
if os.path.exists(OutFileName):
    driver.DeleteDataSource(OutFileName)

OutFile = driver.CreateDataSource(OutFileName)


# Opening Files and Testing if the files exist
g1 = ogr.Open(TestFile)
g2 = ogr.Open(RefFile)

if g1 is None:
    print "Could not open the Test File"
    sys.exit(1)

if g2 is None:
    print "Could not open the Reference File"
    sys.exit(1)

if OutFile is None:
    print "Could not create Output File"
    sys.exit(1)


# Opening Test,Reference and Output layers
TestLayer = g1.GetLayer()
RefLayer = g2.GetLayer()

# Get the number of features in the layer
TestFeatureCount = TestLayer.GetFeatureCount()
RefFeatureCount = RefLayer.GetFeatureCount()
print 'Count of Test Layer Features: ', TestFeatureCount, "\n"
print 'Count of Reference Layer Features: ', RefFeatureCount, "\n"

FtrDefn = TestLayer.GetLayerDefn()
FldCnt = FtrDefn.GetFieldCount()
print "Number of Fields in TestLayer: ", FldCnt, "\n"


# Creating Output layer and setting Geometry type same as Test Layer
OutputLayerName = os.path.splitext(os.path.split(OutFileName)[1])[0]
OutputLayer = OutFile.CreateLayer(OutputLayerName, geom_type=FtrDefn.GetGeomType())
OutputLayerDefn = OutputLayer.GetLayerDefn()


# Creating Field Definitions
copyFields(TestLayer, OutputLayer)
AreaFldDefn = ogr.FieldDefn('Area',ogr.OFTReal)
OutputLayer.CreateField(AreaFldDefn)

##for i in range(10,11):
for i in range(0,19):

    TestFeature = TestLayer.GetFeature(i)
    TestFID = TestFeature.GetFID()
    TestAreaType = TestFeature.GetField('Type')
    TestGeometry = TestFeature.GetGeometryRef()

    RefFID = None
    Match = False
    for j in range(RefFeatureCount):
        RefFeature = RefLayer.GetFeature(j)
        Overlaps = TestGeometry.Distance(RefFeature.GetGeometryRef())==0
##        print j, RefFeature.GetFID()
        if Overlaps==True:
            Match = True
            RefFID,J = RefFeature.GetFID(),j

    print "i:", i
    print "Geometric Match Found: ", Match
    
    if Match == True:
        RefFeature = RefLayer.GetFeature(J)
        OutputGeometry = TestGeometry.Difference(RefFeature.GetGeometryRef())
        OutputArea = OutputGeometry.GetArea()
        print RefFID, J
    elif Match == False:
        OutputGeometry = TestGeometry
        OutputArea = OutputGeometry.GetArea()
    
    featureDefn = OutputLayer.GetLayerDefn()
    OutputFeature = ogr.Feature(featureDefn)
    OutputFeature.SetGeometry(OutputGeometry)
    OutputFeature.SetFID(TestFeature.GetFID())
    OutputFeature.SetField('Type',TestAreaType)
    print "Field Added"
    OutputFeature.SetField('Area',OutputArea)
    print "Field Added"
    
    # destroy the geometry and feature
    OutputLayer.CreateFeature(OutputFeature)
    OutputGeometry.Destroy()
    OutputFeature.Destroy()

# Close the data source
OutFile.Destroy()

Additional information/data available upon request

[skolios]

Input Files

[Even Rouault]

You must remove OutputGeometry.Destroy(). In the non matching case, OutputGeometry = TestGeometry, and TestGeometry = TestFeature.GetGeometryRef(). But geometries returned by GetGeometryRef() should never be destroyed explicitely, since they remain under the ownership of their belonging features. Perhaps something to add in http://trac.osgeo.org/gdal/wiki/PythonGotchas

More generally, assign None if you want to force deletion of the underlying C++ object, and avoid using .Destroy() which is an historical remain of the "old generation" python bindings.