Quite a lot of time has passed since I asked this question. I ran into this problem again this week and found a small Python program in my toolbox that I would like to share. The Python program `getVRTcutline.py` prints the cutline geometry as WKT string in EPSG:3857, the VRT filename, the image filename, a hasCutline flag, the number of bands, the nodata band and a srcExists flag. If we export the output into a CSV file and add the following tab delimited header, we can copy&paste the file into QGIS map canvas. CSV Header: ``` wkt vrt source hascutline bands nodata srcexists ``` OSGeo4W Command line: `python getVRTcutline.py <vrtfile>` getVRTcutline.py: ```python #!/usr/bin/env python ############################################################################### from osgeo import gdal,ogr from osgeo import osr from xml.dom.minidom import parse import xml.dom.minidom import re import os import sys # ============================================================================= def Usage(): print('') print('Read cutline from vrt input file') print('') print('Usage: getVRTcutline.py vrtfile') print('') sys.exit( 1 ) # ============================================================================= infile = None pixel = None line = None # ============================================================================= # Parse command line arguments. # ============================================================================= i = 1 while i < len(sys.argv): arg = sys.argv[i] if infile is None: infile = arg else: Usage() i = i + 1 if infile is None: Usage() # Open XML document using minidom parser DOMTree = xml.dom.minidom.parse(infile) vrtDataset = DOMTree.documentElement hasCutline = "0" gcplist = [] gcps = vrtDataset.getElementsByTagName("GCP") if len(gcps) > 0: for gcp in gcps: gcplist.append([float(gcp.getAttribute("Pixel")),float(gcp.getAttribute("Line")),float(gcp.getAttribute("X")),float(gcp.getAttribute("Y"))]) # Open input dataset indataset = gdal.Open( infile, gdal.GA_ReadOnly ) sourceSrsWkt = vrtDataset.getElementsByTagName("SourceSRS") if len(sourceSrsWkt) > 0: sourceSrsWkt = sourceSrsWkt[0].firstChild.data else: sourceSrsWkt = vrtDataset.getElementsByTagName("SRS")[0].firstChild.data sourceSRS = osr.SpatialReference() sourceSRS.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) sourceSRS.ImportFromWkt(sourceSrsWkt) targetSRS = osr.SpatialReference() targetSRS.SetAxisMappingStrategy(osr.OAMS_TRADITIONAL_GIS_ORDER) targetSRS.ImportFromEPSG(3857) if len(gcplist)>0: gcps=[gdal.GCP(gcp[2],gcp[3],0,gcp[0],gcp[1],'','') for gcp in gcplist] geomatrix = gdal.GCPsToGeoTransform(gcps) else: srcgeotransform = vrtDataset.getElementsByTagName("SrcGeoTransform") if len(srcgeotransform) > 0: srcgeotransform = '('+srcgeotransform[0].firstChild.data+')' geomatrix =eval(srcgeotransform) cutline = vrtDataset.getElementsByTagName("Cutline") if len(cutline) > 0: hasCutline = "1" cutline = cutline[0].firstChild.data cutline = re.sub(r'^[^\(]+[\(\s]+([^\)]+)[\)\s]+',r'POLYGON((\1))',cutline) else: xsize = vrtDataset.getAttribute("rasterXSize") ysize = vrtDataset.getAttribute("rasterYSize") cutline = 'POLYGON((0 0,'+xsize+' 0,'+xsize+' '+ysize+',0 '+ysize+',0 0))' polygon = ogr.CreateGeometryFromWkt(cutline) geom = polygon.GetGeometryRef(0) ring = ogr.Geometry(ogr.wkbLinearRing) for i in range(0,geom.GetPointCount()): if i == geom.GetPointCount(): i = 0 pixel = geom.GetPoint(i)[0] line = geom.GetPoint(i)[1] X = geomatrix[0] + geomatrix[1] * pixel + geomatrix[2] * line Y = geomatrix[3] + geomatrix[4] * pixel + geomatrix[5] * line # Shift to the center of the pixel X += geomatrix[1] / 2.0 Y += geomatrix[5] / 2.0 ring.AddPoint(X,Y) poly = ogr.Geometry(ogr.wkbPolygon) poly.AddGeometry(ring) srcDataset = vrtDataset.getElementsByTagName("SourceDataset")[0].firstChild.data transform = osr.CoordinateTransformation(sourceSRS, targetSRS) poly = ogr.CreateGeometryFromWkt(poly.ExportToWkt()) poly.Transform(transform) rasterBands = str(len(vrtDataset.getElementsByTagName("VRTRasterBand"))) nodata = vrtDataset.getElementsByTagName("NoDataValue") if len(nodata) > 0: nodata = str(nodata[0].firstChild.data) else: nodata = "" colorinterps = vrtDataset.getElementsByTagName("ColorInterp") for colorint in colorinterps: method = colorint.firstChild.data if method == "Alpha": nodata = method break srcExists = "yes" if not os.path.isfile(srcDataset): srcExists = "no" print(poly.ExportToWkt()+"\t"+infile+"\t"+srcDataset+"\t"+hasCutline+"\t"+rasterBands+"\t"+nodata+"\t"+srcExists) ``` I use a Windows batch routine to extract the cutlines of all VRT files inside a directory: ``` for %%i in (c:\Data\VRT\*.vrt) do python getVRTcutline.py "%%i" >>vrtlist.csv ``` ```