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