Converting GDAL command line statements to Python

Question

I am trying to carry out the task of sending each row/polygon from a single shapefile to its own unique raster. I so far use GDAL commands through the OSGeo4W Shell like so after navigating to my directory and entering py3_env:

FOR /L %i IN (1,1,250) DO gdal_rasterize -a Value -l Polygons -where "Polygon_ID=%i"  -tr 100 100 Polygons.shp Polygon_raster%i.tif

What this does is loop through all the polygons (rows) in my Polygons.shp shapefile and send them to their own unique raster, burning in the “Value” values. While this method is nice and straightforward, it does not let me do more complex tasks like multiplying each output raster by a specific value. What I specifically want to do is this:

Each polygon (row) in the shapefile contains both a “Value” column and a “Value2” column. I want to create all of my individual rasters like above, but I want to multiply each raster by the associated “Value2” value for that specific polygon. I do not know how to do this through OSGeo4W command line statements, and would much prefer doing this in python with Jupyter Notebook. However, I am having a lot of trouble directly translating my above gdal command statement into code that would work in Jupyter Notebook with python, since my command statement above uses a Windows for loop.

I would try this:

from osgeo import gdal
shapefile = polygons.shp
Value2 = shapefile[‘Value2’]
for polygon in shapefile:
     raster = gdal_rasterize(‘Value’)
     raster * Value2
     write_shapefile(raster, 'Polygon_raster%i.tif')

I want to do this in Jupyter Notebook using Python, but I am unsure how to directly translate my above command line into the necessary Python code that would work here, and particularly so that each raster would receive its own unique ID number (%i), which is straightforward with the command statement. How could I translate my command line statement into Python code that would let me loop through the rows in my shapefile, rasterize each, and multiply by its corresponding “Value2” value?

Answer 1

Here's a way that loops through polygons, creates a temp layer and rasterizes that to a new tif:

from osgeo import gdal

out_raster_template = "/path/to/test{}.tif"
shape_file = "/path/to/input.shp"

pixel_size = 0.00025  3 must use same units as shape file CRS
nodata = -9999

id_field = 'id'        # field to use to name each tif
value_field = 'VALUE'  # Field in shapefile to burn into each tif

drv = gdal.GetDriverByName("ESRI Shapefile")

shp_ds = gdal.OpenEx(shape_file, gdal.OF_VECTOR)
lyr = shp_ds.GetLayer()

xmin, xmax, ymin, ymax = lyr.GetExtent()
srs = lyr.GetSpatialRef()

feat_def = lyr.GetLayerDefn()

lyr.ResetReading()
for feat in lyr:
    id = feat.GetField(id_field)
    val = feat.GetField(value_field)

    tmp_feat = feat.Clone()  # Copy feature

    out_raster = out_raster_template.format(id)
    tmp_fn = '/vsimem/tmp.shp'    # temp in-memory dataset
    tmp_ds = drv.Create(tmp_fn, 0, 0, 0, gdal.GDT_Unknown )
    tmp_lyr = tmp_ds.CreateLayer(tmp_fn, None, feat_def.GetGeomType())
    tmp_lyr.CreateFeature(tmp_feat)
    tmp_feat, tmp_lyr, tmp_ds = None, None, None

    # tmp_ds = gdal.OpenEx(tmp_fn, gdal.OF_VECTOR)
    out_ds = gdal.Rasterize(out_raster, tmp_fn,
                   outputType=gdal.GDT_Float32, format='GTIFF', creationOptions=["COMPRESS=DEFLATE"],
                   noData=nodata, initValues=nodata,
                   xRes=pixel_size, yRes=-pixel_size, outputBounds=(xmin, ymin, xmax, ymax), outputSRS=srs,
                   allTouched=True, burnValues=val)

    out_ds = None
    drv.Delete(tmp_fn)