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?

• mikewatt: Have you tried to implement this with the actual Python APIs yet? It'd be easier to answer if we had that as a starting point instead of just pseudocode
• user2856: Do the polygons in shapefile overlap? If not, there's no need to create individual output rasters, just create 2 with each having the relevant value from the shapefile, you can then multiply the two output rasters together.
• LostinSpatialAnalysis: @mikewatt By Python API does that mean the python versions of the command line commands? The term API always confused me a bit, but I would assume that would mean using gdal.Rasterize() as opposed to running gdal_rasterize in a command line. Yes I have tried looking through this "GDAL/OGR Python API" resource here: gdal.org/python
• LostinSpatialAnalysis: @user2856 No, the polygons in this shapefile do not overlap. I see what you are suggesting there, and that definitely makes sense. I am just trying to create the individual rasters because the next step, which I did not mention in this post, is to run Euclidean distance on each of the individual rasters.

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)

enter image description here