Combining overlapping rasters by minimum value in python

Question

I would like to combine some 1-band raster images. I would like the overlap to be handled such that the smallest value in a pixel is chosen. All images have the same projection.

I tried looking at gdalwarp and gdal_merge (in command-line), but in the overlap they just use the values from the last image.

I have also seen suggestions to use PIL blend or paste, but these require an alpha layer, which indicates how the overlaps should be blended, which of course I don't have.

Is there anyone who knows how to make the overlap depend on the actual values in the images??

Answer 1

For as far as i know there is no way to do this with the standard command-line utilities. You try 'gdal_calc.py' which usually ships with GDAL, this requires the raster to have the same dimensions, so perhaps some preprocessing is required (with gdalwarp for example). The minimum of two raster can be calculated with something like:

python gdal_calc.py -A file1.tif -B file2.tif --calc="minimum(A,B)" --outfile=res.tif

I'm not quite sure what happens in regions with no overlap, perhaps you need to add the nodata keyword, more info at: http://www.gdal.org/gdal_calc.html

You have to specify the full path to gdal_calc.py.

For future releases:

Starting from GDAL 2.2, there is support for VRT Pixel functions written in Python:

http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_derived_python

You would then be able to make a VRT looking like:

<VRTDataset rasterXSize="20" rasterYSize="20">
  <SRS>EPSG:26711</SRS>
  <GeoTransform>440720,60,0,3751320,0,-60</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTDerivedRasterBand">
    <PixelFunctionType>add</PixelFunctionType>
    <PixelFunctionLanguage>Python</PixelFunctionLanguage>
    <PixelFunctionCode><![CDATA[
import numpy as np
def add(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize,
                   raster_ysize, buf_radius, gt, **kwargs):
    np.round_(np.clip(np.minimum(in_ar, axis = 0, dtype = 'uint16'),0,255),
              out = out_ar)
]]>
    </PixelFunctionCode>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte.tif</SourceFilename>
    </SimpleSource>
    <SimpleSource>
      <SourceFilename relativeToVRT="1">byte2.tif</SourceFilename>
    </SimpleSource>
  </VRTRasterBand>
</VRTDataset>

Also starting from GDAL 2.2, there will be built in pixel functions, which is a great idea. I dont see aggregations like min, max, mean yet, but that should probably be easy to add once the infrastructure is there. http://www.gdal.org/gdal_vrttut.html#gdal_vrttut_derived_c