5

I am trying to rasterize a shapefile, and write values from a specific column of the shapefile into the resulting GTiff.

Here is what I've done so far, but that only creates a GTiff of zeros.

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

ndsm = 'base.tif'
shp = 'polygon.shp'
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
target_ds = gdal.GetDriverByName('GTiff').Create('my.tiff', x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
band = target_ds.GetRasterBand(1)
NoData_value = -999999
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])

How do I get the values of column 'hedgerow' into the GTiff? In addition, I have to use the option 'ALL_TOUCHED=TRUE' with RasterizeLayer. How does the gdal.RasterizeLayer call has to look with both options in it?

Edit: using gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow", "ALL_TOUCHED=TRUE"]) runs without any problem, yet I only get zeros.

  • I think that you only need this line target_ds = None at the botton of your code (see my answer). – xunilk Oct 3 '16 at 11:11
  • What will target_ds = None do? – Thomas Becker Oct 3 '16 at 11:46
  • Close the flow of data. In my case, if I don't use that instruction I get a GTiff of zeros. – xunilk Oct 3 '16 at 11:49
  • By the way, you should use absolute paths to shapefile and raster. The problem can be there. – xunilk Oct 3 '16 at 11:54
  • Hmm, I changed for absolute paths. Still all zeros. – Thomas Becker Oct 3 '16 at 12:28
5

I tried out your code, slightly modified by me for adapting it to mi data (adding your 'hedgerow' field to my shapefile), and it ran well with only add to your code the line: target_ds = None.

from osgeo import gdal
from osgeo import ogr
from osgeo import gdalconst

ndsm = '/home/zeito/pyqgis_data/utah_demUTM2.tif'
shp = '/home/zeito/pyqgis_data/polygon8.shp'
data = gdal.Open(ndsm, gdalconst.GA_ReadOnly)
geo_transform = data.GetGeoTransform()
#source_layer = data.GetLayer()
x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * data.RasterXSize
y_min = y_max + geo_transform[5] * data.RasterYSize
x_res = data.RasterXSize
y_res = data.RasterYSize
mb_v = ogr.Open(shp)
mb_l = mb_v.GetLayer()
pixel_width = geo_transform[1]
output = '/home/zeito/pyqgis_data/my.tif'
target_ds = gdal.GetDriverByName('GTiff').Create(output, x_res, y_res, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
band = target_ds.GetRasterBand(1)
NoData_value = -999999
band.SetNoDataValue(NoData_value)
band.FlushCache()
gdal.RasterizeLayer(target_ds, [1], mb_l, options=["ATTRIBUTE=hedgerow"])

target_ds = None

This is the condition before running the code:

enter image description here

After running the code I got:

enter image description here

Shapefile was adequately rasterized by 'hedgerow' field.

  • Sorry, adding the target_ds = None doesn't change my result. I'm not sure if there can be a problem with the extent or not. The shapefile has a much bigger spatial extent than the GTiff where I retrieve the extent from. However, the GTiff I create using gdal.RasterizeLayer has the right extent. I also checked the layer retrieved from the shapefile, after opening it, and there are more than 300.000 objects in it. I also checked that the column 'hedgerow' is actually there, and it is. Is there a specific format the 'hedgerow' column has to have? – Thomas Becker Oct 3 '16 at 11:44
  • I saw that your raster data type is byte so, my 'hedgerow' column was integer. – xunilk Oct 3 '16 at 11:57

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.