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