[gdal-dev] gdalbuildvrt for multiband files?

[jramm]

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.