2

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?

4
  • 1
    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
    – mikewatt
    Nov 19, 2021 at 19:35
  • 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.
    – user2856
    Nov 26, 2021 at 5:35
  • @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 Nov 26, 2021 at 18:48
  • @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. Nov 26, 2021 at 18:50

1 Answer 1

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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