<div dir="ltr">thanks, <br>That does the trick. <br>Oz. <br><br><div class="gmail_quote">On Mon, Mar 30, 2009 at 7:02 PM, Even Rouault <span dir="ltr"><<a href="mailto:even.rouault@mines-paris.org">even.rouault@mines-paris.org</a>></span> wrote:<br>
<blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">Oz,<br>
<br>
See the hereafter a correction to insert in your code that will hopefully fix<br>
your issue.<br>
<br>
Le Monday 30 March 2009 17:25:13 Oz Nahum, vous avez écrit :<br>
<div><div></div><div class="h5">> Hi,<br>
> I have encountered a weired problem with gdal_rasterize.<br>
><br>
> when I run the script from here<br>
> <a href="http://trac.osgeo.org/gdal/wiki/FAQRaster#HowcanIcreateablankrasterbasedona" target="_blank">http://trac.osgeo.org/gdal/wiki/FAQRaster#HowcanIcreateablankrasterbasedona</a><br>
>vectorfilesextentsforusewithgdal_rasterize<br>
><br>
> When I run:<br>
> #!/usr/bin/env python<br>
><br>
> from osgeo import gdal<br>
> from osgeo import osr<br>
> from osgeo import ogr<br>
> import numpy<br>
> #from numpy import ones<br>
> import stat, sys, os, string, commands<br>
><br>
> shp = 'AOI'<br>
> tiff = 'AOI.tif'<br>
> px = 1<br>
> tiff_width = 80<br>
> tiff_height = 80<br>
><br>
> # Import vector shapefile<br>
> vector = ogr.GetDriverByName('ESRI Shapefile')<br>
> src_ds = vector.Open(shp +'.shp')<br>
> src_lyr = src_ds.GetLayerByIndex(index=0)<br>
> src_extent = src_lyr.GetExtent()<br>
><br>
> # Create new raster layer with 4 bands<br>
> raster = gdal.GetDriverByName('GTiff')<br>
> dst_ds = raster.Create( tiff, tiff_width, tiff_height, 1, gdal.GDT_Byte)<br>
><br>
> # Create raster GeoTransform based on upper left corner and pixel<br>
> resolution raster_transform = [src_extent[0], px, 0.0, src_extent[3], 0.0,<br>
> -px] dst_ds.SetGeoTransform( raster_transform )<br>
><br>
> # Get projection of shapefile and assigned to raster<br>
> srs = osr.SpatialReference()<br>
> srs.ImportFromWkt(src_lyr.GetSpatialRef().__str__())<br>
> dst_ds.SetProjection( srs.ExportToWkt() )<br>
><br>
> # Create blank raster with fully opaque alpha band<br>
> zeros = numpy.zeros( (tiff_height, tiff_width), numpy.uint8 )<br>
> dst_ds.GetRasterBand(1).WriteArray( zeros )<br>
<br>
</div></div># Properly close the dataset before calling gdal_rasterize, so that the TIF<br>
# file is properly written.<br>
dst_ds = None<br>
<div class="im"><br>
><br>
> commandString = 'gdal_rasterize -burn 255 -l AOI AOI.shp AOI.tif'<br>
> commandOutput = commands.getoutput(commandString)<br>
> print commandOutput<br>
> #end of code<br>
><br>
> I get the error this error:<br>
><br>
> augmented to the gdal_rasterize command I see this error. When I run it via<br>
> a shell script, I don't see the error:<br>
> ERROR 4: `AOI.tif' not recognised as a supported file format.<br>
><br>
><br>
> #!/bin/sh<br>
> python create_empty_raster.py<br>
> gdal_rasterize -burn 255 -l AOI AOI.shp AOI.tif<br>
><br>
> Runs ok.<br>
<br>
</div>--> This runs ok as at the end of the create_empty_raster.py script, the GDAL<br>
dataset is automatically closed.<br>
<div><div></div><div class="h5"><br>
><br>
> So I am confused, why I can't call gdal_rasterize from python ?<br>
> Thanks, for any help<br>
><br>
> Oz Nahum<br>
> Universtät Tuebingen<br>
> Applied Environmental Sciences<br>
<br>
<br>
</div></div></blockquote></div><br><br clear="all"><br>-- <br><br>----<br> Imagine there's no countries<br> It isn't hard to do<br> Nothing to kill or die for<br> And no religion too<br>
Imagine all the people<br> Living life in peace <br>----<br> You all must read 'The God Delusion' <br> <a href="http://en.wikipedia.org/wiki/The_God_Delusion">http://en.wikipedia.org/wiki/The_God_Delusion</a><br>
---<br>when one person suffers from a delusion it is called insanity. When many people suffer from a delusion it is called religion."<br>Robert Pirsig, Zen and the Art of Motorcycle Maintenance<br><br>
</div>