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?

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.

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

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.