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