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?

4 Answers 4


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

  • 1
    Yes, this is essentially what Rasterio's rio-stack program does: github.com/mapbox/rasterio/blob/master/rasterio/rio/….
    – sgillies
    Commented Jan 11, 2017 at 7:25
  • 1
    Is it efficient to keep the stack in memory (to perform several functions on the various bands) instead of writing the stacked file out? Or should it be written to a file and then manipulated?
    – Shawn
    Commented Jan 25, 2019 at 10:21
  • alas I get this error "RasterioIOError: '/' not recognized as a supported file format."
    – ilFonta
    Commented Jul 31, 2019 at 19:37
  • @ilFonta, make a new question with a reproducible, minimal code example then. Commented Jan 8, 2020 at 15:41
  • @sgillies I tried a few passes at using rio stack but couldn't figure out how to get it to work in this example. Assuming we have 3 tifs R,G, and B, how would it work? I'm getting squat
    – mmann1123
    Commented Jan 25, 2022 at 1:45

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)
  • if the every file of tifs list is mutli-bands file, and diferent band in every files need to be stacked to outfile, meanwhile sparate is True. This code is not work Commented Jun 14, 2020 at 8:29
  • 1
    That scenario is not part of the question. If you have a new question, then ask it separately.
    – user2856
    Commented Jun 14, 2020 at 10:48
  • i am sorry to rashly comment. my purpose is to extend this question for solving these kind of layer stacking by gdal. I am sorry again Commented Jun 14, 2020 at 11:17

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

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

(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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.