Treat overlapping pixel when doing merge GeoTIFF images using Rasterio

Question

I have a list of GeoTIFF images. I would like to merge these images together.

However, as my understanding, it will take in turn the file and at overlapping pixels will receive the values from the later files.

What do I need? I need the overlapping pixels to receive average values from all inputs.

Any solution to solve it, even with other packages such as GDAL or rioxarray?

import rasterio
from rasterio.merge import merge

input_files = ['img1.tif', 'img2.tif', 'img3.tif']
output_file = "merged.tif"

src_files_to_mosaic = []
for file in input_files:
    src = rasterio.open(file)
    src_files_to_mosaic.append(src)

mosaic, out_trans = merge(src_files_to_mosaic)

out_meta = src.meta.copy()
out_meta.update({"driver": "GTiff",
                 "height": mosaic.shape[1],
                 "width": mosaic.shape[2],
                 "transform": out_trans})

with rasterio.open(output_file, "w", **out_meta) as dest:
    dest.write(mosaic)

Answer 1

The rasterio.merge function has an argument where you can specify the merging method. Availables are:

• first: reverse painting
• last: paint valid new on top of existing
• min: pixel-wise min of existing and new
• max: pixel-wise max of existing and new

You can also give a custom function with signature:

• merged_data (array_like): array to update with new_data
• new_data (array_like): data to merge same shape as merged_data
• merged_mask, new_mask (array_like): boolean masks where merged/new data pixels are invalid same shape as merged_data
• index (int): index of the current dataset within the merged dataset collection
• roff (int): row offset in base array
• coff (int): column offset in base array

For example, here is a debug function for maximum (from here):

def copyto_max(merged_data, new_data, merged_mask, new_mask, **kwargs):
    mask = np.empty_like(merged_mask, dtype="bool")
    np.logical_or(merged_mask, new_mask, out=mask)
    np.logical_not(mask, out=mask)
    merged_data[mask] = np.maximum(merged_data, new_data, where=mask)[mask]
    np.logical_not(new_mask, out=mask)
    np.logical_and(merged_mask, mask, out=mask)
    np.copyto(merged_data, new_data, where=mask, casting="unsafe")

I bet you can construct your averaging function on this model.