python - Using Rasterio or GDAL to stack multiple bands without using subprocess commands

Does anyone have a eloquent way of stacking multiple .tif files into a multiple band stack using Rasterio and/or GDAL?

I am looking for a way to avoid using a subprocess command like gdal_merge.py and rather have it as part of my python script.

I know that both Rasterio and GDAL read the .tif files as arrays, but how do I stack those arrays and write out the result as separate bands?

Answer

Using rasterio you could do

import rasterio

file_list = ['file1.tif', 'file2.tif', 'file3.tif']

# Read metadata of first file
with rasterio.open(file_list[0]) as src0:
    meta = src0.meta


# Update meta to reflect the number of layers
meta.update(count = len(file_list))

# Read each layer and write it to stack
with rasterio.open('stack.tif', 'w', **meta) as dst:
    for id, layer in enumerate(file_list, start=1):
        with rasterio.open(layer) as src1:
            dst.write_band(id, src1.read(1))

It of course assumes that your input layers already share the same extent, resolution and data type