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.
- Create all the individual GeoTIFF, using gdal.RasterizeLayer (No overlap here because you rasterize the individual layers).
- Use gdalwarp to reproject the individual GeoTIFF with the same spatial resolution and extent.
- 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.
- 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