13

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.

Image for Problem- Finished Raster Burns

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])

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?

The files that I was using can be found here.

1
  • Can you provide a link to 'activity_3.shp' and 'AOI_Raster.tif'? I want to see if I can recreate on my end.
    – Rich
    Aug 16, 2012 at 21:47

1 Answer 1

13

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

gdal.RasterizeLayer(raster_dataset, [1], shape_layer, None, None, [1], ['ALL_TOUCHED=TRUE'])
8
  • YES! It's apparently ['ALL_TOUCHED=TRUE'], though unfortunately, that only fixed polygon layers. My point shapefile layers are still being super wonky, and show up one pixel over from where they're placed.
    – Lark
    Aug 17, 2012 at 17:17
  • It ends up looking like this. It's in the same projection as the others, and I was hoping this would somehow magically fix it too, but it just seems to stubbornly burn one pixel off of where it's actually located.
    – Lark
    Aug 17, 2012 at 17:31
  • That certainly looks bug-worthy, where the burn point is offset by dx/2 and dy/2. I wonder if that bug still persists with the latest trunk.
    – Mike T
    Aug 17, 2012 at 20:11
  • It does not! It works in 1.9.0. Thanks so much!
    – Lark
    Aug 20, 2012 at 17:39
  • 1
    There is also quite a good recipe here: gis.stackexchange.com/a/16916/9942
    – j08lue
    Jul 25, 2016 at 8:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.