<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<meta http-equiv="content-type" content="text/html; charset=ISO-8859-1">
</head>
<body text="#000000" bgcolor="#ffffff">
Greetings all,<br>
<br>
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? <br>
<br>
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.<br>
<br>
Any suggestions welcome.<br>
<br>
<br>
Thanks,<br>
Rudi<br>
<br>
========Polygon Layer========<br>
UL: 24.9371 -33.0525<br>
LR: 27.8373 -33.9152<br>
Spatial Reference:<br>
GEOGCS["WGS 84",<br>
DATUM["WGS_1984",<br>
SPHEROID["WGS 84",6378137,298.257223563,<br>
AUTHORITY["EPSG","7030"]],<br>
TOWGS84[0,0,0,0,0,0,0],<br>
AUTHORITY["EPSG","6326"]],<br>
PRIMEM["Greenwich",0,<br>
AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.0174532925199433,<br>
AUTHORITY["EPSG","9108"]],<br>
AUTHORITY["EPSG","4326"]] <br>
<br>
========Mask Raster========<br>
UL: 24.9371 -33.0525<br>
LR: 27.8365142765 -32.1900501578<br>
Spatial Reference:<br>
GEOGCS["WGS 84",<br>
DATUM["WGS_1984",<br>
SPHEROID["WGS 84",6378137,298.257223563,<br>
AUTHORITY["EPSG","7030"]],<br>
TOWGS84[0,0,0,0,0,0,0],<br>
AUTHORITY["EPSG","6326"]],<br>
PRIMEM["Greenwich",0,<br>
AUTHORITY["EPSG","8901"]],<br>
UNIT["degree",0.0174532925199433,<br>
AUTHORITY["EPSG","9108"]],<br>
AUTHORITY["EPSG","4326"]]<br>
Stats before burn: [255.0, 255.0, 255.0, 0.0]<br>
Stats after burn: [255.0, 255.0, 255.0, 0.0]<br>
<br>
<br>
<big><tt>========[burn_test.py]========<br>
import sys<br>
from osgeo import ogr,osr,gdal<br>
from osgeo.gdalconst import *<br>
<br>
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])<br>
pixelWidth = 0.000899322046076<br>
pixelHeight = -0.000899322046076<br>
<br>
#create polygon<br>
ring = ogr.Geometry(ogr.wkbLinearRing)<br>
for point in points:<br>
ring.AddPoint(point[0],point[1])<br>
poly = ogr.Geometry(ogr.wkbPolygon)<br>
poly.AddGeometry(ring)<br>
<br>
#create memory layer for polygon<br>
srWGS84 = osr.SpatialReference()<br>
srWGS84.ImportFromProj4('+proj=longlat +datum=WGS84')<br>
polyds = ogr.GetDriverByName('Memory').CreateDataSource('')<br>
polylyr = polyds.CreateLayer('poly',geom_type=ogr.wkbPolygon,
srs=srWGS84)<br>
feat = ogr.Feature(polylyr.GetLayerDefn())<br>
feat.SetGeometryDirectly(poly)<br>
polylyr.CreateFeature(feat)<br>
<br>
#check polygon layer extent<br>
extent = polylyr.GetExtent()<br>
print '========Polygon Layer========'<br>
print 'UL: ', extent[0], extent[3]<br>
print 'LR: ', extent[1], extent[2]<br>
print 'Spatial Reference:\n',polylyr.GetSpatialRef(),"\n"<br>
x = extent[0]<br>
y = extent[3]<br>
xLR = extent[1]<br>
yLR = extent[2]<br>
xSpan = int((xLR-x)/pixelWidth)<br>
ySpan = int((yLR-y)/pixelHeight)<br>
<br>
<br>
#create mask from polygon<br>
mask = gdal.GetDriverByName('MEM').Create('',xSpan,ySpan,1,GDT_Byte)<br>
transform = [x,pixelWidth,0.0,y,0.0,-pixelHeight]<br>
mask.SetGeoTransform(transform)<br>
mask.SetProjection(srWGS84.ExportToWkt())<br>
projstring = srWGS84.ExportToPrettyWkt()<br>
mask.GetRasterBand(1).Fill(255)<br>
<br>
#check mask extent<br>
transform = mask.GetGeoTransform()<br>
xMask = transform[0] #top left of coverage<br>
yMask = transform[3] #top left of coverage<br>
maskPixelWidth = transform[1]<br>
maskPixelHeight = transform[5]<br>
print "========Mask Raster========"<br>
print "UL: ", xMask, yMask<br>
print "LR: ", xMask+(xSpan*maskPixelWidth), yMask+(ySpan*maskPixelWidth)<br>
print "Spatial Reference:\n", projstring<br>
print "Stats before burn: ", mask.GetRasterBand(1).GetStatistics(0,1)<br>
<br>
#burn polygon onto mask<br>
err = gdal.RasterizeLayer(mask,[1],polylyr,burn_values=[0])<br>
if err != 0:<br>
print(err)<br>
gdaltest.post_reason( 'got non-zero result code from
RasterizeLayer' )<br>
sys.exit(1)<br>
print "Stats after burn: ", mask.GetRasterBand(1).GetStatistics(0,1)<br>
gdal.GetDriverByName('GTiff').CreateCopy('mask.tif',mask)<br>
</tt></big><br>
</body>
</html>