PostGIS 3.6.2dev-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, overridden 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_Int8 : { 'name': 'PT_8BSI', 'id': 3 },
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
232def 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
245def 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
259def fmt2printfmt(fmt):
260 """Returns printf-like formatter for given binary data type specifier."""
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
273def 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
287def 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
296def 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
305def make_sql_value_array(values):
306 sql = "ARRAY["
307 for v in values:
308 if isinstance(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
316def 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
324def 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
329def 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
336def make_sql_drop_table(table):
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
342def make_sql_drop_raster_table(table):
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
354def 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
374def 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
385def 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
424def 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
437def make_sql_create_raster_overviews(options):
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
457def 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
473def make_sql_vacuum(table):
474 sql = 'VACUUM ANALYZE ' + make_sql_full_table_name(table) + ';\n'
475 return sql
476
477################################################################################
478# RASTER OPERATIONS
479
480def 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
500def 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
518def 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
531def 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
543def 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
563def 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
572def 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
587def get_gdal_geotransform(ds):
588 assert ds is not None
589 gt = list(ds.GetGeoTransform())
590 return tuple(gt)
591
592def 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
603def 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
610def calculate_bounding_box(ds, gt):
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
622def 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
631def 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
643def 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
654def 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
670def 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, "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 ### Endianness
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
722def 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
751def 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
831def 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
917def 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
956def 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 isinstance(opts.output, 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:
976 sql = make_sql_create_raster_overviews(opts)
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
1040if __name__ == "__main__":
1041 main()
str
#define str(s)
Definition lwgeom_functions_basic.c:44