1

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?

2 Answers 2

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
1
  • Maybe that's correct for images/rasters that are not clipped. But for my VRT cutline it doesn't work. (gt[0],gt[3]) is not the origin of the original image, it's the origin of the clipped image. And the cutline seems to have the original image coordinates.
    – christoph
    Feb 11, 2016 at 9:45
0

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.