How to rasterize individual feature (polygon) from shapefile using GDAL/OGR in Python?

Question

I would like to convert the individual records (polygons/features) of a shapefile to raster using GDAL/OGR. Here is the example shape and the mask raster. Rasterizing the full layer with all records (two, in this example) at once works:

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

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, 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
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()

# convert the layer with field 'presence' to raster
target_ds = gdal.GetDriverByName('GTiff').Create(r"C:\Temp\test2_rasterized.tif", x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
gdal.RasterizeLayer(target_ds, [1], layer, options=["ATTRIBUTE=presence"])

When I loop over the individual features, I can't find a way to export the feature to raster:

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

raster_file = r"c:\Temp\mask.tif"
shape_file = r"c:\Temp\test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, 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
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    output = os.path.join(r"C:\Temp", "rasterized"+str(feat_id)+".tif")
    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

    # this does not work since I do not pass a layer to the command
    gdal.RasterizeLayer(target_ds, [1], feature, options=["ATTRIBUTE=presence"])

    # e.g. is it possible to achieve it with Rasterize?
    gdal.Rasterize(...) 

Is what I am after at all possible using GDAL/OGR? Or do I need a workaround like converting each feature to an individual layer and then rasterize each layer?

Answer 1

Complete workaround for converting each feature to an individual layer and then rasterize each one it can be observed in following code (where my own paths to layers were used):

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

raster_file = "/home/zeito/pyqgis_data/mask.tif"
shape_file = "/home/zeito/pyqgis_data/test2.shp"

# open raster and obtain extent and properties
data = gdal.Open(raster_file, 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
pixel_width = geo_transform[1]

# open shapefile and create layer
driver = ogr.GetDriverByName("ESRI Shapefile")
dataSource = driver.Open(shape_file, 0)
layer = dataSource.GetLayer()
nr_features = layer.GetFeatureCount()

for feat_id in range(nr_features):
    feature = layer.GetFeature(feat_id)
    id_no = feature.GetField("id_no")
    presence = feature.GetField("presence")
    outShapefile = os.path.join("/home/zeito/pyqgis_data/", "feature" + str(feat_id)+".shp")
    output = os.path.join("/home/zeito/pyqgis_data/", "rasterized" + str(feat_id)+".tif")

    outDriver = ogr.GetDriverByName("ESRI Shapefile")

    # Remove output shapefile if it already exists
    if os.path.exists(outShapefile):
        outDriver.DeleteDataSource(outShapefile)

    # Create the output shapefile
    outDataSource = outDriver.CreateDataSource(outShapefile)
    outLayer = outDataSource.CreateLayer("feature" + str(feat_id), geom_type=ogr.wkbPolygon)

    # Add an ID field
    idField = ogr.FieldDefn("id_no", ogr.OFTInteger)
    outLayer.CreateField(idField)

    # Add a presence field
    presField = ogr.FieldDefn("presence", ogr.OFTInteger)
    outLayer.CreateField(presField)

    # Create the feature and set values
    featureDefn = outLayer.GetLayerDefn()
    feat = ogr.Feature(featureDefn)
    geom = feature.GetGeometryRef()

    feat.SetGeometry(geom)
    feat.SetField("id_no", id_no)
    feat.SetField("presence", presence)
    outLayer.CreateFeature(feat)
    feat = None

    # Save and close DataSource
    outDataSource = None

    target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
    target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))

    new_dataSource = driver.Open(outShapefile, 0)
    new_layer = new_dataSource.GetLayer()

    # this works now 
    gdal.RasterizeLayer(target_ds, [1], new_layer, options=["ATTRIBUTE=presence"])
    target_ds = None

    new_dataSource = None

dataSource = None

After running above code at Python Console of QGIS 3.4, vector and raster layers were obtained as expected.

Initial situation:

Rasterized feature 0:

Rasterized feature 1:

Answer 2

Not sure if anyone is struggling with this but I spent a day on this and the solution was simple. I was itereating through features in a shapefile to export images of a specific feature with an empty mask outside each features geometry. Easier than I thought but it took me a while.

feature = layer.GetFeature(i)
geom = feature.GetGeometryRef()
# gets extents to create image within
minX, maxX, minY, maxY = geom.GetEnvelope() 
# convert geometry to json 
json_coords = geom.ExportToJson() coordinate

# use the json coordinates as culineDSName 
# and cropToCutline as True to create the mask with GdalWarp
output_options = gdal.WarpOptions(outputBounds=[minX, minY, maxX, maxY], 
                                  dstAlpha=True,  
                                  multithread=True, 
                                  cutlineDSName = json_coords,
                                  cropToCutline=True)
gdal.Warp(output_file, Raster, options=output_options)