Two options: 1. Use the xml module in python to build up the VRT file. You would create a VRTRasterBand node for each band of each file. Note you will need to figure out the maximum extent and raster size (in pixels) first, so this will be easiest if you enforce the same pixel size and projection on your inputs
2. Create using the GDAL library. Something along the lines of drv = gdal.GetDriverByName("VRT") dst = drv.Create(....) for src in src_list: for i in range(src.RasterCount): dst.AddBand(...) Again, you would need to work out the size of the final raster beforehand. The start of a very naive python implementation using xml is below. Look at http://www.gdal.org/gdal_vrttut.html to see what further nodes and attributes need creating from xml.etree.ElementTree import * def stack(src_list): # Get the maximum extent ulx = [] uly = [] lrx = [] lry = [] xres = [] yres = [] for src in src_list: gt = src.GetGeoTransform() ulx.append(gt[0]) uly.append(gt[3]) lrx.append(gt[0] + (gt[1] * src.RasterXSize)) lry.append(gt[3] + (gt[5] * src.RasterYSize)) xres.append(xres) yres.append(yres) # Dont handle skewed images assert gt[2] == gt[4] == 0 # Check resolutions are all the same assert xres.count(xres[0]) == len(xres) assert yres.count(yres[0]) == len(yres) dst_ulx = min(ulx) dst_lrx = max(lrx) dst_uly = max(uly) # for north up image dst_lry = min(lry) nxpixels = int((dst_lrx - dst_ulx) / xres[0]) nypixels = int((dst_lry - dst_uly) / yres[0]) root = Element("VRTDataset") root.set("RasterXSize", str(nxpixels)) root.set("RasterYSize", str(nypixels)) gt_node = SubElement(root, "GeoTransform") gt_node.text = ", ".join([dst_ulx, xres[0], 0, dst_uly, 0, yres[0]]) for src in src_list: for i in range(src.RasterCount): bnd_node = SubElement(root, "VRTRasterBand") src_node = SubElement(bnd_node, "SimpleSource") pth_node = SubElement(src_node, "SourceFilename") pth_node.text = src.GetDescription() srcbnd_node = SubElement(src_node, "SourceBand") srcbnd_node.text = str(i+1) A possibly simpler option would be to use gdalbuildvrt to create a VRT for the first band of each dataset, then read the VRT using python and duplicate the relevant nodes to add in new bands, then write it out again. -- View this message in context: http://osgeo-org.1560.x6.nabble.com/gdal-dev-gdalbuildvrt-for-multiband-files-tp5272094p5272472.html Sent from the GDAL - Dev mailing list archive at Nabble.com. _______________________________________________ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev