[GIS] Rasterizing Shapefile with attribute value as pixel value with GDAL in Python

gdalpythonrastershapefile

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

Best Answer

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.

Related Question