I have multiple rasters which overlap eachother geospatially but do not completely occupy the same space. I can combine them with gdalwarp like:
gdalwarp -s_srs EPSG:4326 *.tif larger_output.tif
In this method gdalwarp uses the order of the input files to determine which pixel ends up in the final image in the case that two pixels occupy the same location in the x,y plane.
Is it possible to tell gdalwarp to use the max(intensity(x,y)) of a given pixel that overlaps another pixel in the input as a means of determining which one ends up in the final image?
I presume the code behind gdalwarp works on the images sequentially or at least in separate threads and then combines them one at a time into a final raster? This addition might be as simple as an if statement that before replacing a given pixel in the final raster would compare intensities (if there were already a value there).
Does something like this already exist or should I get into the code to try something out?
Best Answer
I started digging into the gdal source found where I could make the changes but figured that I might not want to do it this way after all. What I decided to do instead (but have yet to do) was:
Iterate through the images, grabbing the bbox of each, then based on each bbox make groups of non-overlapping images. Then use gdalwarp to combine each non-overlapping set (since they dont overlap it doesnt matter what order you send them to gdalwarp) into a single image as above. Whilst looking at each bbox above I calculate a bbox that is the extrema of all bboxes, so that when I do the gdalwarp of each non-overlapping set, I end up with images for each set that are the same dimensions and have the same origin. Finally to combine them taking the maximum intensity at each pixel I use imagemagick's composite with the -compose lighten option. I think I might even be able to forgo making all the images the same size with the same origin if I use imagemagick's -geometry option to specify an offset.