GDAL Rasterization – Rasterizing Shapefiles with GDAL and Python

Question

I am trying to rasterize a shapefile, and write values from a specific column of the shapefile into the resulting GTiff.

Here is what I've done so far, but that only creates a GTiff of zeros.

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

ndsm = 'base.tif'
shp = 'polygon.shp'
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
source_layer = data.GetLayer()
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
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
target_ds = gdal.GetDriverByName('GTiff').Create('my.tiff', x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
band = target_ds.GetRasterBand(1)
NoData_value = -999999
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])

How do I get the values of column 'hedgerow' into the GTiff?
In addition, I have to use the option 'ALL_TOUCHED=TRUE' with RasterizeLayer. How does the gdal.RasterizeLayer call has to look with both options in it?

Edit: using gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow", "ALL_TOUCHED=TRUE"]) runs without any problem, yet I only get zeros.

Answer 1

I tried out your code, slightly modified by me for adapting it to mi data (adding your 'hedgerow' field to my shapefile), and it ran well with only add to your code the line: target_ds = None.

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

ndsm = '/home/zeito/pyqgis_data/utah_demUTM2.tif'
shp = '/home/zeito/pyqgis_data/polygon8.shp'
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
#source_layer = data.GetLayer()
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
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
output = '/home/zeito/pyqgis_data/my.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))
band = target_ds.GetRasterBand(1)
NoData_value = -999999
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])

target_ds = None

This is the condition before running the code:

After running the code I got:

Shapefile was adequately rasterized by 'hedgerow' field.