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