0

I am trying to rasterize a polygon layer currently as a feature class in a file geodatabase using OGR/GDAL (Python script). So far I have managed to do the rasterization with the output tif file having the right projection. However, I don't seem to be able to set the default value (any value outside polygons) to any value I want (it's always 0) and also to compress the output file (the output file produced has 1GB, whereas the rasterized file produced, for instance, by ArcGIS Pro has 16KB).

This is the code I am trying with:

import os
from osgeo import gdal
from osgeo import ogr
from osgeo import osr
from osgeo import gdalconst

fc_path = r'E:\...\Processing\INPUT\DefendedAreas.gdb'
tif_path = r'E:\...\Processing\INPUT\template.tif'

ras = gdal.Open(tif_path, gdalconst.GA_ReadOnly)
geo_transform = ras.GetGeoTransform()

x_min = geo_transform[0]
y_max = geo_transform[3]
x_max = x_min + geo_transform[1] * ras.RasterXSize
y_min = y_max + geo_transform[5] * ras.RasterYSize
x_res = ras.RasterXSize
y_res = ras.RasterYSize
pixel_width = geo_transform[1]

driver = ogr.GetDriverByName("OpenFileGDB")
gdb = driver.Open(fc_path, 0)
polygs = gdb.GetLayer("my_layer")
SpatialRef = polygs.GetSpatialRef()

out_tif = r'E:\...\Processing\INPUT\test.tif'
target = gdal.GetDriverByName('GTiff').Create(out_tif, x_res, y_res, 1, gdal.GDT_Byte)
target.SetGeoTransform((x_min, pixel_width, 0, y_max, 0, -pixel_width))
target.SetProjection(SpatialRef.ExportToWkt())
band = target.GetRasterBand(1)
NoData = -1
band.SetNoDataValue(NoData)
band.FlushCache()

gdal.RasterizeLayer(target, [1], polygs, options=["ATTRIBUTE=SoP", "burnValues=[-1]",  "creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256']"])
target = None

As you can see I have 3 arguments in the options list of gdal.RasterizeLayer:

"ATTRIBUTE=SoP"   (it works)
"burnValues=[-1]"   (it's being ignored)
"creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256']"   (it's being ignored)

NOTE: Does anyone know any good GDAL/OGR documentation? I find "https://gdal.org/python/index.html" and "https://pcjericks.github.io/py-gdalogr-cookbook/index.html" not very helpful.

3
  • 1
    There are thousands of working Python examples in the GDAL autotests in GitHub github.com/OSGeo/gdal/tree/master/autotest. They are not written to be used as tutorial but they are run after each commit to codebase so they must be in good shape.
    – user30184
    Commented Jan 27, 2022 at 10:54
  • Thanks for that link. However, I could not find any example doing what I'm trying to do. I would like to know why my code doesn't do what I would expect it to do. Commented Jan 27, 2022 at 11:40
  • 1
    There is also another route with github.com/OSGeo/gdal/blob/master/autotest/utilities/…. Otherwise I think that you may set TIFF options too late. Don't you have target TIFF created already when you get to gdal.RasterizeLayer(target,?
    – user30184
    Commented Jan 27, 2022 at 11:57

1 Answer 1

0

I think that the compression options need to be passed on to the output dataset creation bit (e.g.):

target = gdal.GetDriverByName('GTiff').Create(
         out_tif, x_res, y_res, 1, gdal.GDT_Byte,
         options=['COMPRESS=LZW'] 
         )

The two options you pass to the rasterise command are exclusive: you either burn the value of the field SoP, or you burn -1. If you want to select the polygons. Or do you want to burn -1 if the value of SoP is something in particular?

The documentation for gdal is... well, could be... better ;-) But this will give you an idea: https://gdal.org/programs/gdal_rasterize.html. the Python bindings are kind of wrappers around the command line utilities.

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.