Segmentation fault (core dumped) on SetField/SetFeature

[malapradej]

I have been trying for the past 2 days to get a python script to complete creating a shapefile of 289995 points with attributes. The points can be created but the script does not complete the attributes. The code runs correctly until the for: for j, p in enumerate(wCoords): loop, when after a while the segmentation fault happens. I have tried to add an if statement which halts the processes at certain intervals to see if I can find the location in the loop cycle when it happens. The loop cycles without fault until 1000 cycles, but not at 10,000 when it it stops without feedback, seeming like an endless loop. The program is to create tree points and then attach tree height attributes to the points. The code is below:

def save_shp(wCoords):
    print 'saving shapefile...'
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists('tree_points.shp'):
        driver.DeleteDataSource('tree_points.shp')
    ds = driver.CreateDataSource('tree_points.shp')
    layer = ds.CreateLayer('trees', geom_type=ogr.wkbPoint)
    layerDefn = layer.GetLayerDefn()
    point = ogr.Geometry(ogr.wkbPoint)
    
    for i, p in enumerate(wCoords):
        point.AddPoint(p[0],p[1])
        featureIndex = i
        feature = ogr.Feature(layerDefn)
        feature.SetGeometry(point)
        feature.SetFID(featureIndex)
        layer.CreateFeature(feature)
    
    fieldDefn = ogr.FieldDefn('tree_hts', ogr.OFTReal)
    layer.CreateField(fieldDefn)
    i = feature.GetFieldIndex('tree_hts')#???
    
    for j, p in enumerate(wCoords):
        
        feature_n = layer.GetFeature(j)
        feature_n.SetField(i, p[2])#???
        layer.SetFeature(feature_n)
    
    try:
        ds.Destroy()
    except:
        print 'still core dumping!'

I don't know enough about gdal/ogr to give you any more info than this. Please help. Jacques

[Even Rouault]

The crash is due to having a 'feature' instance still alive before calling CreateField(). The C++ code assumes that there is no alive native feature object before calling CreateField(), hence the crash if that assumption is not met. And your code was a bit subobtimal. You can both create the geometry and attribute at the same time. See the below corrected snippet :

def save_shp(wCoords):
    print 'saving shapefile...'
    driver = ogr.GetDriverByName('ESRI Shapefile')
    if os.path.exists('tree_points.shp'):
        driver.DeleteDataSource('tree_points.shp')
    ds = driver.CreateDataSource('tree_points.shp')
    layer = ds.CreateLayer('trees', geom_type=ogr.wkbPoint)
    layerDefn = layer.GetLayerDefn()
    point = ogr.Geometry(ogr.wkbPoint)

    fieldDefn = ogr.FieldDefn('tree_hts', ogr.OFTReal)
    layer.CreateField(fieldDefn)
    field_index = layerDefn.GetFieldIndex('tree_hts')

    for i, p in enumerate(wCoords):
        point.AddPoint(p[0],p[1])
        featureIndex = i
        feature = ogr.Feature(layerDefn)
        feature.SetGeometry(point)
        feature.SetFID(featureIndex)
        feature.SetField(field_index, p[2])#???
        layer.CreateFeature(feature)
        feature = None

    ds = None