PostGIS  2.5.7dev-r@@SVN_REVISION@@
window.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 #
4 # Calculates coordinates of window corners of given raster dataset.
5 # It's just a simple helper for testing and debugging WKT Raster.
6 #
7 
24 from osgeo import gdal
25 from osgeo import osr
26 import osgeo.gdalconst as gdalc
27 import sys
28 
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)"
36  sys.exit(0)
37 
39  if gt[0] == 0.0 and gt[1] == 1.0 and gt[3] == 0.0 and gt[5] == 1.0:
40  return False
41  else:
42  return True
43 
44 def calculate_corner(gt, x, y):
45  if is_georeferenced(gt):
46  xgeo = gt[0] + gt[1] * x + gt[2] * y
47  ygeo = gt[3] + gt[4] * x + gt[5] * y
48  else:
49  xgeo = x
50  xgeo = y
51  return (xgeo, ygeo)
52 
53 infile = sys.argv[1]
54 inxoff = int(sys.argv[2])
55 inyoff = int(sys.argv[3])
56 inxsize = int(sys.argv[4])
57 inysize = int(sys.argv[5])
58 print "=== INPUT ==="
59 print "File: %s" % infile
60 print "Window:"
61 print "- upper-left: %d x %d" % (inxoff, inyoff)
62 print "- dimensions: %d x %d" % (inxsize, inysize)
63 
64 ds = gdal.Open(infile, gdalc.GA_ReadOnly);
65 if ds is None:
66  sys.exit('Cannot open input file: ' + str(infile))
67 
68 xsize = ds.RasterXSize
69 ysize = ds.RasterYSize
70 print "=== RASTER ==="
71 print "- dimensions: %d x %d" % (xsize, ysize)
72 
73 if inxsize > xsize or inysize > ysize or inxoff > xsize or inyoff > ysize:
74  print "Invalid size of input window"
75  sys.exit(1)
76 
77 gt = ds.GetGeoTransform()
78 res = ( gt[1], gt[5] ) # X/Y pixel resolution
79 ulp = ( gt[0], gt[3] ) # X/Y upper left pixel corner
80 rot = ( gt[2], gt[4] ) # X-/Y-axis rotation
81 
82 if is_georeferenced(gt):
83  print "- pixel size:", res
84  print "- upper left:", ulp
85  print "- rotation :", rot
86 else:
87  print "No georeferencing is available"
88  sys.exit(1)
89 
90 print "=== WINDOW ==="
91 print "- upper-left :", calculate_corner(gt, inxoff, inyoff)
92 print "- lower-left :", calculate_corner(gt, inxoff, ysize)
93 print "- upper-right:", calculate_corner(gt, xsize, inyoff)
94 print "- lower-right:", calculate_corner(gt, xsize, ysize)
95 print "- center :", calculate_corner(gt, xsize/2, ysize/2)
def is_georeferenced(gt)
Definition: window.py:38
def calculate_corner(gt, x, y)
Definition: window.py:44