gdal RasterizeLayer gives different extent to arcpy PolygonToRaster

Question

I have an existing ArcPy script which uses PolygonToRaster_conversion to convert a polygon shapefile to raster (tiff). I am now exploring GDAL/OGR and have used the cookbook approach to carry out the same task but find two key differences in the output rasters:

1. The extent of the GDAL output raster is different from the arcpy output. Xmin is the same in each. However, the number of rows and cols are one less each in the gdal rasters, and the Y values differ by an odd amount (which is not a multiple of grid pixel size). I have tried the rasterize function with and without the ALL_TOUCHED option to no avail. I suspectthe difference is in calculating the extent, and applying the GeoTransform. Any clues?

2. The NoData value set using gdal is neither correctly set or recognised in ArcGIS or QGIS. Instead of a NoData value, the pixels have a value of zero. Could this be to do with the output raster datatype?

My code using gdal:

    NoDataVal = -9999

    # Open the data source and read in the extent
    shpDS = ogr.Open(inPolygonShp)
    shpLayer = shpDS.GetLayer()

    # Create the destination data source
    xMin, xMax, yMin, yMax = shpLayer.GetExtent()
    xRes = int((xMax - xMin) / inGridSize)
    yRes = int((yMax - yMin) / inGridSize)
    rasterDS = gdal.GetDriverByName('GTiff').Create(outputRaster, xRes, yRes, 1, gdal.GDT_Byte)

    # Define spatial reference
    rasterDS.SetProjection(shpLayer.GetSpatialRef().ExportToWkt())
    rasterDS.SetGeoTransform((xMin, inGridSize, 0, yMax, 0, -inGridSize))

    rBand = rasterDS.GetRasterBand(1)
    rBand.SetNoDataValue(NoDataVal)
    rBand.Fill(NoDataVal)

    # Rasterize
    err = gdal.RasterizeLayer(rasterDS, [1], shpLayer, burn_values=[1], options = ["ALL_TOUCHED=TRUE"])     # http://www.gdal.org/gdal_rasterize.html 

and using arcpy:

    idFld = arcpy.Describe(inPolygonShp).OIDFieldName
    addMsgAndPrint("OID Field used: "+ idFld)
    arcpy.PolygonToRaster_conversion(inPolygonShp, idFld, outputRaster, "CELL_CENTER", "NONE", inGridSize)

Answer 1

I have resolved the extent and y dimension overlap with two modifications to the above. To deal with size difference, I realised the 'int()' call will always truncate the division result, so simply use:

    xRes = int((xMax - xMin) / inGridSize) + 1
    yRes = int((yMax - yMin) / inGridSize) + 1

and to deal with the Y dimension offset, calculate the yMax from the yMin and grid size using:

    rasYMax = yMin + (yRes * inGridSize)
    rasterDS.SetGeoTransform((xMin, inGridSize, 0, rasYMax, 0, -inGridSize))

Now the gdal raster is a perfect match for the arcpy one!

I am still curious about the NoData issue though.

---

Don't use the options = ["ALL_TOUCHED=TRUE"] in the gdal.RasterizeLayer for arcpy equivalence.

Answer 2

You defined your NoData value as -9999 and your raster image is defined as gdal.GDT_Byte. The GDT_Byte specified is an eight bit unsigned integer, which can only contain values between 0 to 255, and therefore the NoData value automatically reverts to 0.