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
Command line: python getVRTcutline.py <vrtfile>
getVRTcutline.py:
#!/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 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