Rasterize polygons using gdal/ogr python API

Question

My understanding so far is that 'gdal' is used for raster data and 'ogr' for vector data. Given this I have been trying to rasterize a polygon data set (output raster metadata is given by a raster template). I have tried with the following piece of code, but I am afraid that there are quite a few things that are probably wrong.

from osgeo import gdal
from osgeo import ogr

#Geodatabase and Polygon Feature Class name
path_gdb = r'Z:\......\Flood_Defences.gdb'
fc = r'FloodDefences'

#Reading polygons with ogr
driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(path_gdb, 0)
flood_def = gdb.GetLayer(fc)
sr = flood_def.GetSpatialRef()

#Path to the tif file where the rasterized polygons will be stored
out_path = r'Z:\.....\raster.tif'


#Template raster and its information
template_path = r'Z:\.......template.tif'
raster = gdal.Open(template_path)
gt = raster.GetGeoTransform()
ndv = raster.GetRasterBand(1).GetNoDataValue()
proj = raster.GetProjection()
driver = raster.GetDriver().ShortName
cols = raster.RasterXSize
rows = raster.RasterYSize
count = raster.RasterCount
dtype = gdal.GetDataTypeName(raster.GetRasterBand(1).DataType)


#Rasterize polygons
gdal.Rasterize(out_path, gdb, format=driver, outputType=gdal.gdalconst.GDT_Float32,                                      
               creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256'],  
               outputSRS='EPSG:29902', layers=fc, width=cols, height=rows, xRes=gt[1], yRes=gt[5],
               noData=ndv, targetAlignedPixels=True, allTouched='YES', burnValues=1)

I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-177-7fcf99466a93> in <module>
      3                outputSRS='EPSG:29902', layers=fc,
      4                width=cols, height=rows, xRes=gt[1], yRes=gt[5], noData=ndv,
----> 5                targetAlignedPixels=True, allTouched='YES', burnValues=1)

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in Rasterize(destNameOrDestDS, srcDS, **kwargs)
   1137 
   1138     if _is_str_or_unicode(destNameOrDestDS):
-> 1139         return wrapper_GDALRasterizeDestName(destNameOrDestDS, srcDS, opts, callback, callback_data)
   1140     else:
   1141         return wrapper_GDALRasterizeDestDS(destNameOrDestDS, srcDS, opts, callback, callback_data)

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in wrapper_GDALRasterizeDestName(*args)
   3488 def wrapper_GDALRasterizeDestName(*args):
   3489     """wrapper_GDALRasterizeDestName(char const * dest, Dataset srcDS, GDALRasterizeOptions options, GDALProgressFunc callback=0, void * callback_data=None) -> Dataset"""
-> 3490     return _gdal.wrapper_GDALRasterizeDestName(*args)
   3491 class GDALBuildVRTOptions(_object):
   3492     """Proxy of C++ GDALBuildVRTOptions class."""

TypeError: in method 'wrapper_GDALRasterizeDestName', argument 2 of type 'GDALDatasetShadow *'

Also, neither sr nor proj seems to give the information I need about the reference system. In the end I had to hard-code it in the function call ('EPSG:29902').

Answer 1

The error refers to the Type of the object you passed to GetDataTypeName, the line in question is:

dtype = gdal.GetDataTypeName(raster.GetRasterBand(1).DataType)

Fot the Rasterize function and its arguments, you can see the help for guessing which Types the function expects:

In [2]: gdal.Rasterize?
Signature: gdal.Rasterize(destNameOrDestDS, srcDS, **kwargs)
Docstring:
Burns vector geometries into a raster
Arguments are :
  destNameOrDestDS --- Output dataset name or object
  srcDS --- a Dataset object or a filename
Keyword arguments are :
  options --- return of gdal.RasterizeOptions(), string or array of strings
  other keywords arguments of gdal.RasterizeOptions()
If options is provided as a gdal.RasterizeOptions() object, other keywords are ignored.
File:      /usr/lib/python3.8/site-packages/osgeo/gdal.py
Type:      function

srource and destination datasets can be passed as objects of filenames, so you can use gdb or path_gdb for source.

Other options are passed to gda.RasterizeOptions:

In [3]: gdal.RasterizeOptions?
Signature:
gdal.RasterizeOptions(
    options=None,
    format=None,
    outputType=0,
    creationOptions=None,

...

  format --- output format ("GTiff", etc...)
  outputType --- output type (gdalconst.GDT_Byte, etc...)

...

  bands --- list of output bands to burn values into

outputType refers to the type of the output pixels, i.e. for float32 use gdal.gdalconst.GDT_Float32

bands take a list of output bands. (not a count)

for other parameters, run gdal.RasterizeOptions? on ipython or jupyter and read the docs.