Getting max feature value with GDAL Python RasterizeLayer

Question

I'd like to know how to use the Python/GDAL function RasterizeLayer to create a raster layer, where each pixel give the maximum value of a integer attribute. According to docs and forums, if I use :

gdal.RasterizeLayer(out_raster_ds,  # output to our new dataset
                             [1],  # output to our new dataset's first band
                             layer,  # rasterize this layer
                             None, None,  # don't worry about transformations since we're in same projection
                             [0],  # burn value 0
                             ['ALL_TOUCHED=TRUE',  # rasterize all pixels touched by polygons
                              'ATTRIBUTE=id']  # put raster values according to the 'id' field values
                             )

This will use the value of "id". But in my case, each raster cell intersects many spatial data (from PostGIS), and I want to keep the highest value.

Any idea ?

Answer 1

I would also be very interested in a nice an easy solution for this problem.

We also had to do something similar and ended using a tricky solution because we could not find any other solution using GDAL. Hope this might help...

In our case we also had to output the individual layers in individual GeoTIFF before merging them together using the maximum value where there is overlap.

1. Create all the individual GeoTIFF, using gdal.RasterizeLayer (No overlap here because you rasterize the individual layers).
2. Use gdalwarp to reproject the individual GeoTIFF with the same spatial resolution and extent.
3. Use gdal_calc.py to create a mosaic using the maximum value (gdal_calc.py -A input1.tif -B input2.tif --output=outfile.tif --calc="maximum(A,B)") where there is overlap, in your case where the are many layer values for a pixel.
4. Use gdal_translate to set 0 value pixels to no data (in my case, 0 is no data).

Note that :

gdal_calc.py takes a maximum of 26 inputs. It is also important to note that gdal_calc.py normally only outputs the overlapping zone of the rasters, this is why we must reproject (with gdalwarp) the GeoTIFF to have them in the same extent.

When reprojecting the rasters with gdalwarp, I use -srcnodata None because gdal_calc.py does not output the raster where there are no data. You can then set your no data value in a final step using gdal_translate if needed.

This solution works, but it's far from been ideal.

If you have gdal 2.2 installed

You could try using vrt build in pixel function. this solution was provided to us by Rutger. See: https://lists.osgeo.org/pipermail/gdal-dev/2017-November/047758.html

Answer 2

I have rasterized the maximum or minimum attribute value of overlapping polygons by first sorting the attribute value of interest and then resetting the index using the python package geopandas. GDAL will rasterize the features in ascending index order using whatever attribute value you specify for the cell value. So if you sort ascending the maximum value of any overlapping polygons will be the output raster cell value.

import geopandas as gp

gdf = gp.read_file(shp_path)
gdf.sort_values(by=['your_attribute'], ascending=True, inplace=True) #ascending=False for min
gdf.reset_index(inplace=True)
gdf.to_file(sorted_shp_path) 

After that, you can read the sorted shapefile back in and use the code in your question to perform the rasterization.