GDAL – Averaging Overlapping Rasters [Closed]

gdalgdal-translategdalwarpoverlayraster

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

Best 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.

<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.

averaged image

Related Question