Using GDAL/Python to stack georeferenced images of different sizes

I am in the process of porting a code I wrote in IDL (interactive data language) to Python but am running into a bit of a problem that I am hoping someone can help me with.

The code goes like this:

1. 
   

   take individual classified Landsat GeoTIFFs (say there are N individual 1-band files per scene, each representing a different day) and further reduce these images to three binary-themed 1-band images (water and not water, land and not land, water/land and not water/land). This will be done by reading the rasters as matrices and replacing values. ** I don't actually need to have these images, so I can save them as memory or just keep them as numpy ndarrays to move to the next step

   
   


2. 
   

   stack these images/arrays to produce 3 different (1 for each 'element') N-band stacks (or a 3-dimensional array-- (samples, lines, N)) for each scene

   
   



3. total the stacks to get a total number of water/land/water&land observations per pixel (produces one 1-band total image for each scene)


4. other stuff

The problem I am running into is when I get to the stacking, as the individual images for each scene vary in size, although they mostly overlap with each other. I originally used an ENVI layer-stacking routine that takes the N different-sized 1-band images for each scene and stacks them into an N-band image with an extent that encompasses all of the images' extents, and then reading the resulting rasters in as 3-D arrays to do the totals. I would like to do something similar with GDAL/Python but am not sure how to go about doing so. I was thinking I would implement GDAL capabilities of GeoTIFFs by using the geotransform info of the images to somehow find the inclusive extent, possibly padding the edges of the images with 0's so they are all the same size, stacking these images/3-d arrays so that they are correctly aligned, then computing the totals. Hopefully there is an easier way, as I'm not sure how to pull that off.

Does anyone have any suggestions or ideas as to what would be the most efficient way (or, any way really), to do what I need to do? I'm open to anything.

Answer

I'm not totally certain on what you need, but I think:

gdalbuildvrt -separate stack.vrt lsat1.tif lsat2.tif ...

Should give you a gdal dataset that is 'stacked and is covers the extent of all images. If you need a tif after that, use

gdal_translate stack.vrt stack.tif