Rasterizing shapefiles with GDAL and Python?
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:
After running the code I got:
Shapefile was adequately rasterized by 'hedgerow' field.
It's my first time writing here, but hope I might give you some insight. I used the code from your post, with some minor changes it worked… kind of. I checked the correlation matrices that I get by comparing some of my classification data with data rasterized in the QGIS tool and raster generated by this code. Raster generated by your code had kappa coefficient near 0.4, while that generated by QGIS was near 0.9. I think the problem lies somewhere in
target_ds.SetGeoTransform((x_min, pixel_width, 0, y_min, 0, pixel_width))
Using .GetGeoTransform() solved the problem for me (got same correlation matrices as from QGIS data)
target_ds.SetGeoTransform(data.GetGeTransform())