PostGIS  2.1.10dev-r@@SVN_REVISION@@
raster2pgsql.py
Go to the documentation of this file.
1 #! /usr/bin/env python
2 #
3 # $Id: raster2pgsql.py 12060 2013-10-28 19:44:03Z dustymugs $
4 #
5 # This is a simple utility used to dump GDAL dataset into HEX WKB stream.
6 # It's considered as a prototype of raster2pgsql tool planned to develop
7 # in future.
8 # For more details about raster2pgsql tool, see Specification page:
9 # http://trac.osgeo.org/postgis/wiki/WKTRaster
10 #
11 # The script requires Python bindings for GDAL.
12 # Available at http://trac.osgeo.org/gdal/wiki/GdalOgrInPython
13 #
14 ################################################################################
15 # Copyright (C) 2009-2010 Mateusz Loskot <mateusz@loskot.net>
16 # Copyright (C) 2009-2011 Pierre Racine <pierre.racine@sbf.ulaval.ca>
17 # Copyright (C) 2009-2010 Jorge Arevalo <jorge.arevalo@deimos-space.com>
18 #
19 # This program is free software; you can redistribute it and/or
20 # modify it under the terms of the GNU General Public License
21 # as published by the Free Software Foundation; either version 2
22 # of the License, or (at your option) any later version.
23 #
24 # This program is distributed in the hope that it will be useful,
25 # but WITHOUT ANY WARRANTY; without even the implied warranty of
26 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
27 # GNU General Public License for more details.
28 #
29 # You should have received a copy of the GNU General Public License
30 # along with this program; if not, write to the Free Software Foundation,
31 # Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
32 #
33 ################################################################################
34 #
35 from osgeo import gdal
36 from osgeo import osr
37 import osgeo.gdalconst as gdalc
38 from optparse import OptionParser, OptionGroup
39 import binascii
40 import glob
41 import math
42 import numpy
43 import os
44 import sys
45 
46 ################################################################################
47 # CONSTANTS - DO NOT CHANGE
48 
49 # Endianness enumeration
50 NDR = 1 # Little-endian
51 XDR = 0 # Big-endian
52 
53 # Default version of WKTRaster protocol.
54 # WARNING: Currently, this is the only valid value
55 # and option -w, --raster-version is ignored, if specified.
56 g_rt_version = 0
57 
58 # Default format of binary output is little-endian (NDR)
59 # WARNING: Currently, big-endian (XDR) output is not supported
60 # and option -e, --endian is ignored, if specified.
61 g_rt_endian = NDR
62 
63 # Default name of column, overriden with -f, --field option.
64 g_rt_column = 'rast'
65 
66 g_rt_catalog = ''
67 g_rt_schema = 'public'
68 
69 ################################################################################
70 # UTILITIES
71 VERBOSE = False
72 SUMMARY = []
73 
74 def is_nan(x):
75  if sys.hexversion < 0x02060000:
76  return str(float(x)).lower() == 'nan'
77  else:
78  return math.isnan(x)
79 
81  """Collects, parses and validates command line arguments."""
82 
83  prs = OptionParser(version="%prog $Revision: 12060 $")
84 
85  # Mandatory parameters
86  grp0 = OptionGroup(prs, "Source and destination",
87  "*** Mandatory parameters always required ***")
88  grp0.add_option("-r", "--raster", dest="raster", action="append", default=None,
89  help="append raster to list of input files, at least one raster "
90  "file required. You may use wildcards (?,*) for specifying multiple files.")
91  grp0.add_option("-t", "--table", dest="table", action="store", default=None,
92  help="raster destination in form of [<schema>.]<table> or base raster table for overview level>1, required")
93  prs.add_option_group(grp0);
94 
95  # Optional parameters - raster manipulation
96  grp_r = OptionGroup(prs, "Raster processing",
97  "Optional parameters used to manipulate input raster dataset")
98  grp_r.add_option("-s", "--srid", dest="srid", action="store", type="int", default=-1,
99  help="assign output raster with specified SRID")
100  grp_r.add_option("-b", "--band", dest="band", action="store", type="int", default=None,
101  help="specify number of the band to be extracted from raster file")
102  grp_r.add_option("-k", "--block-size", dest="block_size", action="store", default=None,
103  help="cut raster(s) into tiles to be inserted one by table row."
104  "BLOCK_SIZE is expressed as WIDTHxHEIGHT. Incomplete tiles are completed with nodata values")
105  grp_r.add_option("-R", "--register", dest="register", action="store_true", default=False,
106  help="register the raster as a filesystem (out-db) raster")
107  grp_r.add_option("-l", "--overview-level", dest="overview_level", action="store", type="int", default=1,
108  help='create overview tables named as o_<LEVEL>_<RASTER_TABLE> and '
109  'populate with GDAL-provided overviews (regular blocking only)')
110  prs.add_option_group(grp_r);
111 
112  # Optional parameters - database/table manipulation
113  grp_t = OptionGroup(prs, 'Database processing',
114  'Optional parameters used to manipulate database objects')
115  grp_t.add_option('-c', '--create', dest='create_table', action='store_true', default=False,
116  help='create new table and populate it with raster(s), this is the default mode')
117  grp_t.add_option('-a', '--append', dest='append_table', action='store_true', default=False,
118  help='append raster(s) to an existing table')
119  grp_t.add_option("-d", "--drop", dest="drop_table", action="store_true", default=False,
120  help="drop table, create new one and populate it with raster(s)")
121  grp_t.add_option("-f", "--field", dest="column", action="store", default=g_rt_column,
122  help="specify name of destination raster column, default is 'rast'")
123  grp_t.add_option("-F", "--filename", dest="filename", action="store_true", default=False,
124  help="add a column with the name of the file")
125  grp_t.add_option("-I", "--index", dest="index", action="store_true", default=False,
126  help="create a GiST index on the raster column")
127  grp_t.add_option("-M", "--vacuum", dest="vacuum", action="store_true", default=False,
128  help="issue VACUUM command against all generated tables")
129  grp_t.add_option('-V', '--create-raster-overviews', dest='create_raster_overviews_table',
130  action='store_true', default=False,
131  help='create RASTER_OVERVIEWS table used to store overviews metadata')
132  prs.add_option_group(grp_t);
133 
134  # Other optional parameters
135  grp_u = OptionGroup(prs, "Miscellaneous", "Other optional parameters")
136  grp_u.add_option("-e", "--endian", dest="endian", action="store", type="int", default=g_rt_endian,
137  help="control endianness of generated binary output of raster; "
138  "specify 0 for XDR and 1 for NDR (default); "
139  "only NDR output is supported now")
140  grp_u.add_option("-w", "--raster-version", dest="version",
141  action="store", type="int", default=g_rt_version,
142  help="specify version of raster protocol, default is 0; "
143  "only value of zero is supported now")
144  grp_u.add_option("-o", "--output", dest="output", action="store", default=sys.stdout,
145  help="specify output file, otherwise send to stdout")
146  grp_u.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
147  help="verbose mode. Useful for debugging")
148  prs.add_option_group(grp_u);
149 
150  (opts, args) = prs.parse_args()
151 
152  # Validate options
153  if opts.create_table and opts.drop_table and opts.append_table:
154  prs.error("options -c, -a and -d are mutually exclusive")
155  if opts.create_table and opts.drop_table:
156  prs.error("options -c and -d are mutually exclusive")
157  if opts.create_table and opts.append_table:
158  prs.error("options -c and -a are mutually exclusive")
159  if opts.append_table and opts.drop_table:
160  prs.error("options -a and -d are mutually exclusive")
161  if (not opts.create_table and not opts.drop_table and not opts.append_table) or opts.drop_table:
162  opts.create_table = True
163 
164  if opts.raster is None:
165  prs.error("use option -r to specify at least one input raster. Wildcards (?,*) are accepted.")
166 
167  if opts.block_size is not None and len(opts.raster) != 1:
168  prs.error("regular blocking supports single-raster input only")
169 
170  if opts.block_size is not None:
171  if len(opts.block_size.split('x')) != 2 and len(opts.block_size.split('X')) != 2:
172  prs.error("invalid format of block size, expected WIDTHxHEIGHT")
173 
174  if opts.overview_level > 1 and opts.block_size is None:
175  prs.error("regular blocking mode required to enable overviews support (level > 1)")
176 
177  if opts.create_raster_overviews_table and opts.overview_level <= 1:
178  prs.error('create table for RASTER_OVERVIEWS available only if overviews import requested')
179 
180  # XXX: Now, if --band=Nth, then only Nth band of all specified rasters is dumped/imported
181  # This behavior can be changed to support only single-raster input if --band option used.
182  #if opts.band is not None and len(opts.raster) > 1:
183  # prs.error("option -b requires single input raster only, specify -r option only once")
184 
185  if opts.table is None:
186  prs.error("use option -t to specify raster destination table")
187  if len(opts.table.split('.')) > 2:
188  prs.error("invalid format of table name specified with option -t, expected [<schema>.]table")
189 
190  if opts.output is None:
191  prs.error("failed to initialise output file, try to use option -o explicitly")
192 
193  if opts.version is not None:
194  if opts.version != g_rt_version:
195  prs.error("invalid version of WKT Raster protocol specified, only version 0 is supported")
196  else:
197  prs.error("use option -w to specify version of WKT Raster protocol")
198 
199  if opts.endian is not None:
200  if opts.endian != NDR and opts.endian != XDR:
201  prs.error("invalid endianness value, valid ones are 0 for NDR or 1 for XDR")
202  else:
203  prs.error("use option -e to specify endianness of binary output")
204 
205  return (opts, args)
206 
207 
208 def logit(msg):
209  """If verbose mode requested, sends extra progress information to stderr"""
210  if VERBOSE is True:
211  sys.stderr.write(msg)
212 
213 
214 def gdt2pt(gdt):
215  """Translate GDAL data type to WKT Raster pixel type."""
216  pixtypes = {
217  gdalc.GDT_Byte : { 'name': 'PT_8BUI', 'id': 4 },
218  gdalc.GDT_Int16 : { 'name': 'PT_16BSI', 'id': 5 },
219  gdalc.GDT_UInt16 : { 'name': 'PT_16BUI', 'id': 6 },
220  gdalc.GDT_Int32 : { 'name': 'PT_32BSI', 'id': 7 },
221  gdalc.GDT_UInt32 : { 'name': 'PT_32BUI', 'id': 8 },
222  gdalc.GDT_Float32 : { 'name': 'PT_32BF', 'id': 10 },
223  gdalc.GDT_Float64 : { 'name': 'PT_64BF', 'id': 11 }
224  }
225 
226  # XXX: Uncomment these logs to debug types translation
227  #logit('MSG: Input GDAL pixel type: %s (%d)\n' % (gdal.GetDataTypeName(gdt), gdt))
228  #logit('MSG: Output WKTRaster pixel type: %(name)s (%(id)d)\n' % (pixtypes.get(gdt, 13)))
229 
230  return pixtypes.get(gdt, 13)
231 
232 def pt2numpy(pt):
233  """Translate GDAL data type to NumPy data type"""
234  ptnumpy = {
235  gdalc.GDT_Byte : numpy.uint8,
236  gdalc.GDT_Int16 : numpy.int16,
237  gdalc.GDT_UInt16 : numpy.uint16,
238  gdalc.GDT_Int32 : numpy.int32,
239  gdalc.GDT_UInt32 : numpy.uint32,
240  gdalc.GDT_Float32: numpy.float32,
241  gdalc.GDT_Float64: numpy.float64
242  }
243  return ptnumpy.get(pt, numpy.uint8)
244 
245 def pt2fmt(pt):
246  """Returns binary data type specifier for given pixel type."""
247  fmttypes = {
248  4: 'B', # PT_8BUI
249  5: 'h', # PT_16BSI
250  6: 'H', # PT_16BUI
251  7: 'i', # PT_32BSI
252  8: 'I', # PT_32BUI
253  10: 'f', # PT_32BF
254  11: 'd' # PT_64BF
255  }
256  return fmttypes.get(pt, 'x')
257 
258 
259 def fmt2printfmt(fmt):
260  """Returns printf-like formatter for given binary data type sepecifier."""
261  fmttypes = {
262  'B': '%d', # PT_8BUI
263  'h': '%d', # PT_16BSI
264  'H': '%d', # PT_16BUI
265  'i': '%d', # PT_32BSI
266  'I': '%d', # PT_32BUI
267  'f': '%.15f', # PT_32BF
268  'd': '%.15f', # PT_64BF
269  's': '%s'
270  }
271  return fmttypes.get(fmt, 'f')
272 
273 def parse_block_size(options):
274  assert options is not None
275  assert options.block_size is not None
276 
277  wh = options.block_size.split('x')
278  if len(wh) != 2:
279  wh = options.block_size.split('X')
280 
281  assert len(wh) == 2, "invalid format of specified block size"
282  return ( int(wh[0]), int(wh[1]) )
283 
284 ################################################################################
285 # SQL OPERATIONS
286 
287 def quote_sql_value(value):
288  assert value is not None, "None value given"
289 
290  if len(value) > 0 and value[0] != "'" and value[:-1] != "'":
291  sql = "'" + str(value) + "'"
292  else:
293  sql = value
294  return sql
295 
296 def quote_sql_name(name):
297  assert name is not None, "None name given"
298 
299  if name[0] != "\"" and name[:-1] != "\"":
300  sql = "\"" + str(name) + "\""
301  else:
302  sql = name
303  return sql
304 
306  sql = "ARRAY["
307  for v in values:
308  if type(v) == str:
309  sql += quote_sql_value(v) + ","
310  else:
311  sql += str(v) + ','
312  sql = sql[:-1] # Trim comma
313  sql += "]"
314  return sql
315 
316 def make_sql_schema_table_names(schema_table):
317  st = schema_table.split('.')
318  if len(st) == 1:
319  # TODO: Should we warn user that public is used implicitly?
320  st.insert(0, "public")
321  assert len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
322  return (st[0], st[1])
323 
324 def make_sql_full_table_name(schema_table):
325  st = make_sql_schema_table_names(schema_table)
326  table = "\"%s\".\"%s\"" % (st[0], st[1])
327  return table
328 
329 def make_sql_table_name(schema_table):
330  st = schema_table.split('.')
331  assert len(st) == 1 or len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
332  if len(st) == 2:
333  return st[1]
334  return st[0]
335 
337  sql = "DROP TABLE IF EXISTS %s CASCADE;\n" \
338  % make_sql_full_table_name(table)
339  logit("SQL: %s" % sql)
340  return sql
341 
343  st = make_sql_schema_table_names(table)
344 
345  if len(st[0]) == 0:
346  target = "'', '%s'" % st[1]
347  else:
348  target = "'%s', '%s'" % (st[0], st[1])
349  sql = "SELECT DropRasterTable(%s);\n" % target
350  logit("SQL: %s" % sql)
351  return sql
352 
353 
354 def make_sql_create_table(options, table = None, is_overview = False):
355 
356  if table is None:
357  table = options.table
358 
359  if options.filename:
360  sql = "CREATE TABLE %s (rid serial PRIMARY KEY, \"filename\" text);\n" \
361  % (make_sql_full_table_name(table))
362  else:
363  if options.overview_level > 1 and is_overview:
364  sql = "CREATE TABLE %s (rid serial PRIMARY KEY, %s RASTER);\n" \
365  % (make_sql_full_table_name(table), quote_sql_name(options.column))
366  else:
367  sql = "CREATE TABLE %s (rid serial PRIMARY KEY);\n" \
368  % (make_sql_full_table_name(table))
369 
370  logit("SQL: %s" % sql)
371  return sql
372 
373 
374 def make_sql_create_gist(table, column):
375  gist_table = make_sql_table_name(table)
376  target_table = make_sql_full_table_name(table)
377  logit("MSG: Create GiST spatial index on %s\n" % gist_table)
378 
379  sql = "CREATE INDEX \"%s_%s_gist_idx\" ON %s USING GIST (st_convexhull(%s));\n" % \
380  (gist_table, column, target_table, column)
381  logit("SQL: %s" % sql)
382  return sql;
383 
384 
385 def make_sql_addrastercolumn(options, pixeltypes, nodata, pixelsize, blocksize, extent):
386  assert len(pixeltypes) > 0, "No pixel types given"
387  ts = make_sql_schema_table_names(options.table)
388  pt = make_sql_value_array(pixeltypes)
389 
390  nd = 'null'
391  if nodata is not None and len(nodata) > 0:
392  nd = make_sql_value_array(nodata)
393 
394  odb = 'false'
395  if options.register:
396  odb = 'true'
397 
398  rb = 'false'
399  extgeom = 'null'
400  bs = ( 'null', 'null' )
401  # Check if regular blocking mode requested
402  if options.block_size is not None:
403  assert pixelsize is not None, "Pixel size is none, but regular blocking requested"
404  assert blocksize is not None, "Block size is none, but regular blocking requested"
405  assert extent is not None, "Extent is none, but regular blocking requested"
406  assert len(pixelsize) == 2, "Invalid pixel size, two values expected"
407  assert len(blocksize) == 2, "Invalid block size, two values expected"
408  assert len(extent) == 4, "Invalid extent, four coordinates expected"
409  assert len(extent[0]) == len(extent[3]) == 2, "Invalid extent, pair of X and Y pair expected"
410  rb = 'true'
411  bs = ( blocksize[0], blocksize[1] )
412  extgeom = "ST_Envelope(ST_SetSRID('POLYGON((%.15f %.15f,%.15f %.15f,%.15f %.15f,%.15f %.15f,%.15f %.15f))'::geometry, %d))" % \
413  (extent[0][0], extent[0][1], extent[1][0], extent[1][1],
414  extent[2][0], extent[2][1], extent[3][0], extent[3][1],
415  extent[0][0], extent[0][1], options.srid)
416 
417  sql = "SELECT AddRasterColumn('%s','%s','%s',%d, %s, %s, %s, %s, %.15f, %.15f, %s, %s, %s);\n" % \
418  (ts[0], ts[1], options.column, options.srid, pt, odb, rb, nd,
419  pixelsize[0], pixelsize[1], bs[0], bs[1], extgeom)
420 
421  logit("SQL: %s" % sql)
422  return sql
423 
424 def make_sql_insert_raster(table, rast, hexwkb, insert_filename, file):
425  if insert_filename:
426  assert file is not None, "Missing filename, but insert_filename requested"
427  sql = "INSERT INTO %s ( filename, %s ) VALUES ( (\'%s\')::text, (\'%s\')::raster );\n" \
428  % (make_sql_full_table_name(table), rast, file, hexwkb)
429  else:
430  sql = "INSERT INTO %s ( %s ) VALUES ( (\'%s\')::raster );\n" \
431  % (make_sql_full_table_name(table), rast, hexwkb)
432 
433  # NOTE: Usually, no need for such detailed verbosity
434  #logit("SQL: %s" % sql)
435  return sql
436 
438  schema = make_sql_schema_table_names(options.table)[0]
439  table = make_sql_full_table_name(schema + '.raster_overviews')
440  sql = 'CREATE TABLE ' + table + ' ( ' \
441  'o_table_catalog character varying(256) NOT NULL, ' \
442  'o_table_schema character varying(256) NOT NULL, ' \
443  'o_table_name character varying(256) NOT NULL, ' \
444  'o_column character varying(256) NOT NULL, ' \
445  'r_table_catalog character varying(256) NOT NULL, ' \
446  'r_table_schema character varying(256) NOT NULL, ' \
447  'r_table_name character varying(256) NOT NULL, ' \
448  'r_column character varying(256) NOT NULL, ' \
449  'out_db boolean NOT NULL, ' \
450  'overview_factor integer NOT NULL, ' \
451  'CONSTRAINT raster_overviews_pk ' \
452  'PRIMARY KEY (o_table_catalog, o_table_schema, o_table_name, o_column, overview_factor));\n'
453 
454  return sql
455 
456 
457 def make_sql_register_overview(options, ov_table, ov_factor):
458  assert len(ov_table) > 0
459  assert ov_factor > 0
460 
461  catalog = quote_sql_value('')
462  schema = make_sql_schema_table_names(options.table)[0]
463  r_table = make_sql_table_name(options.table)
464 
465  sql = "INSERT INTO public.raster_overviews( " \
466  "o_table_catalog, o_table_schema, o_table_name, o_column, " \
467  "r_table_catalog, r_table_schema, r_table_name, r_column, out_db, overview_factor) " \
468  "VALUES ('%s', '%s', '%s', '%s', '%s', '%s', '%s', '%s', FALSE, %d);\n" % \
469  (catalog, schema, ov_table, options.column, catalog, schema, r_table, options.column, ov_factor)
470 
471  return sql
472 
473 def make_sql_vacuum(table):
474  sql = 'VACUUM ANALYZE ' + make_sql_full_table_name(table) + ';\n'
475  return sql
476 
477 ################################################################################
478 # RASTER OPERATIONS
479 
480 def calculate_overviews(ds, band_from = None, band_to = None):
481  assert ds is not None
482 
483  if band_from is None:
484  band_from = 0
485  if band_to is None:
486  band_to = ds.RasterCount
487 
488  assert band_to <= ds.RasterCount,'Failed checking band_to=%d <= RasterCount=%d' % (band_to,ds.RasterCount)
489  assert band_from <= band_to
490 
491  nov = 0
492  for i in range(band_from, band_to + 1):
493  n = ds.GetRasterBand(i).GetOverviewCount()
494  if 0 == nov:
495  nov = n
496  assert n == nov, 'Number of raster overviews is not the same for all bands'
497 
498  return nov
499 
500 def calculate_overview_factor(ds, overview):
501  assert ds is not None
502 
503 
504  # Assume all bands have same layout of overviews
505  band = ds.GetRasterBand(1)
506  assert band is not None
507  assert overview < band.GetOverviewCount()
508 
509  ov_band = band.GetOverview(overview)
510  assert ov_band is not None
511 
512  ovf = int(0.5 + ds.RasterXSize / float(ov_band.XSize))
513  logit('MSG: Overview factor = %d\n' % ovf)
514 
515  return ovf
516 
517 
518 def collect_pixel_types(ds, band_from, band_to):
519  """Collect pixel types of bands in requested range.
520  Use names of pixel types in format as returned
521  by rt_core function rt_pixtype_name()"""
522 
523  pt =[]
524  for i in range(band_from, band_to):
525  band = ds.GetRasterBand(i)
526  pixel_type = gdt2pt(band.DataType)['name'][3:]
527  pt.append(pixel_type)
528 
529  return pt
530 
531 def collect_nodata_values(ds, band_from, band_to):
532  """Collect nodata values of bands in requested range"""
533 
534  nd = []
535  for i in range(band_from, band_to):
536  band = ds.GetRasterBand(i)
537  nodata = band.GetNoDataValue()
538  if nodata is not None and not is_nan(nodata):
539  nd.append(nodata)
540 
541  return nd
542 
543 def calculate_block_size(ds, band_from, band_to):
544  """Size of natural block reported by GDAL for bands of given dataset"""
545 
546  block_dims = None
547  for i in range(band_from, band_to):
548  band = ds.GetRasterBand(i)
549  assert band is not None, "Cannot access raster band %d" % i
550  dims = band.GetBlockSize()
551 
552  # Assume bands with common block size
553  if i == band_from:
554  block_dims = dims
555 
556  # Validate dimensions of bands block
557  if block_dims != dims:
558  logit("MSG: Block sizes don't match: %s != %s\n" % (str(block_dims), str(dims)))
559 
560  assert block_dims is not None, "Failed to calculate block size"
561  return (int(block_dims[0]), int(block_dims[1]))
562 
563 def calculate_grid_size(raster_size, block_size):
564  """Dimensions of grid made up with blocks of requested size"""
565 
566  # Exact number of grid dimensions
567  nx = float(raster_size[0]) / float(block_size[0])
568  ny = float(raster_size[1]) / float(block_size[1])
569 
570  return ( int(math.ceil(nx)), int(math.ceil(ny)))
571 
572 def calculate_block_pad_size(band, xoff, yoff, block_size):
573  """Calculates number of columns [0] and rows [1] of padding"""
574  assert band is not None
575 
576  xpad = 0
577  ypad= 0
578  block_bound = ( xoff + block_size[0], yoff + block_size[1] )
579 
580  if block_bound[0] > band.XSize:
581  xpad = block_bound[0] - band.XSize
582  if block_bound[1] > band.YSize:
583  ypad = block_bound[1] - band.YSize
584 
585  return (xpad, ypad)
586 
588  assert ds is not None
589  gt = list(ds.GetGeoTransform())
590  return tuple(gt)
591 
592 def calculate_geoxy(gt, xy):
593  """Calculate georeferenced coordinate from given x and y"""
594  assert gt is not None
595  assert xy is not None
596  assert len(xy) == 2
597 
598  xgeo = gt[0] + gt[1] * xy[0] + gt[2] * xy[1];
599  ygeo = gt[3] + gt[4] * xy[0] + gt[5] * xy[1];
600 
601  return (xgeo, ygeo)
602 
603 def calculate_geoxy_level(gt, xy, level):
604 
605  # Update pixel resolution according to overview level
606  newgt = ( gt[0], gt[1] * float(level), gt[2], gt[3], gt[4], gt[5] * float(level) )
607 
608  return calculate_geoxy(newgt, xy)
609 
611  """Calculate georeferenced coordinates of spatial extent of raster dataset"""
612  assert ds is not None
613 
614  # UL, LL, UR, LR
615  dim = ( (0,0),(0,ds.RasterYSize),(ds.RasterXSize,0),(ds.RasterXSize,ds.RasterYSize) )
616 
617  ext = (calculate_geoxy(gt, dim[0]), calculate_geoxy(gt, dim[1]),
618  calculate_geoxy(gt, dim[2]), calculate_geoxy(gt, dim[3]))
619 
620  return ext
621 
622 def check_hex(hex, bytes_size = None):
623  assert hex is not None, "Error: Missing hex string"
624  size = len(hex)
625  assert size > 0, "Error: hex string is empty"
626  assert size % 2 == 0, "Error: invalid size of hex string"
627  if bytes_size is not None:
628  n = int(size / 2)
629  assert n == bytes_size, "Error: invalid number of bytes %d, expected %d" % (n, bytes_size)
630 
631 def dump_block_numpy(pixels):
632  assert pixels.ndim == 2
633  print 'BEGIN BLOCK SCANLINES (numpy): (%d, %d)' %(len(pixels[0]), len(pixels))
634 
635  i = 0
636  for row in range (0, len(pixels)):
637  s = binascii.hexlify(pixels[row])
638  print '%d (%d)\t%s' % (i, (len(s) / 2), s)
639  i = i + 1
640 
641  print 'END BLOCK SCANLINES'
642 
643 def fetch_band_nodata(band, default = 0):
644  assert band is not None
645 
646  nodata = default
647  if band.GetNoDataValue() is not None:
648  nodata = band.GetNoDataValue()
649  else:
650  logit("WARNING: No nodata flagged in raster_columns metadata. "
651  "In serialized raster, nodata bytes will have value of 0.\n")
652  return nodata
653 
654 def wkblify(fmt, data):
655  """Writes raw binary data into HEX-encoded string using binascii module."""
656  import struct
657 
658  # Binary to HEX
659  fmt_little = '<' +fmt
660  hexstr = binascii.hexlify(struct.pack(fmt_little, data)).upper()
661 
662  # String'ify raw value for log
663  valfmt = '\'' + fmt2printfmt(fmt[len(fmt) - 1]) + '\''
664  val = valfmt % data
665  logit('HEX (\'fmt=%s\', bytes=%d, val=%s):\t\t%s\n' \
666  % (fmt, len(hexstr) / 2, str(val), hexstr))
667 
668  return hexstr
669 
670 def wkblify_raster_header(options, ds, level, ulp, xsize = None, ysize = None):
671  """Writes WKT Raster header based on given GDAL into HEX-encoded WKB."""
672  assert ds is not None, "Error: Missing GDAL dataset"
673  assert level >= 1
674  assert len(ulp) == 2 is not None, "Error: invalid upper-left corner"
675 
676  if xsize is None or ysize is None:
677  assert xsize is None and ysize is None
678  xsize = ds.RasterXSize
679  ysize = ds.RasterYSize
680 
681  # Collect GeoReference information
682  gt = get_gdal_geotransform(ds)
683  ul = calculate_geoxy(gt, (ulp[0], ulp[1]))
684  rt_ip = ( ul[0], ul[1] )
685  rt_skew = ( gt[2], gt[4] )
686  rt_scale = ( gt[1] * level, gt[5] * level )
687 
688  # TODO: Any way to lookup for SRID based on SRS in WKT?
689  #srs = osr.SpatialReference()
690  #srs.ImportFromWkt(ds.GetProjection())
691 
692  # Burn input raster as WKTRaster WKB format
693  hexwkb = ''
694  ### Endiannes
695  hexwkb += wkblify('B', options.endian)
696  ### Version
697  hexwkb += wkblify('H', options.version)
698  ### Number of bands
699  if options.band is not None and options.band > 0:
700  hexwkb += wkblify('H', 1)
701  else:
702  hexwkb += wkblify('H', ds.RasterCount)
703  check_hex(hexwkb, 5)
704  ### Georeference
705  hexwkb += wkblify('d', rt_scale[0])
706  hexwkb += wkblify('d', rt_scale[1])
707  hexwkb += wkblify('d', rt_ip[0])
708  hexwkb += wkblify('d', rt_ip[1])
709  hexwkb += wkblify('d', rt_skew[0])
710  hexwkb += wkblify('d', rt_skew[1])
711  hexwkb += wkblify('i', options.srid)
712  check_hex(hexwkb, 57)
713  ### Number of columns and rows
714  hexwkb += wkblify('H', xsize)
715  hexwkb += wkblify('H', ysize)
716  check_hex(hexwkb, 61)
717 
718  logit("MSG: Georeference: px = %s -> ul = %s \tresolution = %s \trotation = %s\n" \
719  % (str(ulp), str(rt_ip), str(rt_scale), str(rt_skew)))
720  return hexwkb
721 
722 def wkblify_band_header(options, band):
723  """Writes band header into HEX-encoded WKB"""
724  assert band is not None, "Error: Missing GDAL raster band"
725 
726  hexwkb = ""
727 
728  first4bits = 0
729 
730  # If the register option is enabled, set the first bit to 1
731  if options.register:
732  first4bits = 128
733 
734  nodata = band.GetNoDataValue()
735  # If there is no nodata value, set it to 0. Otherwise set the HasNodata bit to 1
736  if nodata is not None:
737  first4bits += 64
738  else:
739  nodata = 0
740 
741  # Encode pixel type
742  pixtype = gdt2pt(band.DataType)['id']
743  hexwkb += wkblify('B', pixtype + first4bits)
744 
745  # Encode nodata value (or Zero, if nodata unavailable)
746  hexwkb += wkblify(pt2fmt(pixtype), nodata)
747 
748  check_hex(hexwkb)
749  return hexwkb
750 
751 def wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size, infile, bandidx):
752  """Writes band of given GDAL dataset into HEX-encoded WKB for WKT Raster output."""
753  assert band is not None, "Error: Missing GDAL raster band"
754 
755  hexwkb = ''
756 
757  if options.register:
758  # Off-db raster
759  # TODO: Do we want to handle options.overview_level? --mloskot
760  # ANSWER:
761  # TODO: Where bandidx and ds come from? --mloskot
762  # ANSWER: Provided by caller method --jorgearevalo
763  hexwkb += wkblify('B', bandidx - 1)
764  filepath = os.path.abspath(infile.replace('\\', '\\\\'))
765  logit('MSG: Out-db raster path=%s\n' % filepath)
766  hexwkb += wkblify(str(len(filepath)) + 's', filepath)
767  hexwkb += wkblify('B', 0)
768  else:
769  # In-db raster
770 
771  # Right most column and bottom most row of blocks have
772  # portions that extend beyond the raster
773  read_padding_size = calculate_block_pad_size(band, xoff, yoff, read_block_size)
774  valid_read_block_size = ( read_block_size[0] - read_padding_size[0],
775  read_block_size[1] - read_padding_size[1] )
776 
777 
778  if read_padding_size[0] > 0 or read_padding_size[1] > 0:
779  target_block_size = (valid_read_block_size[0] / level, valid_read_block_size[1] / level)
780  target_padding_size = (read_padding_size[0] / level, read_padding_size[1] / level)
781  else:
782  target_block_size = block_size
783  target_padding_size = ( 0, 0 )
784 
785  logit('MSG: Normalize read_block=%s for level=%d to valid_read_block=%s with padding=%s\n' % \
786  (read_block_size, level, valid_read_block_size, read_padding_size))
787  logit('MSG: Normalize target_block=%s for level=%d to valid_target_block=%s with padding=%s\n' % \
788  (block_size, level, target_block_size, target_padding_size))
789  logit('MSG: ReadAsArray( %d, %d, %s, %s)\n' % \
790  (xoff, yoff, str(valid_read_block_size), str(target_block_size)))
791 
792  assert valid_read_block_size[0] > 0 and valid_read_block_size[1] > 0
793  assert target_block_size[0] > 0 and target_block_size[1] > 0
794 
795  pixels = band.ReadAsArray(xoff, yoff, valid_read_block_size[0], valid_read_block_size[1],
796  target_block_size[0], target_block_size[1])
797 
798  # XXX: Use for debugging only
799  #dump_block_numpy(pixels)
800 
801  out_pixels = numpy.zeros((block_size[1], block_size[0]), pt2numpy(band.DataType))
802 
803  logit('MSG: Read valid source:\t%d x %d\n' % (len(pixels[0]), len(pixels)))
804  logit('MSG: Write into block:\t%d x %d\n' % (len(out_pixels[0]), len(out_pixels)))
805 
806  if target_padding_size[0] > 0 or target_padding_size[1] > 0:
807 
808  ysize_read_pixels = len(pixels)
809  nodata_value = fetch_band_nodata(band)
810 
811  # Apply columns padding
812  pad_cols = numpy.array([nodata_value] * target_padding_size[0])
813  for row in range (0, ysize_read_pixels):
814  out_line = numpy.append(pixels[row], pad_cols)
815  out_pixels[row] = out_line
816 
817  # Fill rows padding with nodata value
818  for row in range(ysize_read_pixels, ysize_read_pixels + target_padding_size[1]):
819  out_pixels[row].fill(nodata_value)
820  else:
821  out_pixels = pixels
822 
823  # XXX: Use for debugging only
824  #dump_block_numpy(out_pixels)
825 
826  hexwkb = binascii.hexlify(out_pixels)
827 
828  check_hex(hexwkb)
829  return hexwkb
830 
831 def wkblify_raster_level(options, ds, level, band_range, infile, i):
832  assert ds is not None
833  assert level >= 1
834  assert len(band_range) == 2
835 
836  band_from = band_range[0]
837  band_to = band_range[1]
838 
839  # Collect raster and block dimensions
840  raster_size = ( ds.RasterXSize, ds.RasterYSize )
841  if options.block_size is not None:
842  block_size = parse_block_size(options)
843  read_block_size = ( block_size[0] * level, block_size[1] * level)
844  grid_size = calculate_grid_size(raster_size, read_block_size)
845  else:
846  block_size = raster_size # Whole raster as a single block
847  read_block_size = block_size
848  grid_size = (1, 1)
849 
850  logit("MSG: Processing raster=%s using read_block_size=%s block_size=%s of grid=%s in level=%d\n" % \
851  (str(raster_size), str(read_block_size), str(block_size), str(grid_size), level))
852 
853  # Register base raster in RASTER_COLUMNS - SELECT AddRasterColumn();
854  if level == 1:
855  if i == 0 and options.create_table:
856  gt = get_gdal_geotransform(ds)
857  pixel_size = ( gt[1], gt[5] )
858  pixel_types = collect_pixel_types(ds, band_from, band_to)
859  nodata_values = collect_nodata_values(ds, band_from, band_to)
860  extent = calculate_bounding_box(ds, gt)
861  sql = make_sql_addrastercolumn(options, pixel_types, nodata_values,
862  pixel_size, block_size, extent)
863  options.output.write(sql)
864  gen_table = options.table
865 
866  else:
867  # Create overview table and register in RASTER_OVERVIEWS
868 
869  # CREATE TABLE o_<LEVEL>_<NAME> ( rid serial, options.column RASTER )
870  schema_table_names = make_sql_schema_table_names(options.table)
871  level_table_name = 'o_' + str(level) + '_' + schema_table_names[1]
872  level_table = schema_table_names[0] + '.' + level_table_name
873  if i == 0:
874  sql = make_sql_create_table(options, level_table, True)
875  options.output.write(sql)
876  sql = make_sql_register_overview(options, level_table_name, level)
877  options.output.write(sql)
878  gen_table = level_table
879 
880  # Write (original) raster to hex binary output
881  tile_count = 0
882  hexwkb = ''
883 
884  for ycell in range(0, grid_size[1]):
885  for xcell in range(0, grid_size[0]):
886 
887  xoff = xcell * read_block_size[0]
888  yoff = ycell * read_block_size[1]
889 
890  logit("MSG: --------- CELL #%04d\tindex = %d x %d\tdim = (%d x %d)\t(%d x %d) \t---------\n" % \
891  (tile_count, xcell, ycell, xoff, yoff, xoff + read_block_size[0], yoff + read_block_size[1]))
892 
893  if options.block_size is not None:
894  hexwkb = '' # Reset buffer as single INSERT per tile is generated
895  hexwkb += wkblify_raster_header(options, ds, level, (xoff, yoff),
896  block_size[0], block_size[1])
897  else:
898  hexwkb += wkblify_raster_header(options, ds, level, (xoff, yoff))
899 
900  for b in range(band_from, band_to):
901  band = ds.GetRasterBand(b)
902  assert band is not None, "Missing GDAL raster band %d" % b
903  logit("MSG: Band %d\n" % b)
904 
905  hexwkb += wkblify_band_header(options, band)
906  hexwkb += wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size, infile, b)
907 
908  # INSERT INTO
909  check_hex(hexwkb) # TODO: Remove to not to decrease performance
910  sql = make_sql_insert_raster(gen_table, options.column, hexwkb, options.filename, infile)
911  options.output.write(sql)
912 
913  tile_count = tile_count + 1
914 
915  return (gen_table, tile_count)
916 
917 def wkblify_raster(options, infile, i, previous_gt = None):
918  """Writes given raster dataset using GDAL features into HEX-encoded of
919  WKB for WKT Raster output."""
920 
921  assert infile is not None, "Input file is none, expected file name"
922  assert options.version == g_rt_version, "Error: invalid WKT Raster protocol version"
923  assert options.endian == NDR, "Error: invalid endianness, use little-endian (NDR) only"
924  assert options.srid >= -1, "Error: do you really want to specify SRID = %d" % options.srid
925 
926  # Open source raster file
927  ds = gdal.Open(infile, gdalc.GA_ReadOnly);
928  if ds is None:
929  sys.exit('Error: Cannot open input file: ' + str(infile))
930 
931  # By default, translate all raster bands
932 
933  # Calculate range for single-band request
934  if options.band is not None and options.band > 0:
935  band_range = ( options.band, options.band + 1 )
936  else:
937  band_range = ( 1, ds.RasterCount + 1 )
938 
939  # Compare this px size with previous one
940  current_gt = get_gdal_geotransform(ds)
941  if previous_gt is not None:
942  if previous_gt[1] != current_gt[1] or previous_gt[5] != current_gt[5]:
943  sys.exit('Error: Cannot load raster with different pixel size in the same raster table')
944 
945  # Generate requested overview level (base raster if level = 1)
946  summary = wkblify_raster_level(options, ds, options.overview_level, band_range, infile, i)
947  SUMMARY.append( summary )
948 
949  # Cleanup
950  ds = None
951 
952  return current_gt
953 
954 ################################################################################
955 
956 def main():
957 
958  (opts, args) = parse_command_line()
959 
960  global VERBOSE
961  VERBOSE = opts.verbose
962 
963  global SUMMARY
964  SUMMARY = []
965 
966  saved_out = sys.stdout
967  if type(opts.output) is str:
968  filename = opts.output
969  opts.output = open(filename, "w")
970 
971  # BEGIN
972  opts.output.write('BEGIN;\n')
973 
974  # If overviews requested, CREATE TABLE raster_overviews
975  if opts.create_raster_overviews_table:
977  opts.output.write(sql)
978 
979  # Base raster schema
980  if opts.overview_level == 1:
981  # DROP TABLE
982  if opts.drop_table:
983  sql = make_sql_drop_raster_table(opts.table)
984  opts.output.write(sql)
985 
986  # CREATE TABLE
987  if opts.create_table and opts.overview_level == 1:
988  sql = make_sql_create_table(opts)
989  opts.output.write(sql)
990 
991  # INSERT
992  i = 0
993 
994  # Burn all specified input raster files into single WKTRaster table
995  gt = None
996  for infile in opts.raster:
997  filelist = glob.glob(infile)
998  assert len(filelist) > 0, "No input raster files found for '" + str(infile) + "'"
999 
1000  for filename in filelist:
1001  logit("MSG: Dataset #%d: %s\n" % (i + 1, filename))
1002 
1003  # Write raster data to WKB and send it to opts.output
1004  gt = wkblify_raster(opts, filename.replace( '\\', '/') , i, gt)
1005  i += 1
1006 
1007  # INDEX
1008  if opts.index and SUMMARY is not None:
1009  sql = make_sql_create_gist(SUMMARY[0][0], opts.column)
1010  opts.output.write(sql)
1011 
1012  # COMMIT
1013  opts.output.write('END;\n')
1014 
1015  # VACUUM
1016  if opts.vacuum and SUMMARY is not None:
1017  sql = make_sql_vacuum(SUMMARY[0][0])
1018  opts.output.write(sql)
1019 
1020  # Cleanup
1021  if opts.output != sys.stdout:
1022  sys.stdout = saved_out
1023 
1024  print "------------------------------------------------------------"
1025  print " Summary of GDAL to PostGIS Raster processing:"
1026  print "------------------------------------------------------------"
1027  if i == 1:
1028  m = '%d (%s)' % (i, infile)
1029  else:
1030  m = '%d' % i
1031  print "Number of processed raster files: " + m
1032  print "List of generated tables (number of tiles):"
1033  i = 0
1034  for s in SUMMARY:
1035  i += 1
1036  print "%d\t%s (%d)" % (i, s[0], s[1])
1037 
1038 ################################################################################
1039 
1040 if __name__ == "__main__":
1041  main()
def make_sql_full_table_name(schema_table)
def calculate_geoxy(gt, xy)
def parse_command_line()
Definition: raster2pgsql.py:80
def gdt2pt(gdt)
def pt2fmt(pt)
def calculate_geoxy_level(gt, xy, level)
def dump_block_numpy(pixels)
def make_sql_register_overview(options, ov_table, ov_factor)
def calculate_block_pad_size(band, xoff, yoff, block_size)
def make_sql_drop_raster_table(table)
def collect_nodata_values(ds, band_from, band_to)
def make_sql_table_name(schema_table)
def calculate_block_size(ds, band_from, band_to)
def wkblify(fmt, data)
def fmt2printfmt(fmt)
def collect_pixel_types(ds, band_from, band_to)
def quote_sql_value(value)
SQL OPERATIONS.
def make_sql_value_array(values)
def parse_block_size(options)
def pt2numpy(pt)
def wkblify_raster_header
def fetch_band_nodata
def make_sql_create_table
def get_gdal_geotransform(ds)
def make_sql_vacuum(table)
def make_sql_create_raster_overviews(options)
def wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size, infile, bandidx)
def calculate_bounding_box(ds, gt)
def logit(msg)
def wkblify_band_header(options, band)
def wkblify_raster_level(options, ds, level, band_range, infile, i)
def quote_sql_name(name)
def make_sql_insert_raster(table, rast, hexwkb, insert_filename, file)
def calculate_grid_size(raster_size, block_size)
def make_sql_schema_table_names(schema_table)
def make_sql_addrastercolumn(options, pixeltypes, nodata, pixelsize, blocksize, extent)
def make_sql_create_gist(table, column)
def make_sql_drop_table(table)
def is_nan(x)
Definition: raster2pgsql.py:74
def calculate_overviews
RASTER OPERATIONS.
def calculate_overview_factor(ds, overview)