mirror of
https://git.osgeo.org/gitea/postgis/postgis
synced 2024-10-25 01:22:47 +00:00
dc8b342f29
- Modified PostGIS 'configure.ac' script to add support for new raster type (driven by '--with-raster' configure option. - Added raster build options to 'GNUMakefile' script. - PostGIS Python scripts moved to a 'python' subdirectory inside 'raster/scripts' directory. - References to "WKT" deleted from source code. git-svn-id: http://svn.osgeo.org/postgis/trunk@6019 b70326c6-7e19-0410-871a-916f4a2858ee
86 lines
2.7 KiB
Python
Executable file
86 lines
2.7 KiB
Python
Executable file
#! /usr/bin/env python
|
|
#
|
|
# $Id$
|
|
#
|
|
# Calculates coordinates of window corners of given raster dataset.
|
|
# It's just a simple helper for testing and debugging WKT Raster.
|
|
#
|
|
##############################################################################
|
|
# Copyright (C) 2009 Mateusz Loskot <mateusz@loskot.net>
|
|
#
|
|
# This is free software; you can redistribute and/or modify it under
|
|
# the terms of the GNU General Public Licence. See the COPYING file.
|
|
##############################################################################
|
|
from osgeo import gdal
|
|
from osgeo import osr
|
|
import osgeo.gdalconst as gdalc
|
|
import sys
|
|
|
|
if len(sys.argv) != 6:
|
|
print "Usage: window.py <raster> <x> <y> <xsize> <ysize>"
|
|
print "\traster - GDAL supported dataset"
|
|
print "\tx - column - 1..N where N is raster X dimension"
|
|
print "\ty - row - 1..N where N is raster Y dimension"
|
|
print "\txsize - x-dimension of requested window (xsize <= xsize of raster - x)"
|
|
print "\tysize - y-dimension of requested window (ysize <= ysize of raster - y)"
|
|
sys.exit(0)
|
|
|
|
def is_georeferenced(gt):
|
|
if gt[0] == 0.0 and gt[1] == 1.0 and gt[3] == 0.0 and gt[5] == 1.0:
|
|
return False
|
|
else:
|
|
return True
|
|
|
|
def calculate_corner(gt, x, y):
|
|
if is_georeferenced(gt):
|
|
xgeo = gt[0] + gt[1] * x + gt[2] * y
|
|
ygeo = gt[3] + gt[4] * x + gt[5] * y
|
|
else:
|
|
xgeo = x
|
|
xgeo = y
|
|
return (xgeo, ygeo)
|
|
|
|
infile = sys.argv[1]
|
|
inxoff = int(sys.argv[2])
|
|
inyoff = int(sys.argv[3])
|
|
inxsize = int(sys.argv[4])
|
|
inysize = int(sys.argv[5])
|
|
print "=== INPUT ==="
|
|
print "File: %s" % infile
|
|
print "Window:"
|
|
print "- upper-left: %d x %d" % (inxoff, inyoff)
|
|
print "- dimensions: %d x %d" % (inxsize, inysize)
|
|
|
|
ds = gdal.Open(infile, gdalc.GA_ReadOnly);
|
|
if ds is None:
|
|
sys.exit('Cannot open input file: ' + str(infile))
|
|
|
|
xsize = ds.RasterXSize
|
|
ysize = ds.RasterYSize
|
|
print "=== RASTER ==="
|
|
print "- dimensions: %d x %d" % (xsize, ysize)
|
|
|
|
if inxsize > xsize or inysize > ysize or inxoff > xsize or inyoff > ysize:
|
|
print "Invalid size of input window"
|
|
sys.exit(1)
|
|
|
|
gt = ds.GetGeoTransform()
|
|
res = ( gt[1], gt[5] ) # X/Y pixel resolution
|
|
ulp = ( gt[0], gt[3] ) # X/Y upper left pixel corner
|
|
rot = ( gt[2], gt[4] ) # X-/Y-axis rotation
|
|
|
|
if is_georeferenced(gt):
|
|
print "- pixel size:", res
|
|
print "- upper left:", ulp
|
|
print "- rotation :", rot
|
|
else:
|
|
print "No georeferencing is available"
|
|
sys.exit(1)
|
|
|
|
print "=== WINDOW ==="
|
|
print "- upper-left :", calculate_corner(gt, inxoff, inyoff)
|
|
print "- lower-left :", calculate_corner(gt, inxoff, ysize)
|
|
print "- upper-right:", calculate_corner(gt, xsize, inyoff)
|
|
print "- lower-right:", calculate_corner(gt, xsize, ysize)
|
|
print "- center :", calculate_corner(gt, xsize/2, ysize/2)
|