[gdal-dev] VRT creation from Python, adding existing geotiff as bands

[Even Rouault]

Rémi,

> Hey dear list,
> this is my first post here,
> so sorry if I'm asking something obvious.
> 
> (I noticed the 2 previous mail thread regarding python and VRT).
> 
> I wrote python functions to convert a numpy structured array to several
> gdal rasters.
> 
> A numpy structured array is a 2D pixel array, but each pixel has in fact
> several values of possibly different type (like band).
> For instance a pixel dtype can be [float32,uint16,Float64].
> 
> I want to convert this array to raster(s).
> First I tried to create one raster with several bands,
> but GeoTiff cannot handle bands with different pixel type (or I didn't find
> how).

This is a limitation of the TIFF format, and most formats actually.

> 
> I concluded that I must create several rasters, and put them all as bands
> in a VRT.
> 
> So far I have been able to create several geoTiff raster of various type
> corresponding to the numpy structured array.
> 
> Now, I don't know how to create from python the VRT that unites all this
> rasters.
> 
> When I try to write in a VRT raster band I get errors.

"write" what exactly : the VRT XML/structure or the pixel values ?

You can partially create VRT files with a mix of API and XML. See for example 
the following extract from autotest/gdrivers/vrtderived.py

    filename = 'tmp/derived.vrt'
    vrt_ds = gdal.GetDriverByName('VRT').Create(filename, 50, 50, 0)

    options = [
        'subClass=VRTDerivedRasterBand',  # ==> You would likely need 
VRTSourcedRasterBand heere
    ]
    vrt_ds.AddBand(gdal.GDT_Byte, options)

    simpleSourceXML = '''    <SimpleSource>
      <SourceFilename>data/byte.tif</SourceFilename>
      <SourceBand>1</SourceBand>
    </SimpleSource>'''

    md = {}
    md['source_0'] = simpleSourceXML

    vrt_ds.GetRasterBand(1).SetMetadata(md, 'vrt_sources')
    md_read = vrt_ds.GetRasterBand(1).GetMetadata('vrt_sources')
    vrt_ds = None


This is a bit documented in the VRT tutorial. See 
http://gdal.org/gdal_vrttut.html (search "In this example a virtual dataset is 
created with the Create() method, and adding bands and sources 
programmatically" in the page)

Note that datasets with bands of different datatypes are not always well 
supported in the various algorithms, utilities or API (mostly due to the lack 
of formats that support it and the ensuing lack of testing of such situations)

Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com