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.
<VRTDataset rasterXSize="1620" rasterYSize="1386">
<SRS>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"]]</SRS>
<GeoTransform> -3.0531428271702840e+01, 3.7890083929483308e-02, 0.0000000000000000e+00, 6.7079735828607269e+01, 0.0000000000000000e+00, -3.7890083929483308e-02</GeoTransform>
<VRTRasterBand dataType="Byte" band="1">
<NoDataValue>0</NoDataValue>
<ColorInterp>Gray</ColorInterp>
<ComplexSource resampling="average">
<SourceFilename relativeToVRT="1">a.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="1272" RasterYSize="791" DataType="Byte" BlockXSize="1272" BlockYSize="6" />
<SrcRect xOff="0" yOff="0" xSize="1272" ySize="791" />
<DstRect xOff="183.541791108252" yOff="0" xSize="1436.01175091236" ySize="892.991584097231" />
<NODATA>0</NODATA>
</ComplexSource>
<ComplexSource resampling="average">
<SourceFilename relativeToVRT="1">b.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="1166" RasterYSize="1007" DataType="Byte" BlockXSize="1166" BlockYSize="7" />
<SrcRect xOff="0" yOff="0" xSize="1166" ySize="1007" />
<DstRect xOff="0" yOff="508.697635340442" xSize="1015.655894997" ySize="877.157363861048" />
<NODATA>0</NODATA>
</ComplexSource>
</VRTRasterBand>
</VRTDataset>
Raster after adding pixel function.
<VRTDataset rasterXSize="1620" rasterYSize="1386">
<SRS>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"]]</SRS>
<GeoTransform> -3.0531428271702840e+01, 3.7890083929483308e-02, 0.0000000000000000e+00, 6.7079735828607269e+01, 0.0000000000000000e+00, -3.7890083929483308e-02</GeoTransform>
<VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
<PixelFunctionType>average</PixelFunctionType>
<PixelFunctionLanguage>Python</PixelFunctionLanguage>
<PixelFunctionCode><![CDATA[
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)
]]>
</PixelFunctionCode>
<NoDataValue>0</NoDataValue>
<ColorInterp>Gray</ColorInterp>
<ComplexSource>
<SourceFilename relativeToVRT="1">a.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="1166" RasterYSize="1007" DataType="Byte" BlockXSize="1166" BlockYSize="7" />
<SrcRect xOff="0" yOff="0" xSize="1166" ySize="1007" />
<DstRect xOff="0" yOff="508.697635340442" xSize="1015.655894997" ySize="877.157363861048" />
<NODATA>0</NODATA>
</ComplexSource>
<ComplexSource>
<SourceFilename relativeToVRT="1">b.tif</SourceFilename>
<SourceBand>1</SourceBand>
<SourceProperties RasterXSize="1272" RasterYSize="791" DataType="Byte" BlockXSize="1272" BlockYSize="6" />
<SrcRect xOff="0" yOff="0" xSize="1272" ySize="791" />
<DstRect xOff="183.541791108252" yOff="0" xSize="1436.01175091236" ySize="892.991584097231" />
<NODATA>0</NODATA>
</ComplexSource>
</VRTRasterBand>
</VRTDataset>
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.
