2

Trying to figure out how to rasterize in GDAL.

I have two problems (see example code below):

1) Even creating the raster layer with the gdal.GDT_Float32 type option, the resulting layer ('s bands) all end up being integer

2) In spite all efforts regarding noData setting, cells outside the polyons end up being set to zero.

# Create the destination data source
tempSource = rasterDriver.Create(tempRasterPath + tileName + "_Building.tif", tileHdl.RasterXSize, tileHdl.RasterYSize, gdal.GDT_Float32)
tempSource.SetGeoTransform(tileGeoTransformationParams)
tempTile = tempSource.GetRasterBand(1)
tempTile.SetNoDataValue(-999)
tempTile.Fill(-999)

# Rasterize
gdal.RasterizeLayer(tempSource, [1], buildingPolys.layer, options=["ATTRIBUTE=%s" % meanZAttributeName])
tempSource = None
3
  • 2
    Perhaps the same error with unset band count than here osgeo-org.1560.x6.nabble.com/…
    – user30184
    Nov 17 '17 at 17:16
  • Your issues are because you missed band parameter in created raster. Please, see my answer.
    – xunilk
    Nov 17 '17 at 22:40
  • Brillant! Thanks a lot. It works - both interms of the floats, and the NoData. By the way, what would be your recommendation for a comprehensive documentation of the attributes and methods of the OGR/GDAL classes? Besides the pydocs, of course.... I got Erik Westra’s book, but it do not a full coverage. Web resources, the GIT etc. seems rather unexplained... Hans
    – Hans SP
    Nov 20 '17 at 12:47
5

Your issues are because you missed band parameter in created raster. So, I tried following code out:

from osgeo import gdal, ogr

tileHdl = "/home/zeito/pyqgis_data/aleatorio.tif"
tempRasterPath = "/home/zeito/pyqgis_data/"
tileName = "aleatorio"

vector_layer = "/home/zeito/pyqgis_data/buildingPolys.shp"

# open the raster layer and get its relevant properties
tileHdl = gdal.Open(raster_layer, gdal.GA_ReadOnly)

tileGeoTransformationParams = tileHdl.GetGeoTransform()
projection = tileHdl.GetProjection()

rasterDriver = gdal.GetDriverByName('GTiff')

buildingPolys_ds = ogr.Open(vector_layer)
buildingPolys = buildingPolys_ds.GetLayer()

# Create the destination data source
tempSource = rasterDriver.Create(tempRasterPath + tileName + "_Building.tif", 
                                 tileHdl.RasterXSize, 
                                 tileHdl.RasterYSize,
                                 1, #missed parameter (band)
                                 gdal.GDT_Float32)

tempSource.SetGeoTransform(tileGeoTransformationParams)
tempSource.SetProjection(projection)
tempTile = tempSource.GetRasterBand(1)
tempTile.Fill(-999)
tempTile.SetNoDataValue(-999)

gdal.RasterizeLayer(tempSource, [1], buildingPolys, options=["ATTRIBUTE=value"])
tempSource = None

with following layers:

enter image description here

and after running it, I got rasterized vector layer as expected.

enter image description here

As it can be observed at Metadata Layer Properties, nodata values were adequately set as -999 (it is also Float32 data type) .

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.