root/spike/wktraster/scripts/gdal2wktraster.py

Revision 5713, 37.7 KB (checked in by jorgearevalo, 20 months ago)

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