Is there a way to perform the same task as the gdalbuildvrt utility using GDAL Python bindings? So far I have not found any way to do this other than creating a vrt of a single dataset and manually editing the xml. I would like to create a vrt from multiple rasters (essentially performing a mosaic). Is this possible using pure Python? My other option is to use subprocess to simply call gdalbuildvrt.
Python GDAL – Equivalent of gdalbuildvrt for Raster Data
gdalgdalbuildvrtpythonrastervrt
Related Solutions
Not sure if it helps you, but this was my workflow to tile 2GB of 300 Dutch Topo maps to OSM compatible zoomlevel 15 and 16:
- Create vrts for each tif and expand indexed colours to RGBA:
for %%N in (D:\Karten\gdal\gdal2tiles\NL25*.tif) DO gdal_translate -of vrt -expand rgba %%N D:\Karten\gdal\gdal2tiles\NL25\%%~nN.vrt
- Create an index vrt for all files:
gdalbuildvrt -allow_projection_difference index25.vrt NL25*.vrt
- Use gdal2tiles to define source CRS, create zoomlevel 16, and mosaic that to zoomlevel 15:
gdal2tiles --s_srs EPSG:28992 --zoom 15-16 index25.vrt tiles
I had to make some minor changes to gdal2tiles to enable correct tile numbering, as described here: GDAL2Tiles: MapTiles from BSB/KAP are Switched
It took several days, but python just is not that quick...
Update:
As of GDAL 2.1, you can use the gdal.Translate() function and pass the bands in the bandList argument (list of arguments).
from osgeo import gdal
infn='in.tif'
outfn='out.tif'
gdal.Translate(outfn, infn, bandList=[1,2,3], **other_kwargs)
Original answer:
gdal_translate does this by building a VRT from scratch, copying all the metadata from the source dataset (unless everything matches, then it just uses CreateCopy).
You can do this in python but it's a bit fiddly. An easier way is use the VRT driver to create a VRT copy, then modify the VRT XML to remove the bands you aren't interested in before writing the modified VRT out to new file. See example below:
from osgeo import gdal
from lxml import etree
import os,sys
def read_vsimem(fn):
'''Read GDAL vsimem files'''
vsifile = gdal.VSIFOpenL(fn,'r')
gdal.VSIFSeekL(vsifile, 0, 2)
vsileng = gdal.VSIFTellL(vsifile)
gdal.VSIFSeekL(vsifile, 0, 0)
return gdal.VSIFReadL(1, vsileng, vsifile)
def write_vsimem(fn,data):
'''Write GDAL vsimem files'''
vsifile = gdal.VSIFOpenL(fn,'w')
size = len(data)
gdal.VSIFWriteL(data, 1, size, vsifile)
return gdal.VSIFCloseL(vsifile)
infn='test.tif'
outfn='test2.tif'
#List of bands to retain, output bands will be reordered to match
bands=[4,3,2]
ds = gdal.Open(os.path.abspath(infn)) #Ensure path is absolute not relative path
vrtdrv=gdal.GetDriverByName('VRT')
tifdrv=gdal.GetDriverByName('GTIFF')
#Create an in memory VRT copy of the input raster
vfn='/vsimem/test.vrt'
vrtds=vrtdrv.CreateCopy(vfn,ds)
#Read the XML from memory and parse it
#Could also write the VRT to a temp file on disk instead of /vsimem
#and used etree.parse(vrtfilepath) instead
#i.e.
#vfn='some/path/on/disk/blah.vrt'
#vrtds=vrtdrv.CreateCopy(vfn,ds)
#tree=etree.parse(vfn)
#os.unlink(vfn)
vrtds=vrtdrv.CreateCopy(vfn,ds)
vrtxml=read_vsimem(vfn)
tree=etree.fromstring(vrtxml)
#Ignore bands not in band list
#And put bands we want to keep in the correct order
for band in tree.findall('VRTRasterBand'):
try:
i=bands.index(int(band.attrib['band']))
band.attrib['band']=str(i+1)
bands[i]=band
except ValueError:pass
finally:
tree.remove(band)
for band in bands:
tree.append(band)
#Get the XML string from the tree
vrtxml=etree.tostring(tree, pretty_print=True)
#Write it to the in memory file
#Could also just write to a file on disk
write_vsimem(vfn,vrtxml)
#Open the VRT containing only the selected bands
vrtds=gdal.Open(vfn)
#Write to a new raster
tifds=tifdrv.CreateCopy(outfn,vrtds, 0, ['COMPRESS=JPEG','PHOTOMETRIC=YCBCR','TFW=YES'])
#Close dataset to ensure cache is properly flushed
del tifds
Best Answer
The answer of @rcoup only worked for me, if modify it as follows:
Otherwise, the file is not written to disk.