2

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

4

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.

  • Yes. Thank you for your help. But I changed it and got the same result.. – yeti Dec 16 '16 at 17:42
  • I have edited the post with the correct type. – yeti Dec 16 '16 at 17:44
3

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

  • Thanks for the advice. The for loop seems to started working thanks to you, because I got an error. Explained in the post edit. – yeti Dec 19 '16 at 9:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.