Usando Rasterio o GDAL para apilar múltiples bandas sin usar comandos de subproceso

¿Alguien tiene una forma elocuente de apilar múltiples archivos .tif en una pila de múltiples bandas usando Rasterio y / o GDAL?

Estoy buscando una forma de evitar el uso de un comando de subproceso como gdal_merge.py y más bien tenerlo como parte de mi script de Python.

Sé que tanto Rasterio como GDAL leen los archivos .tif como matrices, pero ¿cómo apilo esas matrices y escribo el resultado como bandas separadas?

Usando rasteriousted podría hacer

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

Por supuesto, supone que sus capas de entrada ya comparten la misma extensión, resolución y tipo de datos

Si usa GDAL 2.1+ es tan simple como gdal.BuildVRTentonces gdal.Translate:

from osgeo import gdal
outvrt = '/vsimem/stacked.vrt' #/vsimem is special in-memory virtual "directory"
outtif = '/tmp/stacked.tif'
tifs = ['a.tif', 'b.tif', 'c.tif', 'd.tif'] 
#or for all tifs in a dir
#import glob
#tifs = glob.glob('dir/*.tif')

outds = gdal.BuildVRT(outvrt, tifs, separate=True)
outds = gdal.Translate(outtif, outds)