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
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_rasterize
in 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(...)
No comments:
Post a Comment