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?
Best Answer
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):
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: