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.