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

Question

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 1

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

Answer 2

If using GDAL 2.1+ it's as simple as gdal.BuildVRT then 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)

Answer 3

With gdal and numpy the following is possible assuming the images are the same size (rows and columns) and have the same CRS (e.g. 10m bands from a Sentinel 2 scene).

import gdal
import numpy as np

file_list = ['band1.tif', 'band2.tif', 'band3.tif']

array_list = []

# Read arrays
for file in file_list:
    src = gdal.Open(file)
    geotransform = src.GetGeoTransform() # Could be done more elegantly outside the for loop
    projection = src.GetProjectionRef()
    array_list.append(src.ReadAsArray())
    src = None

# Stack arrays
stacked_array = np.stack(array_list, axis=0)
array_list = None

# Write to disk
driver = gdal.GetDriverByName('GTiff')
n, rows, cols = stacked_array.shape
dataset = driver.Create('output_file_name.tif', cols, rows, n,
                                gdal.GDT_Uint16)
dataset.SetGeoTransform(geotransform)
dataset.SetProjection(projection)

for b in range(1,n+1):
    band = dataset.GetRasterBand(b) # GetRasterBand is not zero indexed
    band.WriteArray(stacked_array[b-1]) # Numpy is zero indexed

dataset = None
stacked_array = None

Answer 4

(thanks to https://gis.stackexchange.com/a/299218/155670, https://gis.stackexchange.com/a/223915/155670)

import osgeo.gdal
from osgeo.gdal import Dataset

vrt_internal_name: str = '' # if non-empty and not the same as target_tif, vrt will be created
vrt: Dataset = osgeo.gdal.BuildVRT(vrt_internal_name, ['r_0.tif', 'r_1.tif', 'r_2.tif'])

target_tif: str = 'r.tif'
osgeo.gdal.Translate(target_tif, vrt)

For me to make BuildVRT work all paths (inputs and outputs) were needed to be resolved