Using GDAL in python to stack Landsat bands

Question

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.

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',
             'LC80430342013291LGN00_B2.TIF',
             'LC80430342013291LGN00_B3.TIF',
             'LC80430342013291LGN00_B4.TIF',
             'LC80430342013291LGN00_B5.TIF',
             'LC80430342013291LGN00_B6.TIF',
             'LC80430342013291LGN00_B7.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/