I have a vector layer that contains many overlapping polygons, each containing a "Value" in the attribute table. I am trying to rasterize each polygon as an individual raster layer, and then run the raster calculator to find the average of all of the new raster layers. Each new raster layer would have the "Value" burned into the shape of the polygon, where the rest of the pixels in the layer are Null/No data. Each new raster layer would have the same extent. The reason for this is so that I can find the average of all the overlapping regions. Since I am working from a single polygons layer that contains dozens of polygons, I am using python, given the ability loop through files.

Here is what I am trying to code out in Python:

rasters_list = []
For row in geodataframe
    append to rasters_list

raster_calculator(Average, data = rasters_list)
send output raster to new file

Does this look to be on the right track conceptually? What this is saying is to create a geodataframe from my polygons layer, where each row is a separate polygon. For each row, create a new raster layer. Append all those new raster layers to a list. Then run the raster calculator, taking the average of all the raster layers in that list, and send the output to a new raster file.

How might I code this out with the correct functions and syntax? I am unsure how to actually iterate through each raster layer and send them to the raster calculator as inputs for the average.

  • Are you trying to create a heatmap of where your polygons overlap? If so, you could use gdal.Rasterize(output_raster, vector_input, options=options) where options = gdal.RasterizeOptions(format="GTiff', bands=[1], attribute="attribute_field", add=True) and skip the for loop entirely. See gdal.org/python/index.html for the docs Sep 23 '21 at 16:58

Your Answer

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

Browse other questions tagged or ask your own question.