PostGIS  2.1.10dev-r@@SVN_REVISION@@
window.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 # $Id: window.py 12060 2013-10-28 19:44:03Z dustymugs $
4 #
5 # Calculates coordinates of window corners of given raster dataset.
6 # It's just a simple helper for testing and debugging WKT Raster.
7 #
8 ##############################################################################
9 # Copyright (C) 2009 Mateusz Loskot <mateusz@loskot.net>
10 #
11 # This program is free software; you can redistribute it and/or
12 # modify it under the terms of the GNU General Public License
13 # as published by the Free Software Foundation; either version 2
14 # of the License, or (at your option) any later version.
15 #
16 # This program is distributed in the hope that it will be useful,
17 # but WITHOUT ANY WARRANTY; without even the implied warranty of
18 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
19 # GNU General Public License for more details.
20 #
21 # You should have received a copy of the GNU General Public License
22 # along with this program; if not, write to the Free Software Foundation,
23 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
24 #
25 ##############################################################################
26 from osgeo import gdal
27 from osgeo import osr
28 import osgeo.gdalconst as gdalc
29 import sys
30 
31 if len(sys.argv) != 6:
32  print "Usage: window.py <raster> <x> <y> <xsize> <ysize>"
33  print "\traster - GDAL supported dataset"
34  print "\tx - column - 1..N where N is raster X dimension"
35  print "\ty - row - 1..N where N is raster Y dimension"
36  print "\txsize - x-dimension of requested window (xsize <= xsize of raster - x)"
37  print "\tysize - y-dimension of requested window (ysize <= ysize of raster - y)"
38  sys.exit(0)
39 
41  if gt[0] == 0.0 and gt[1] == 1.0 and gt[3] == 0.0 and gt[5] == 1.0:
42  return False
43  else:
44  return True
45 
46 def calculate_corner(gt, x, y):
47  if is_georeferenced(gt):
48  xgeo = gt[0] + gt[1] * x + gt[2] * y
49  ygeo = gt[3] + gt[4] * x + gt[5] * y
50  else:
51  xgeo = x
52  xgeo = y
53  return (xgeo, ygeo)
54 
55 infile = sys.argv[1]
56 inxoff = int(sys.argv[2])
57 inyoff = int(sys.argv[3])
58 inxsize = int(sys.argv[4])
59 inysize = int(sys.argv[5])
60 print "=== INPUT ==="
61 print "File: %s" % infile
62 print "Window:"
63 print "- upper-left: %d x %d" % (inxoff, inyoff)
64 print "- dimensions: %d x %d" % (inxsize, inysize)
65 
66 ds = gdal.Open(infile, gdalc.GA_ReadOnly);
67 if ds is None:
68  sys.exit('Cannot open input file: ' + str(infile))
69 
70 xsize = ds.RasterXSize
71 ysize = ds.RasterYSize
72 print "=== RASTER ==="
73 print "- dimensions: %d x %d" % (xsize, ysize)
74 
75 if inxsize > xsize or inysize > ysize or inxoff > xsize or inyoff > ysize:
76  print "Invalid size of input window"
77  sys.exit(1)
78 
79 gt = ds.GetGeoTransform()
80 res = ( gt[1], gt[5] ) # X/Y pixel resolution
81 ulp = ( gt[0], gt[3] ) # X/Y upper left pixel corner
82 rot = ( gt[2], gt[4] ) # X-/Y-axis rotation
83 
84 if is_georeferenced(gt):
85  print "- pixel size:", res
86  print "- upper left:", ulp
87  print "- rotation :", rot
88 else:
89  print "No georeferencing is available"
90  sys.exit(1)
91 
92 print "=== WINDOW ==="
93 print "- upper-left :", calculate_corner(gt, inxoff, inyoff)
94 print "- lower-left :", calculate_corner(gt, inxoff, ysize)
95 print "- upper-right:", calculate_corner(gt, xsize, inyoff)
96 print "- lower-right:", calculate_corner(gt, xsize, ysize)
97 print "- center :", calculate_corner(gt, xsize/2, ysize/2)
def calculate_corner(gt, x, y)
Definition: window.py:46
def is_georeferenced(gt)
Definition: window.py:40