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
PolyGeo
  • 65,136
  • 29
  • 109
  • 338
Hans Skov-Petersen
  • 37
  • 1
  • 2
  • 2
    Perhaps the same error with unset band count than here http://osgeo-org.1560.x6.nabble.com/gdal-dev-problem-with-datatype-in-python-td5222518.html – 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 – user109393 Nov 20 '17 at 12:47

1 Answers1

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) .

xunilk
  • 29,891
  • 4
  • 41
  • 80