Rasterizing Shapefile with attribute value as pixel value with GDAL in Python

Question

Using suggestions and this question, the code now looks like this:

from osgeo import gdal, ogr
from sys import getsizeof
import os
# Define pixel_size and NoData value of new raster
pixel_size = 1
NoData_value = -9999
# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'
# Filename of input OGR file
vector_fn = 'Layer.shp'
# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
feature = source_layer.GetNextFeature()
# Array with distinct values of field 'Source' (Classes)
field_vals =set([feature.GetFieldAsString('Source') for feature in source_layer])
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, len(field_vals), gdal.GDT_Float32) 
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
#for each value of field_vals, make a band, on the band rasterize polygons with each value. Pixels have values of object's area 
counter = 0
for i in range(source_layer.GetFeatureCount()):
 feature = source_layer.GetFeature(i)
 if feature.GetFieldAsString('Source') in field_vals:
 band = target_ds.GetRasterBand(counter+1)
 band.SetNoDataValue(NoData_value)
 # Rasterize
 # Counter +1, because there is no 0 band in GTiff raster 
 gdal.RasterizeLayer(target_ds, [counter+1], source_layer, burn_values=(255,255,255), options = ["ATTRIBUTE=%s" % AREA])
 counter += 1
target_ds = None 

I think the for loop just started to work, because now I am getting an error:

gdal.RasterizeLayer(target_ds, [counter+1], source_layer, burn_values=(255,255,255), options = ["ATTRIBUTE=%s" % AREA])

NameError: name 'AREA' is not defined

I don't know how I have to deliver 'AREA' to method as a value of the field from attribute table. In mentioned question, _RASTERIZE_COLOR_FIELD_ was defined with value "color". I am just starting coding Python and maybe the problem is in the fact that I don't get some Python things.

Second thing, I don't understand is that I have to deliver to method _burn_values_ even if I want pixels to have Area's value.

Original question:

I want to rasterize a polygon Shapefile. Shapefile consists of polygons from 28 other polygon Shapefiles. Shapefile's atrribute table has two fields: 'Source' with the name of source Shapefile and 'Area'. My target is to generate a raster with number of bands corresponding with number of sources, 1 band for each source. On each band - polygons from each source should be rasterized with values of their areas - each pixel covered by polygon should have its area as a value.

Here is my code:

from osgeo import gdal, ogr
from sys import getsizeof
import os
# Define pixel_size and NoData value of new raster
pixel_size = 1
NoData_value = -9999
# Filename of the raster Tiff that will be created
raster_fn = 'test.tif'
# Filename of input OGR file
vector_fn = 'Layer.shp'
# Open the data source and read in the extent
source_ds = ogr.Open(vector_fn)
source_layer = source_ds.GetLayer()
x_min, x_max, y_min, y_max = source_layer.GetExtent()
feature = source_layer.GetNextFeature()
# Array with distinct values of field 'Source' (Classes)
field_vals = []
while feature:
 field_vals.append(feature.GetFieldAsString('Source'))
 feature = source_layer.GetNextFeature()
# Create the destination data source
x_res = int((x_max - x_min) / pixel_size)
y_res = int((y_max - y_min) / pixel_size)
target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, len(set(field_vals)), gdal.GDT_Float32) #without parsing field_vals to set getting number of all features from layer
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
counter = 0
for val in field_vals:
 while feature:
 if (feature.GetFieldAsString('Source') == (set(field_vals)[counter])):
 band = target_ds.GetRasterBand(counter+1)
 band.SetNoDataValue(NoData_value)
 # Rasterize
 # Counter +1, because there is no 0 band in GTiff raster 
 gdal.RasterizeLayer(target_ds, [counter+1], source_layer, burn_values, options = ["BURN_VALUE_FROM=Area"])
 feature = source_layer.GetNextFeature()
 counter += 1
target_ds = None 

But I only get black raster with 28 canals and without data. There are also no pixels with value no data. What am I doing wrong?

Attribute table of Shapefile Raster result

• radouxju: you cannot create a single multiband raster with a different pixel types. you must create different raster of each type.

Answer 1

Because your Area field contains floating point values, you need to account for that when creating the raster. As it stands, you are creating a raster that can only hold a single value:

target_ds = gdal.GetDriverByName('GTiff').Create(raster_fn, x_res, y_res, len(set(field_vals)), gdal.GDT_Byte)

Instead, change gdal.GDT_Byte to gdal.GDT_Float32.

• yeti: Yes. Thank you for your help. But I changed it and got the same result..
• yeti: I have edited the post with the correct type.

Answer 2

You can greatly reduce the relative complexity of your script with List Comprehension

For example, if you want the (unique) Source values in the shapefile:

field_vals =[feature.GetFieldAsString('Source') for feature in source_layer]
print field_values
['rubble', 'rubble', 'foundations', 'foundations']
# unique values
field_vals = set(field_vals)
print field_vals
set(['foundations', 'rubble'])
# or in one line
field_vals =set([feature.GetFieldAsString('Source') for feature in source_layer])
print field_vals
set(['foundations', 'rubble'])

And you can simply use (without while or GetNextFeature())

counter = 0
for i in range(source_layer.GetFeatureCount()):
 feature = source_layer.GetFeature(i)
 if feature.GetFieldAsString('Source') in field_vals:
 band = target_ds.GetRasterBand(counter+1)
 ....
 counter += 1

With the correction of Aaron

• yeti: Thanks for the advice. The for loop seems to started working thanks to you, because I got an error. Explained in the post edit.

Answer 3

Instead of

gdal.RasterizeLayer(target_ds, [counter+1], source_layer, burn_values=(255,255,255), options = ["ATTRIBUTE=%s" % AREA])

use

gdal.RasterizeLayer(target_ds, [counter+1], source_layer, options = ["ATTRIBUTE=Area"])

You don't want to burn any values except your Attibutes from the columns Area, so leave that. And you tell RasterizeLayer the Column name directly.