[GIS] Writing array to vrt

gdalosgeopythonvrt

I would like to write an array from a single band image to a 3 band vrt where each band has the ColorInterpretation() tag of red, green, blue respectively.

My code below;

from osgeo import gdal
import os

#IM = "/path/to/image.tif"
IM = im

### read image ###
ds = gdal.Open(IM)
X = ds.RasterXSize
Y = ds.RasterYSize
band = ds.GetRasterBand(1).ReadAsArray()

### write to 3 bands ###
driver = gdal.GetDriverByName("vrt")
outPath = os.path.join(os.path.split(IM)[0], "test_image.vrt")
outIM = driver.Create(outPath, X, Y, 0, gdal.GDT_Int16)
for i in range(1, 4):
    outIM.AddBand()
    outIM.GetRasterBand(i).SetRasterColorInterpretation(2 + i)
    outIM.GetRasterBand(i).WriteArray(band)
    print outIM.GetRasterBand(i).ReadAsArray()
    print outIM.GetRasterBand(i).GetRasterColorInterpretation()
outIM = None

The output is:

[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
3
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
4
[[0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 ..., 
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]
 [0 0 0 ..., 0 0 0]]
5

So setting the ColorInterpretation() (unlike when using GTiff, see here), but all values are zero. Why?!

#### UPDATE ####

I have been running this through ipython and just noticed this error on the command line:

Error 1: Writing through VRTSourcedRasterBand is not supported

Best Answer

You can't write data to a VRT. It's a "Virtual Raster", and is actually just an XML file. Try something like:

gdal_translate -of VRT -b1 -b1 -b1 in.tif out.vrt

Then open the out.vrt in python and set the colour interpretation.

You can write the vrt xml in python instead of using the commandline, but it's easier to use the subprocess module to run the gdal_translate command.

Related Question