gdalwarp - Averaging overlapping rasters with gdal

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 take max or average 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

Answer

The following approach worked pretty well.

First I build virtual raster.

gdalbuildvrt raster.vrt -srcnodata 0 -input_file_list paths.txt

paths.txt is file with following content:

a.tif
b.tif

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.

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
   -3.0531428271702840e+01,  3.7890083929483308e-02,  0.0000000000000000e+00,  6.7079735828607269e+01,  0.0000000000000000e+00, -3.7890083929483308e-02
  
    0
    Gray
    

      a.tif
      1
      
      
      
      0
    
    
      b.tif
      1

      
      
      
      0
    
  

Raster after adding pixel function.

  GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0],UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]]
   -3.0531428271702840e+01,  3.7890083929483308e-02,  0.0000000000000000e+00,  6.7079735828607269e+01,  0.0000000000000000e+00, -3.7890083929483308e-02
  
    average
    Python
    import numpy as np

def average(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,raster_ysize, buf_radius, gt, **kwargs):
    div = np.zeros(in_ar[0].shape)

    for i in range(len(in_ar)):
        div += (in_ar[i] != 0)
    div[div == 0] = 1

    y = np.sum(in_ar, axis = 0, dtype = 'uint16')
    y = y / div

    np.clip(y,0,255, out = out_ar)
]]>
        

    0
    Gray
    
      a.tif
      1
      
      
      
      0
    

    
      b.tif
      1
      
      
      
      0
    
  

And finally, transform it to raster using gdal_translate and gdal python option set to 'yes':

gdal_translate --config GDAL_VRT_ENABLE_PYTHON YES raster.vrt raster.tif

A result image for this example.