24 from osgeo
import gdal
26 import osgeo.gdalconst
as gdalc
29 if len(sys.argv) != 6:
30 print "Usage: window.py <raster> <x> <y> <xsize> <ysize>"
31 print "\traster - GDAL supported dataset"
32 print "\tx - column - 1..N where N is raster X dimension"
33 print "\ty - row - 1..N where N is raster Y dimension"
34 print "\txsize - x-dimension of requested window (xsize <= xsize of raster - x)"
35 print "\tysize - y-dimension of requested window (ysize <= ysize of raster - y)"
39 if gt[0] == 0.0
and gt[1] == 1.0
and gt[3] == 0.0
and gt[5] == 1.0:
46 xgeo = gt[0] + gt[1] * x + gt[2] * y
47 ygeo = gt[3] + gt[4] * x + gt[5] * y
54 inxoff = int(sys.argv[2])
55 inyoff = int(sys.argv[3])
56 inxsize = int(sys.argv[4])
57 inysize = int(sys.argv[5])
59 print "File: %s" % infile
61 print "- upper-left: %d x %d" % (inxoff, inyoff)
62 print "- dimensions: %d x %d" % (inxsize, inysize)
64 ds = gdal.Open(infile, gdalc.GA_ReadOnly);
66 sys.exit(
'Cannot open input file: ' + str(infile))
68 xsize = ds.RasterXSize
69 ysize = ds.RasterYSize
70 print "=== RASTER ==="
71 print "- dimensions: %d x %d" % (xsize, ysize)
73 if inxsize > xsize
or inysize > ysize
or inxoff > xsize
or inyoff > ysize:
74 print "Invalid size of input window"
77 gt = ds.GetGeoTransform()
78 res = ( gt[1], gt[5] )
79 ulp = ( gt[0], gt[3] )
80 rot = ( gt[2], gt[4] )
83 print "- pixel size:", res
84 print "- upper left:", ulp
85 print "- rotation :", rot
87 print "No georeferencing is available"
90 print "=== WINDOW ==="
def calculate_corner(gt, x, y)