Uso de GDAL en python para apilar bandas de Landsat

Question

Estoy procesando una gran cantidad de datos del Landsat. Las escenas de Landsat se envían como una colección de tiffs, un archivo para cada banda espectral, además de un montón de máscaras. Necesito apilar algunas de las bandas en un solo tiff para su posterior procesamiento. Puedo hacer esto fácilmente en varios programas, pero me gustaría automatizar el proceso utilizando GDAL en python. Eventualmente necesitaré automatizar el apilamiento de cientos o miles de imágenes, por lo que la eficiencia será importante.

He escrito un script que es simple y perfectamente funcional pero se siente extremadamente hacky- abro los tiffs deseados, los convierto en arrays de numpy, apilo los arrays, y guardo el array apilado en un nuevo geotiff con metadatos tomados de la entrada. Aquí está una versión simplificada del código:

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"))

Si es posible, me gustaría evitar la conversión de los datos a arrays numpy, y simplemente apilar los objetos gdal.Dataset (b1 y b2). Esto haría el código más limpio y presumiblemente más eficiente. Parece que debería ser sencillo de hacer, pero no he podido encontrar ninguna información sobre cómo hacerlo.

Los enfoques alternativos que he encontrado en la investigación sería utilizar gdal_merge.py o vrts (como se sugiere aquí ). Ninguno de los dos me parece ideal; gdal_merge está diseñado para tareas mucho más complicadas que las mías y presumiblemente implica mucha sobrecarga, mientras que los vrts sólo se pueden utilizar en la línea de comandos. Sin embargo, si el enfoque que he esbozado más arriba no es factible, ¿podría uno de estos aprestos ser más eficiente/pitónico que el hack que ya he escrito?

Edición: mi código debe ejecutarse en un sistema Windows.

Answer 1

Aunque requeriría otra biblioteca (que actualmente sólo está disponible en OS X y Linux) podría utilizar RSGISLib ( http://rsgislib.org ), que está construido sobre GDAL para hacer esto. Hay una función para apilar bandas, como ejemplo:

#!/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)

También permite pasar una lista de nombres de bandas.

Sin embargo, si se va a instalar RSGISLib para apilar las bandas del Landsat, se podría instalar ARCSI, que apilaría las bandas y las convertiría en radiancia o reflectancia (si se quiere hacer esto). Más detalles están aquí: http://spectraldifferences.wordpress.com/2014/05/27/arcsi/