PostGIS  2.1.10dev-r@@SVN_REVISION@@
rtrowdump.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 # $Id$
4 #
5 # Brute-force dump of single row from WKT Raster table as GeoTIFF.
6 # This utility is handy for debugging purposes.
7 #
8 # WARNING: Tha main purpose of this program is to test and
9 # debug WKT Raster implementation. It is NOT supposed to be an
10 # efficient performance killer, by no means.
11 #
12 ###############################################################################
13 # Copyright (C) 2009 Mateusz Loskot <mateusz@loskot.net>
14 #
15 # This program is free software; you can redistribute it and/or
16 # modify it under the terms of the GNU General Public License
17 # as published by the Free Software Foundation; either version 2
18 # of the License, or (at your option) any later version.
19 #
20 # This program is distributed in the hope that it will be useful,
21 # but WITHOUT ANY WARRANTY; without even the implied warranty of
22 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
23 # GNU General Public License for more details.
24 #
25 # You should have received a copy of the GNU General Public License
26 # along with this program; if not, write to the Free Software Foundation,
27 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
28 #
29 ###############################################################################
30 import rtreader
31 import numpy
32 import osgeo.gdalconst
33 from osgeo import gdal
34 from optparse import OptionParser
35 import sys
36 
37 def logit(msg):
38  if VERBOSE is True:
39  sys.stderr.write("LOG - " + str(msg) + "\n")
40 
41 def pt2gdt(pt):
42  """Translate WKT Raster pixel type to GDAL type"""
43  pixtypes = {
44  '8BUI' : osgeo.gdalconst.GDT_Byte,
45  '16BSI' : osgeo.gdalconst.GDT_Int16,
46  '16BUI' : osgeo.gdalconst.GDT_UInt16,
47  '32BSI' : osgeo.gdalconst.GDT_Int32,
48  '32BUI' : osgeo.gdalconst.GDT_UInt32,
49  '32BF' : osgeo.gdalconst.GDT_Float32,
50  '64BF' : osgeo.gdalconst.GDT_Float64
51  }
52  return pixtypes.get(pt, 'UNKNOWN')
53 
54 def pt2numpy(pt):
55  """Translate WKT Raster pixel type to NumPy data type"""
56  numpytypes = {
57  '8BUI' : numpy.uint8,
58  '16BSI' : numpy.int16,
59  '16BUI' : numpy.uint16,
60  '32BSI' : numpy.int32,
61  '32BUI' : numpy.uint32,
62  '32BF' : numpy.float32,
63  '64BF' : numpy.float64
64  }
65  return numpytypes.get(pt, numpy.uint8)
66 
67 ###############################################################################
68 try:
69 
70  prs = OptionParser(version="%prog $Revision: 4037 $",
71  usage="%prog -d <DB> -t <TABLE> [-c <COLUMN>]",
72  description="Brute-force dump of single row from WKT Raster table as GeoTIF")
73  prs.add_option("-d", "--db", dest="db", action="store", default=None,
74  help="PostgreSQL database connection string, required")
75  prs.add_option("-t", "--table", dest="table", action="store", default=None,
76  help="table with raster column [<schema>.]<table>, required")
77  prs.add_option("-c", "--column", dest="column", action="store", default="rast",
78  help="raster column, optional, default=rast")
79  prs.add_option("-w", "--where", dest="where", action="store", default="",
80  help="SQL WHERE clause to filter record")
81  prs.add_option("-o", "--output", dest="output", action="store", default=None,
82  help="GeoTIFF output file for pixel data read from WKT Raster table")
83  prs.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
84  help="be excessively verbose and useful for debugging")
85 
86  (opts, args) = prs.parse_args()
87 
88  if opts.db is None:
89  prs.error("use -d option to specify database connection string")
90  if opts.table is None:
91  prs.error("use -t option to specify raster table")
92  if opts.column is None:
93  prs.error("use -c option to specify raster column in raster table")
94  if opts.output is None:
95  prs.error("use -o option to specify raster output file")
96 
97  global VERBOSE
98  VERBOSE = opts.verbose
99 
100  rt = rtreader.RasterReader(opts.db, opts.table, opts.column, opts.where)
101  if VERBOSE is True:
102  rt.logging = True
103 
104  logit("Connected to %s" % opts.db)
105  logit("Source WKT raster:")
106  logit("\trow=%s" % opts.where)
107  logit("\twidth=%d, height=%d, bands=%d, pixel types=%s" \
108  %(rt.width, rt.height, rt.num_bands, str(rt.pixel_types)))
109  logit("Target GeoTIFF: %s" % opts.output)
110 
111  out_format = "GTiff"
112  out_driver = gdal.GetDriverByName(out_format)
113  out_data_type = pt2gdt(rt.pixel_types[0])
114  out_ds = out_driver.Create(opts.output, rt.width, rt.height, rt.num_bands, out_data_type)
115 
116 
117  for b in range(1, rt.num_bands +1):
118  logit("--- BAND %d ---------------------------------" % b)
119 
120  ### Be careful!!
121  ### Zeros function's input parameter can be a (height x width) array,
122  ### not (width x height): http://docs.scipy.org/doc/numpy/reference/generated/numpy.zeros.html?highlight=zeros#numpy.zeros
123  raster = numpy.zeros((rt.height, rt.width), pt2numpy(out_data_type))
124  for width_index in range(0, rt.width):
125  for height_index in range(0, rt.height):
126  pixel = rt.get_value(b, width_index + 1, height_index + 1)
127  raster[height_index, width_index] = pixel
128 
129  logit(str(raster))
130 
131  band = out_ds.GetRasterBand(b)
132  assert band is not None
133  band.WriteArray(raster)
134 
135 except rtreader.RasterError, e:
136  print "ERROR - ", e
def pt2numpy(pt)
Definition: rtrowdump.py:54
def logit(msg)
Definition: rtrowdump.py:37
def pt2gdt(pt)
Definition: rtrowdump.py:41
RASTER driver (read-only)
Definition: rtreader.py:38