Unexpected output when using gdal_rasterize?

Question

I used the following code for converting my layer into raster format but I get a tif file that is completely black or white...

from osgeo import gdal, ogr

# Define pixel_size and NoData value of new raster
pixel_size = 1
NoData_value = 100

# Filename of input OGR file
# vector_fn = '/Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp'
vector_fn = '/Users/PeterVanvoorden/Desktop/GRBgis_322/New_Shapefiles/Sbne40179.shp'

# Filename of the raster Tiff that will be created
raster_fn = '/Users/PeterVanvoorden/Desktop/test.tif'

# 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()
print x_min, x_max, y_min, y_max

# 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, 1, gdal.GDT_Byte)
target_ds.SetGeoTransform((x_min, pixel_size, 0, y_max, 0, -pixel_size))
band = target_ds.GetRasterBand(1)
band.SetNoDataValue(NoData_value)

# Rasterize
gdal.RasterizeLayer(target_ds, [1], source_layer, burn_values=[0])

STARTING QUESTION:

I am trying to rasterize vector data using the GDAL_rasterize function. (Started here: syntax error when using gdal in python console QGIS) But with all methods they told me were possible, I always get something I don't need.

Below the 3 methods I tried...

---

Directly in python IDLE:

>>> import os
>>> os.system('gdal_rasterize -tr 7999.8 9999.53 -te 170000.2 168000.27 177999.8 177999.8 -a Pixelkost /Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp /Users/PeterVanvoorden/Desktop/Test.tif')

Then I get the value: 32512 as result but it should output a tif file on my desktop...

---

When I do the same in a python file and run it through IDLE:

import os

shapefile = '/Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp'   

(shapefilefilepath, shapefilename) = os.path.split(shapefile)
(shapefileshortname, extension) = os.path.splitext(shapefilename)            

outrastername = shapefileshortname + '.tif'
outraster = '/Users/PeterVanvoorden/Desktop/' + outrastername

os.system('gdal_rasterize -tr 7999.8 9999.53 -te 170000.2 168000.27 177999.8 177999.8 -a Pixelkost' +  shapefile + ' ' + outraster) 

In this case nothing happens at all.

---

Running it in terminal (OSX):

mbp-van-peter:~ PeterVanvoorden$ gdal_rasterize -tr 7999.8 9999.53 -te 170000.2 168000.27 177999.8 177999.8 -a Pixelkost /Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp /Users/PeterVanvoorden/Desktop/Test.tif

Than I get: -bash: gdal_rasterize: command not found

---

I assume, looking at the last method, that there is a problem with my gdal package. I reinstalled it after having a few problems and installed GDAL 1.10 using my terminal just with a easy install like this: sudo easy_install /Users/PeterVanvoorden/Downloads/GDAL-1.10.0

Answer 1

If you use the terminal (OSX):

Try fixing the GDAL Path. If you use the GDAL Framework of KyngChaos, the complete Path of the command is /Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize or .../Versions/1.9/...or...

So:

1) in the Terminal:

export PATH=/Library/Frameworks/GDAL.framework/Versions/1.10/Programs

and you can use gdal_rasterizein the terminal during the session.

2) put in the /Users/you/.bash_profile file the line:

export PATH=/Library/Frameworks/GDAL.framework/Versions/1.10/Programs:$PATH

and you can use gdal_rasterize in the terminal all the time

3) in all other cases (script, through IDLE, QGIS, etc.) , you must use the complete Path:

/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize

It is therefore normal that os.system('gdal_rasterize... does not work, you need

os.system('/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize...

If you use Python:

1) os.system is deprecated, you must use the subprocess module which you can control the stdout and stderr in a flexible fashion (look at 17.1.4.3. Replacing os.system(), subprocess.check_output or Calling an external command in Python). You can use subprocess.call(), subprocess.check_call()or subprocess.Popen():

args = ['/Library/Frameworks/GDAL.framework/Versions/1.10/Programs/gdal_rasterize', '-tr', '7999.8', '9999.53', '-te', '170000.2', '168000.27', '177999.8', '177999.8' '-a Pixelkost /Users/PeterVanvoorden/Desktop/GRBgis_322/Shapefile/Adp322.shp', '/Users/PeterVanvoorden/Desktop/Test.tif']
proc = subprocess.Popen(args, stdout=subprocess.PIPE)
proc.wait() # or proc.communicate()[0] or proc.stdout.read() if you want to read the errors

2) if you use the Python osgeo.gdal module, it has the gdal.RasterizeLayer() function ( Python GDAL/OGR Cookbook,Python: Rasterizing a GDAL layer, Making Rasters From Vectors or GDAL RasterizeLayer Doesn't Burn All Polygons to Raster) and you don't need to call gdal_rasterize.

from osgeo import gdal
raster = gdal.RasterizeLayer(...)