Extract a cutline from VRT file in world coordinates?

Question

I'm trying to extract the world coordinates of a cutline (clipping polygon) in a GDAL VRT file. Unfortunately, the cutline is stored in image coordinates (pixels), and i could not find the right method to transform it into world coordinates.

The following python functions didn't help:

def Pixel2world(geoMatrix, x, y):
  ulX = geoMatrix[0]
  ulY = geoMatrix[3]
  xDist = geoMatrix[1]
  yDist = geoMatrix[5]
  coorX = (ulX + (x * xDist))
  coorY = (ulY + (y * yDist))
  return (coorX, coorY)

def Pixel2coord(geoMatrix, x, y):
  """Returns global coordinates from pixel x, y coords"""
  xoff, a, b, yoff, d, e = geoMatrix
  xp = a * x + b * y + xoff
  yp = d * x + e * y + yoff
  return(xp, yp)

raster=r'virtual_raster.vrt'
ds=gdal.Open(raster)

gt=ds.GetGeoTransform()
print Pixel2coord(gt,1000,1000) #=> Returns wrong coordinates
print Pixel2world(gt,1000,1000) #=> Returns wrong coordinates

BTW: due to a coordinate transformation inside the VRT from EPSG:31255 to EPSG:3857, the resulting image is shifted and rotated. And (geoMatrix[0],geoMatrix[3]) does contain the upper left corner of the cutted image and not the uncutted one.

Has anyone successfully extracted a VRT cutline in world coordinates?

Answer 1

The image/raster coordinates of a point x,y (world coordinates) are

gt = (263104.72544800001, 10.002079999999999, 0.0, 155223.647811, 0.0, -10.002079999999999)
x,y  = (263220.5,155110.6)
rasterx = (x - gt[0]) / gt[1]
rastery = (y - gt[3]) / gt[5]
print rasterx,rastery
11.5750475901 11.3024301945

This allows to extract the value of the pixel at coordinate x, y (a = number of band 1,2, 3)

layer.GetRasterBand(a).ReadAsArray(int(rasterx),int(rastery), 1, 1)

Therefore, if you have the image coordinates (rasterx, rastery)

x = rasterx * gt[1] + gt[0]
y = rastery * gt[5] + gt[3]
print x, y
263220.5 155110.6

Answer 2

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:

#!/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