Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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

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.

share|improve this question
    
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
add comment

1 Answer 1

up vote 4 down vote accepted

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'])
share|improve this answer
    
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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