I am working on creating a class which will merge several georeferenced rasters into one using different strategies, essentially taking average
, max
, min
where the images are overlapping.
So far I've tried using gdalwarp
with --resample
parameter set to average
.
gdalwarp -srcnodata 0 -r average a.tif b.tif output.tif
But gdalwarp
just overlaps the images. I've tried other approaches with gdal_merge.py
and gdalbuildvrt
but they also simply overlap images, without taking average.
Reading gdal
dev list I've seen people taking following approach:
- reproject images to same dimensions, filling the rest with no data value
- filling no-data values with zeroes
- using
gdal-calc
to takemax
oraverage
on images
I wanted to try this approach but stumbled on a problem of changing dimensions of image with adding no-data value, i. e. the following command changed the whole image, instead of just inserting extra no-data pixels.
gdalwarp -ts 1591 1859 a.tif r1.tif
So my question are:
- Is there any other approach on how this averaging could be done?
- Is there any utility, preferably with GDAL, that could change dimension of the image by adding no-data value pixels to it?
Note: you can find sample files here https://drive.google.com/drive/folders/1cm8Y4WX03wn4XrNKOifYBhd13GqVNGdb?usp=sharing
Best Answer
The following approach worked pretty well.
First I build virtual raster.
paths.txt
is file with following content:Then I add a pixel function to it, as showed here https://lists.osgeo.org/pipermail/gdal-dev/2016-September/045134.html. Pixel function is written using numpy, basically it sums all images and divides each pixel by the number of overlapping images for that particular pixel.
Raster before adding pixel function.
Raster after adding pixel function.
And finally, transform it to raster using
gdal_translate
and gdal python option set to 'yes':A result image for this example.