no feature ID when polygonizing raster with python GDAL

Question

I'm trying to polygonize the following binary raster:

The output is a polygon layer (ESRI shapefile) that looks good, however, doesn't have FIDs for the features and thus I am not able to iterate through them in the later part of my script.

Is there any way to add feature IDs to the shapefile's attribute table while polygonizing the raster or afterwards?

I'm using the following code:

from osgeo import gdal, ogr, osr


raster = gdal.Open(path_raster)

band = raster.GetRasterBand(1)

# determine EPSG
sp_ref = osr.SpatialReference()
sp_ref.SetFromUserInput('EPSG:3857')

drv = ogr.GetDriverByName('ESRI Shapefile') # file type: shapefile
outfile = drv.CreateDataSource(path_vec_results) #file path
layer = outfile.CreateLayer(vector_file, geom_type = ogr.wkbMultiPolygon, srs = sp_ref) 
layer.CreateField(ogr.FieldDefn('DN', ogr.OFTInteger)) #add DN column for raster valules
layer.CreateField(ogr.FieldDefn('ID', ogr.OFTInteger)) #ID column to fill later?

# polygonize
vector_result = gdal.Polygonize(band, band, layer, 0, []) 
outfile = None

Answer 1

In case someone else is having the same issue - the following code worked for me to create an FID column for the vector file's attribute table.

Happy to hear about better / shorter solutions, always looking to improve my clumsy coding.

from osgeo import gdal, ogr, osr

# load raster and get band
raster = gdal.Open(path_raster, 1)
band = raster.GetRasterBand(1)

# determine EPSG
sp_ref = osr.SpatialReference()
sp_ref.SetFromUserInput('EPSG:3857')


drv = ogr.GetDriverByName('ESRI Shapefile') # file type: shapefile
outfile = drv.CreateDataSource(path_vec_results) #path_vec_results is path to results file including filename

layer = outfile.CreateLayer(vector_file, geom_type = ogr.wkbMultiPolygon, srs = sp_ref) #file name
layer.CreateField(ogr.FieldDefn('DN', ogr.OFTInteger)) #add DN column for raster valules


# transform raster to polygon and save
vector_poly = gdal.Polygonize(band, band, layer, 0, [])

outfile.Destroy() #important!!


# opening newly created vector file using ogr
vector_final = ogr.Open(path_vec_results)

# checking polygon file - it prints feature IDs, so there are feature IDs!
for layer in vector_final:

    for feature in layer:
        print('Feature ID', feature.GetFID())    


drv = ogr.GetDriverByName("ESRI Shapefile")
datasource = drv.Open(path_vec_results, 1) # 1 for writing and not just reading
vector_final = datasource.GetLayer()

# create FID field
id_field = ogr.FieldDefn("FID", ogr.OFTInteger)
vector_final.CreateField(id_field)


for feature in vector_final:
    fid = feature.GetFID()
    feature.SetField('FID', fid)
    vector_final.SetFeature(feature)

datasource.Destroy()