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?