Home > python > GDAL RasterizeLayer Doesn't Burn All Polygons to Raster

GDAL RasterizeLayer Doesn't Burn All Polygons to Raster

August 15Hits:2
Advertisement

I am trying to burn a shapefile to a raster using GDAL's RasterizeLayer. I pre-create an area of interest raster from a different shapefile, given a specific pixel size. This AOI then serves as a base for all following rasterizations (same number of collumns and rows, same projection and geotransform).

The issue occurs, however, when I go to burn the shapes to their own raster, based on the same pixel size and projections. The link below (don't have enough rep to post the image), shows the original shapefile in tan, and the dark pink where RasterizeLayer has burned data. The light pink is nodata values for the dark pink raster data. The grey is the AOI based on which the shapefile burn was completed.

Given the extents of the shapefile polygons, I would expect to see burn values in the bottom two corners, as well as the two pixels underneath the data that does show. Obviously, however, this is not the case.

GDAL RasterizeLayer Doesn't Burn All Polygons to Raster

As follows is the code that I used to generate these. All of the shapes were created using QGIS, and were all created in the same projection. (It should be noted that the gridding in the picture shown was just to give an idea of the pixel size I was using.)

from osgeo import ogr from osgeo import gdal  aoi_uri = 'AOI_Raster.tif' aoi_raster = gdal.Open(aoi_uri)  def new_raster_from_base(base, outputURI, format, nodata, datatype):      cols = base.RasterXSize     rows = base.RasterYSize     projection = base.GetProjection()     geotransform = base.GetGeoTransform()     bands = base.RasterCount      driver = gdal.GetDriverByName(format)      new_raster = driver.Create(str(outputURI), cols, rows, bands, datatype)     new_raster.SetProjection(projection)     new_raster.SetGeoTransform(geotransform)      for i in range(bands):         new_raster.GetRasterBand(i + 1).SetNoDataValue(nodata)         new_raster.GetRasterBand(i + 1).Fill(nodata)      return new_raster  shape_uri = 'activity_3.shp' shape_datasource = ogr.Open(shape_uri) shape_layer = shape_datasource.GetLayer()  raster_out = 'new_raster.tif'  raster_dataset = new_raster_from_base(aoi_raster, raster_out, 'GTiff',                                 -1, gdal.GDT_Int32) band = raster_dataset.GetRasterBand(1) nodata = band.GetNoDataValue()  band.Fill(nodata)  gdal.RasterizeLayer(raster_dataset, [1], shape_layer, burn_values=[1]) 

My question therefore, is: Is this a bug in GDAL? Or is RasterizeLayer burning data based on something other than just the presence or lack of a polygon within a specified pixel area?

Edit: The files that I was using can be found here.

Answers

I've been playing with GDALRasterizeLayers this week and have a pretty good idea of what it is doing. By default, it will rasterise a pixel if the pixel centre is within the polygon. If there is nothing in the centre, it won't be rasterised, even if there are parts of a polygon within the pixel limits. To allow the rasterising to work the way you intend, try the "ALL_TOUCHED" option:

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])

Related Articles

Copyright (C) 2017 ceus-now.com, All Rights Reserved. webmaster#ceus-now.com 14 q. 0.389 s.