GDAL RasterizeLayer does not Burn All Polygons to Raster?

Question

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.

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.

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 '12 at 21:47

Mike T · Accepted Answer · 2023-07-24T03:32:33.100

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

edited Jul 24 '23 at 03:32

answered Aug 17 '12 at 04:06

Mike T

42,095
10
126
187

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 '12 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 '12 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 '12 at 20:11
It does not! It works in 1.9.0. Thanks so much! – Lark Aug 20 '12 at 17:39
1

There is also quite a good recipe here: http://gis.stackexchange.com/a/16916/9942 – j08lue Jul 25 '16 at 08:42
Dear @MikeT, I have one doubt regarding RasterizeLayer. What is the return type of this function and does it return the rasterized dataset similar to C++ implementation of Rasterize – Navneet Srivastava Jul 23 '23 at 16:29
@NavneetSrivastava the Python function returns a status int, where zero is normal exit. It is wrapper function to the C implementation GDALRasterizeLayers. – Mike T Jul 24 '23 at 03:35
@MikeT Thanks. Is there any similar wrapper for GDALRasterize in python ? – Navneet Srivastava Jul 25 '23 at 04:11

GDAL RasterizeLayer does not Burn All Polygons to Raster?

1 Answers1

Linked