PostGIS  2.2.8dev-r@@SVN_REVISION@@
rtrowdump.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 #
4 # Brute-force dump of single row from WKT Raster table as GeoTIFF.
5 # This utility is handy for debugging purposes.
6 #
7 # WARNING: Tha main purpose of this program is to test and
8 # debug WKT Raster implementation. It is NOT supposed to be an
9 # efficient performance killer, by no means.
10 #
11 
28 import rtreader
29 import numpy
30 import osgeo.gdalconst
31 from osgeo import gdal
32 from optparse import OptionParser
33 import sys
34 
35 def logit(msg):
36  if VERBOSE is True:
37  sys.stderr.write("LOG - " + str(msg) + "\n")
38 
39 def pt2gdt(pt):
40  """Translate WKT Raster pixel type to GDAL type"""
41  pixtypes = {
42  '8BUI' : osgeo.gdalconst.GDT_Byte,
43  '16BSI' : osgeo.gdalconst.GDT_Int16,
44  '16BUI' : osgeo.gdalconst.GDT_UInt16,
45  '32BSI' : osgeo.gdalconst.GDT_Int32,
46  '32BUI' : osgeo.gdalconst.GDT_UInt32,
47  '32BF' : osgeo.gdalconst.GDT_Float32,
48  '64BF' : osgeo.gdalconst.GDT_Float64
49  }
50  return pixtypes.get(pt, 'UNKNOWN')
51 
52 def pt2numpy(pt):
53  """Translate WKT Raster pixel type to NumPy data type"""
54  numpytypes = {
55  '8BUI' : numpy.uint8,
56  '16BSI' : numpy.int16,
57  '16BUI' : numpy.uint16,
58  '32BSI' : numpy.int32,
59  '32BUI' : numpy.uint32,
60  '32BF' : numpy.float32,
61  '64BF' : numpy.float64
62  }
63  return numpytypes.get(pt, numpy.uint8)
64 
65 
66 try:
67 
68  prs = OptionParser(version="%prog $Revision: 4037 $",
69  usage="%prog -d <DB> -t <TABLE> [-c <COLUMN>]",
70  description="Brute-force dump of single row from WKT Raster table as GeoTIF")
71  prs.add_option("-d", "--db", dest="db", action="store", default=None,
72  help="PostgreSQL database connection string, required")
73  prs.add_option("-t", "--table", dest="table", action="store", default=None,
74  help="table with raster column [<schema>.]<table>, required")
75  prs.add_option("-c", "--column", dest="column", action="store", default="rast",
76  help="raster column, optional, default=rast")
77  prs.add_option("-w", "--where", dest="where", action="store", default="",
78  help="SQL WHERE clause to filter record")
79  prs.add_option("-o", "--output", dest="output", action="store", default=None,
80  help="GeoTIFF output file for pixel data read from WKT Raster table")
81  prs.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
82  help="be excessively verbose and useful for debugging")
83 
84  (opts, args) = prs.parse_args()
85 
86  if opts.db is None:
87  prs.error("use -d option to specify database connection string")
88  if opts.table is None:
89  prs.error("use -t option to specify raster table")
90  if opts.column is None:
91  prs.error("use -c option to specify raster column in raster table")
92  if opts.output is None:
93  prs.error("use -o option to specify raster output file")
94 
95  global VERBOSE
96  VERBOSE = opts.verbose
97 
98  rt = rtreader.RasterReader(opts.db, opts.table, opts.column, opts.where)
99  if VERBOSE is True:
100  rt.logging = True
101 
102  logit("Connected to %s" % opts.db)
103  logit("Source WKT raster:")
104  logit("\trow=%s" % opts.where)
105  logit("\twidth=%d, height=%d, bands=%d, pixel types=%s" \
106  %(rt.width, rt.height, rt.num_bands, str(rt.pixel_types)))
107  logit("Target GeoTIFF: %s" % opts.output)
108 
109  out_format = "GTiff"
110  out_driver = gdal.GetDriverByName(out_format)
111  out_data_type = pt2gdt(rt.pixel_types[0])
112  out_ds = out_driver.Create(opts.output, rt.width, rt.height, rt.num_bands, out_data_type)
113 
114 
115  for b in range(1, rt.num_bands +1):
116  logit("--- BAND %d ---------------------------------" % b)
117 
118 
121  raster = numpy.zeros((rt.height, rt.width), pt2numpy(out_data_type))
122  for width_index in range(0, rt.width):
123  for height_index in range(0, rt.height):
124  pixel = rt.get_value(b, width_index + 1, height_index + 1)
125  raster[height_index, width_index] = pixel
126 
127  logit(str(raster))
128 
129  band = out_ds.GetRasterBand(b)
130  assert band is not None
131  band.WriteArray(raster)
132 
133 except rtreader.RasterError, e:
134  print "ERROR - ", e
def pt2numpy(pt)
Definition: rtrowdump.py:52
def logit(msg)
Definition: rtrowdump.py:35
def pt2gdt(pt)
Definition: rtrowdump.py:39
RASTER driver (read-only)
Definition: rtreader.py:37