Add field to a geojson and populate it - python

Question

Took from Python GDAL/OGR Cookbook 1.0 documentation, I insert a "in memory" geojson in string format and a gdal "in memory" raster as input_zone_polygon in loop_zonal_stats, then I receive a int value from zonal_stats() for each feature. This value is stored in a dictionary by FID.

Original code:

def loop_zonal_stats(input_zone_polygon, input_value_raster):

    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()
    featList = range(lyr.GetFeatureCount())
    statDict = {}

    for FID in featList:
        feat = lyr.GetFeature(FID)
        maxValue = zonal_stats(feat, input_zone_polygon, input_value_raster)
        statDict[FID] = int(maxValue)
    return statDict

print statDict 
{0: 5, 1: 5, 2: 5, 3: 5, ...}

The objective is to add a field named Densitykm2 and populate it with the value acquired from zonal_stats() directly in the geojson "in memory". This updated geojson will be the return of the method.

So here is the modified code:

def loop_zonal_stats(input_zone_polygon, input_value_raster):

    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()
    featList = range(lyr.GetFeatureCount())
    maxFldDef = ogr.FieldDefn('Densitykm2', ogr.OFTInteger)
    lyr.CreateField(maxFldDef)

    for FID in featList:
        feat = lyr.GetFeature(FID)
        maxValue = zonal_stats(feat, input_zone_polygon, input_value_raster)
        feat.SetField('Densitykm2', int(maxValue))
        lyr.SetFeature(feat)
    return input_zone_polygon

What ever, I tried I received two kinds of error :

1. ERROR 6: Update from inline definition not supported
2. A list ERROR 1: Invalid index : -1

Maybe there is another possibility with some sort of merging the "in memory" geojson with returned dictionary of the original code to Add field and populate it. Maybe with geojson lib?

-- UPDATE --

Here is a solution, I code:

def loop_zonal_stats(input_zone_polygon, input_value_raster):
    # Load Polygon
    shp = ogr.Open(input_zone_polygon)
    lyr = shp.GetLayer()
    # Density field Definition
    maxFldDef = ogr.FieldDefn('Densitykm2', ogr.OFTInteger)
    # Create Temporary memory layer
    outdriver = ogr.GetDriverByName('MEMORY')
    source = outdriver.CreateDataSource('memData')
    outdriver = outdriver.Open('MEM', 1)
    # Copy of Polygon layer
    source.CopyLayer(lyr, 'temp', ['OVERWRITE=YES'])
    # Load temp Temporary memory layer
    memlyr = source.GetLayer()
    # Density Field creation
    memlyr.CreateField(maxFldDef)
    # List of ID to enumerate on
    featList = range(lyr.GetFeatureCount())
    print "Number of Wetlands in process  = {0}".format(len(featList))
    # input SpatialReference
    inSpatialRef = osr.SpatialReference()
    inSpatialRef.ImportFromEPSG(2959)
    # input SpatialReference
    outSpatialRef = osr.SpatialReference()
    outSpatialRef.ImportFromEPSG(4326)
    # create the CoordinateTransformation
    coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)
    # Empty list to be filled with Geojon Feature Objet = []
    GeojonFeatObj = []
    # Loop over features
    for FID in featList:
        FeatIn = lyr.GetFeature(FID)
        # Get feature of Polygon
        featOut = memlyr.GetFeature(FID)
        # get ZonalStats
        maxValue = zonal_stats(FeatIn, input_zone_polygon, input_value_raster)
        # Add density to feature
        featOut.SetField('Densitykm2', maxValue)
        print "Wetlands features proceed FID = {0} / {1}".format(FID, len(featList))
        # write the feature to the layer
        memlyr.SetFeature(featOut)
        # Get geometry
        geom = featOut.GetGeometryRef()
        # Transform geometry (UTM18N ---> WGS84)
        geom.Transform(coordTrans)
        # update feature geometry
        featOut.SetGeometry(geom)
        # Export the feature to a feature geojson
        GeojonFeatObj.append(featOut.ExportToJson())

    return GeojonFeatObj

Note that you need : - build a FeatureCollection for the list of feature object - Add crs as prefix - Clean ORG string mess

fc = geojson.FeatureCollection(tupple)
    outputGeojson = geojson.dumps(fc)
    #outputGeojson.replace(,'') }, {"geometry":
    prefix = '{"crs": {"type": "link", "properties": {"href": "http://spatialreference.org/ref/epsg/4326/", "type": "proj4"}}, '
    output = prefix+outputGeojson[1:]
    output = output.replace('\\', '').replace('}"]', '}]').replace('["{', '[{').replace('", "{"geometry":', ', {"geometry":')

Alex