rasterize_features_field_values.py

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

#This raster is the model for our output (CRS, extent)
ndsm = 'geotif_output_extent_proj.tif'
#This shapefile contains the features we want to burn
shp = 'shapefile_features_fields.shp'

data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)

geo_transform = data.GetGeoTransform()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize


vec = ogr.Open(shp)
lyr = vec.GetLayer(0)
pixel_width = geo_transform[1]

output = 'output.tif'
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_UInt32)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
band = target_ds.GetRasterBand(1)
NoData_value = 0
band.SetNoDataValue(NoData_value)
band.FlushCache()


driver = ogr.GetDriverByName('ESRI Shapefile')
for feat in lyr:
    burn_value = int(feat.GetField("FIELD_NAME"))
    datasource = driver.CreateDataSource('temp.shp')
    layer = datasource.CreateLayer('temp', lyr.GetSpatialRef(), geom_type=ogr.wkbPolygon)
    layer.CreateFeature(feat)
    gdal.RasterizeLayer(target_ds, [1], layer, burn_values=[burn_value], options=["ALL_TOUCHED=TRUE"])
    datasource.Destroy()

#This makes raster to write to disk
target_ds = None