[gdal-dev] Help with python gdal.RasterizeLayer

[Rudi von Staden]

Greetings all,

I am trying to create a raster mask using a polygon template. I can 
create the raster datasource to be the right extent, and using the same 
projection, but somehow gdal.RasterizeLayer does not do anything to the 
mask. Am I missing something in the implementation?

Below I've included the output I'm getting (followed by my python 
script). The RasterizeLayer command should be burning the polygon onto 
the single band raster (initialized to have all cells = 255) using a 
value of 0. However, the stats show that nothing has happened, and if I 
check the resulting GTif, all cells are still 255.

Any suggestions welcome.


Thanks,
Rudi

========Polygon Layer========
UL:  24.9371 -33.0525
LR:  27.8373 -33.9152
Spatial Reference:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.257223563,
             AUTHORITY["EPSG","7030"]],
         TOWGS84[0,0,0,0,0,0,0],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.0174532925199433,
         AUTHORITY["EPSG","9108"]],
     AUTHORITY["EPSG","4326"]]

========Mask Raster========
UL:  24.9371 -33.0525
LR:  27.8365142765 -32.1900501578
Spatial Reference:
GEOGCS["WGS 84",
     DATUM["WGS_1984",
         SPHEROID["WGS 84",6378137,298.257223563,
             AUTHORITY["EPSG","7030"]],
         TOWGS84[0,0,0,0,0,0,0],
         AUTHORITY["EPSG","6326"]],
     PRIMEM["Greenwich",0,
         AUTHORITY["EPSG","8901"]],
     UNIT["degree",0.0174532925199433,
         AUTHORITY["EPSG","9108"]],
     AUTHORITY["EPSG","4326"]]
Stats before burn:  [255.0, 255.0, 255.0, 0.0]
Stats after burn:  [255.0, 255.0, 255.0, 0.0]


========[burn_test.py]========
import sys
from osgeo import ogr,osr,gdal
from osgeo.gdalconst import *

points = ([24.937100000000001, -33.837400000000002], [25.2897, 
-33.671399999999998], [26.1754, -33.415300000000002], 
[27.837299999999999, -33.052500000000002], [27.693300000000001, 
-33.151899999999998], [25.475999999999999, -33.915199999999999], 
[25.238299999999999, -33.880899999999997], [24.937100000000001, 
-33.837400000000002])
pixelWidth = 0.000899322046076
pixelHeight = -0.000899322046076

#create polygon
ring = ogr.Geometry(ogr.wkbLinearRing)
for point in points:
     ring.AddPoint(point[0],point[1])
poly = ogr.Geometry(ogr.wkbPolygon)
poly.AddGeometry(ring)

#create memory layer for polygon
srWGS84 = osr.SpatialReference()
srWGS84.ImportFromProj4('+proj=longlat +datum=WGS84')
polyds = ogr.GetDriverByName('Memory').CreateDataSource('')
polylyr = polyds.CreateLayer('poly',geom_type=ogr.wkbPolygon, srs=srWGS84)
feat = ogr.Feature(polylyr.GetLayerDefn())
feat.SetGeometryDirectly(poly)
polylyr.CreateFeature(feat)

#check polygon layer extent
extent = polylyr.GetExtent()
print '========Polygon Layer========'
print 'UL: ', extent[0], extent[3]
print 'LR: ', extent[1], extent[2]
print 'Spatial Reference:\n',polylyr.GetSpatialRef(),"\n"
x = extent[0]
y = extent[3]
xLR = extent[1]
yLR = extent[2]
xSpan = int((xLR-x)/pixelWidth)
ySpan = int((yLR-y)/pixelHeight)


#create mask from polygon
mask = gdal.GetDriverByName('MEM').Create('',xSpan,ySpan,1,GDT_Byte)
transform = [x,pixelWidth,0.0,y,0.0,-pixelHeight]
mask.SetGeoTransform(transform)
mask.SetProjection(srWGS84.ExportToWkt())
projstring = srWGS84.ExportToPrettyWkt()
mask.GetRasterBand(1).Fill(255)

#check mask extent
transform = mask.GetGeoTransform()
xMask = transform[0]          #top left of coverage
yMask = transform[3]          #top left of coverage
maskPixelWidth = transform[1]
maskPixelHeight = transform[5]
print "========Mask Raster========"
print "UL: ", xMask, yMask
print "LR: ", xMask+(xSpan*maskPixelWidth), yMask+(ySpan*maskPixelWidth)
print "Spatial Reference:\n", projstring
print "Stats before burn: ", mask.GetRasterBand(1).GetStatistics(0,1)

#burn polygon onto mask
err = gdal.RasterizeLayer(mask,[1],polylyr,burn_values=[0])
if err != 0:
     print(err)
     gdaltest.post_reason( 'got non-zero result code from RasterizeLayer' )
     sys.exit(1)
print "Stats after burn: ", mask.GetRasterBand(1).GetStatistics(0,1)
gdal.GetDriverByName('GTiff').CreateCopy('mask.tif',mask)

-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/gdal-dev/attachments/20100807/ccb87f8e/attachment-0001.html