source: spike/wktraster/scripts/gdal2wktraster.py @ 5522

Last change on this file since 5522 was 5522, checked in by pracine, 6 years ago

Fix for ticket 348. ST_Envelope and ST_Box2d are now SQL functions of ST_ConvexHull. Removed useless rt_raster_get_envelope function from core.

FIX for ticket 478. Negate Y GDAL pixelsize when the image is not georeferenced.

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