How to merge overlapping raster data sets in GDAL for Python to select maximum value at each pixel

Question

I am trying to merge a variable number of geotiffs into a single, clipped output geotiff, using GDAL 's Python bindings (or via command line).

Requirements:

• there will be 2 to N input files each time
• the inputs are all in a common projection, with different (but overlapping) extents
• the input files contain no_data pixels
• I want to keep the maximum input value at each each output pixel
• I want to clip the final image to a known extent
• the final output will be in the same projection as the input files

Ideally, I'd like to use something like gdal.Warp on a VRT that is built using the input files, to avoid having to create intermediate files/arrays. However, my experience and research indicate that a VRT will use pixels based on the order in which the input files are specified, which is also the case for gdal_merge and gdal_translate.

This seems like it should be a common task, but I've been unable to find a solution that works. The most likely looking approach I've found is to use a PixelFunction, based on this example: Combining overlapping rasters by minimum value in python. Failing so, far though. I have implemented this:

    ########### Create the VRT ###########
    vrt_file = 'mosaic.vrt'
    cmd = ['gdalbuildvrt', vrt_file]
    for file in file_list:
        cmd.append(file) # add all the files to be mosaicked

    proc = subprocess.Popen(cmd)
    stdout,stderr=proc.communicate()
    exit_code=proc.wait()

    if exit_code:
        raise RuntimeError(stderr)
    else:print(stdout)

    ########### Add the Pixel Function ###########
    try_pixel_function = True
    gdal.SetConfigOption('GDAL_VRT_ENABLE_PYTHON', 'YES')

    if try_pixel_function: # FAILS WHEN THIS BLOCK IS RUN
        # Adapted from https://gis.stackexchange.com/questions/350233/how-can-i-modify-a-vrtrasterband-sub-class-etc-from-python
        from lxml import etree
        vrt_tree = etree.parse(vrt_file)
        vrt_root = vrt_tree.getroot()
        vrt_band = vrt_root.findall(".//VRTRasterBand[@band='1']")[0]
        vrt_band.set('subClass', 'VRTDerivedRasterBand')
        pixelFunctionType = etree.SubElement(vrt_band, 'PixelFunctionType')
        pixelFunctionType.text = 'find_max'
        pixelFunctionLanguage = etree.SubElement(vrt_band, 'PixelFunctionLanguage')
        pixelFunctionLanguage.text = 'Python'
        pixelFunctionCode = etree.SubElement(vrt_band, 'PixelFunctionCode')
        pixelFunctionCode.text = etree.CDATA("""
        import numpy as np

        def find_max(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
            np.amax(in_ar, axis=0, initial=255, out=out_ar)
        """)
        vrt_tree.write(vrt_file)

    ########### get vrt's projection info ###########
    ds=gdal.Open(vrt_file)
    prj=ds.GetProjection()
    print(prj)

    ########### create mosaic ###########
    outfile_name = 'scratch.tif'
    cmd = ['gdalwarp', '-s_srs', prj, '-t_srs', 'EPSG:3857', '-r', 'max', '-of', 'VRT', vrt_file, outfile_name, '-te_srs', 'EPSG:4326', '-te', str(min_lon), str(min_lat), str(max_lon), str(max_lat), '-ts', str(resolution), str(resolution), '-et', '0', '-ot', 'Int16', '-ovr', 'NONE', '-srcnodata', '0', '-overwrite']

    proc = subprocess.Popen(cmd)
    stdout,stderr=proc.communicate()
    exit_code=proc.wait()

    if exit_code:
        raise RuntimeError(stderr)
    else:print(stdout)

    ########### retrieve the mosaicked data ###########
    ds = gdal.Open(outfile_name)
    data = ds.GetRasterBand(1).ReadAsArray()
    gdal.SetConfigOption('GDAL_VRT_ENABLE_PYTHON', 'NO')

This works fine when try_pixel_function is false (except that the last file overwrites everything else, as mentioned above), but fails on the last line, data = ds.GetRasterBand(1).ReadAsArray(), when try_pixel_function is true, with the error:

ERROR 1: Couldn't compile code: <class 'IndentationError'>, ('unexpected indent', ('gdal_vrt_module_0x287682af0', 6, 8, '
import numpy as np\n', 6, -1)) ERROR 1: scratch.tif, band 1: IReadBlock failed at X offset 0, Y offset 0: Couldn't compile code: <class 'IndentationError'>, ('unexpected indent', ('gdal_vrt_module_0x287682af0', 6, 8, ' import numpy as np\n', 6, -1))

I'm guessing the indent error refers to the code I'm trying to add using pixelFunctionCode.text = etree.CDATA(...

Any suggestions on what I'm doing wrong?

Other research: I also looked at using gdal_calc, as in this example from https://gdal.org/programs/gdal_calc.html

gdal_calc.py -A input1.tif -B input2.tif -C input3.tif --outfile=result.tif --calc="numpy.max((A,B,C),axis=0)

This seems doable but not super performant. I may have up to about 20 input files. I would have to either run on the inputs (which could involve wasting a lot of computation on the overlap areas that will be outside the desired final extent), OR clip each input individually to the final output extent before running this max operation.

(Most relevant results I've found:)

• Merge overlapping rasters using GDAL
• Combining overlapping rasters by minimum value in python

Answer 1

I think the error comes in this part:

        pixelFunctionCode.text = etree.CDATA("""
        import numpy as np

        def find_max(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
            np.amax(in_ar, axis=0, initial=255, out=out_ar)
        """)

whereas it should be:

        pixelFunctionCode.text = etree.CDATA("""import numpy as np
def find_max(in_ar, out_ar, xoff, yoff, xsize, ysize, raster_xsize, raster_ysize, buf_radius, gt, **kwargs):
    np.amax(in_ar, axis=0, initial=255, out=out_ar)""")

Note: Because of browser I indented with 4 * " " (whitespace). This might lead to problems, so check this if it doesnt work.