I am processing a large amount of landsat data. Landsat scenes ship as a collection of tiffs- one file for each spectral band, plus a bunch of masks. I need to stack some of the bands into a single tiff for subsequent processing. I can easily do this in various programs, but would like to automate the process using GDAL in python. I will eventually need to automate the stacking of hundreds or thousands of images, so efficiency will be important.

I have written a script which is simple and perfectly functional but feels extremely hacky- I open the desired tiffs, convert them to numpy arrays, stack the arrays, and save the stacked array to a new geotiff with metadata taken from the input. Here is a simplified version of the code:

from osgeo import gdal, gdal_array
import numpy as np
b1 = gdal.Open("LT50250232011160PAC01_sr_band1.tif")
b2 = gdal.Open("LT50250232011160PAC01_sr_band2.tif")
array1 = b1.ReadAsArray()
array2 = b2.ReadAsArray()
stacked = np.array([array1,array2])
gdal_array.SaveArray(stacked.astype("int"), "b12_stacked.tif", "GTiff", gdal.Open("LT50250232011160PAC01_sr_band1.tif"))

If possible, I'd like to avoid converting the data to numpy arrays, and just stack the gdal.Dataset objects (b1 and b2). This would make the code cleaner and presumably more efficient. It seems like it should be straightforward to do, but I have not been able to find any information on how to do it.

Alternate approaches which I have come across in research would be to use gdal_merge.py or vrts (as suggested here). Neither of these strikes me as ideal; gdal_merge is designed for tasks much more complicated than mine and presumably involves a lot of overhead, while vrts are only usable in the command line. However, if the approach I outlined above is not feasible, would one of these appraoches be more efficient/pythonic than the hack I have already written?

Edit: my code needs to be run on a Windows system.

  • What do you mean that vrts are only usable in command line? Gdalbuildvrt gdal.org/gdalbuildvrt.html should be able to stack your bands to vrt within a second and after that you can use the vrt file just as if it was a multiband tiffs. Howerver, I have no experience on mask bands, only regular 7-band Landsat scenes.
    – user30184
    Commented Sep 29, 2014 at 17:40
  • I mean that there is no gdalbuildvrt python function- I would need to use a subprocess, os.system or something like this. Which is not to say it's not an option- stacking gdal.Dataset objects just seems simpler, if it's possible.
    – Joe
    Commented Sep 29, 2014 at 18:20
  • 1
    There used to be gdal_vrtmerge.py. Could it be of some use for you? Some version is here trac.osgeo.org/gdal/attachment/ticket/1985/gdal_vrtmerge.py
    – user30184
    Commented Sep 29, 2014 at 18:30
  • 2
    my suggestion is to use the code you have! maybe clean it up, add some options and have it accept commandline arguments and you have yourself a perfectly reasonable script to use over and over again. As you say, it does what you need. If your source files are produced by another process of your making, you can streamline the process by perhaps saving them as compressed numpy arrays (only if they're intermediate data, that is) or hdr/flt rasters (binary and fast to load). Commented Sep 29, 2014 at 20:55
  • @user30184, gdal_vrtmerge.py appears to be writing a vrt, and I do need a tiff at the end bot for subsequent processing and passing to collaborators. My python-fu is not great enough to pull it apart without considerable headache, but it could be useful to others with the same problem.
    – Joe
    Commented Sep 30, 2014 at 13:58

1 Answer 1


Although it would require another library (which currently is only available on OS X and Linux) you could use RSGISLib (http://rsgislib.org), which is built on top of GDAL to do this. There is a function to stack bands, as an example:

#!/usr/bin/env python
import rsgislib
from rsgislib import imageutils
# Create list of images
imageList = ['LC80430342013291LGN00_B1.TIF',
# Create list of band names           
bandNamesList = ['Coastal','Blue','Green','Red','NIR','SWIR1','SWIR2']
# Set output image
outputImage = 'LC80430342013291LGN00_stack.tif'
# Set format and type
gdalFormat = 'GTiff'
dataType = rsgislib.TYPE_16UINT
# Stack
imageutils.stackImageBands(imageList, bandNamesList, outputImage, None, 0, gdalFormat, dataType)

It also allows passing in a list of band names.

However, if you were going to install RSGISLib to stack the Landsat bands you could just install ARCSI, which would stack the bands and convert to radiance or reflectance (if you were looking to do this). More details are here: http://spectraldifferences.wordpress.com/2014/05/27/arcsi/

  • These both look like good resources, but unfortunately all of the computers in my lab run Windows. I could set up a virtual box, but I need to be able to pass this code around to my advisors and labmates, so I'm not sure that's wise- we run into enough compatibility/consistency issues as it is.
    – Joe
    Commented Sep 30, 2014 at 13:52

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.