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.9371001, -33.8374002], [25.2897,
-33.6713998], [26.1754, -33.4153002],
[27.8372999, -33.0525002], [27.6933001,
-33.1518998], [25.475, -33.9151999],
[25.2382999, -33.8808997], [24.9371001,
-33.8374002])
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)
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
