[postgis-commits] svn - r3855 - spike/wktraster/scripts
postgis-commits at postgis.refractions.net
postgis-commits at postgis.refractions.net
Wed Mar 11 10:35:22 PDT 2009
Author: mloskot
Date: 2009-03-11 10:35:21 -0700 (Wed, 11 Mar 2009)
New Revision: 3855
Added:
spike/wktraster/scripts/gdal2wktraster.py
Removed:
spike/wktraster/scripts/gdal2wktraster.py
Log:
Setting executable attribute on gdal2wktraster.py
Deleted: spike/wktraster/scripts/gdal2wktraster.py
===================================================================
--- spike/wktraster/scripts/gdal2wktraster.py 2009-03-11 05:31:53 UTC (rev 3854)
+++ spike/wktraster/scripts/gdal2wktraster.py 2009-03-11 17:35:21 UTC (rev 3855)
@@ -1,349 +0,0 @@
-#! /usr/bin/env python
-#
-# $Id$
-#
-# This is a simple utility used to dump GDAL dataset into HEX WKB stream.
-# It's considered as a prototype of raster2pgsql tool planned to develop in future.
-# For more details about raster2pgsql tool, see Specification page:
-# http://www.cef-cfr.ca/index.php?n=Membres.PierreRacineWKTRasterSpecifications
-#
-# For this script to work you'll have to install the GDAL Python library
-# from http://pypi.python.org/pypi/GDAL/
-#
-# Copyright (C) 2009 Mateusz Loskot <mateusz at loskot.net>
-#
-# This program is free software; you can redistribute it and/or modify
-# it under the terms of the GNU General Public License as published by
-# the Free Software Foundation; either version 3 of the License, or
-# (at your option) any later version.
-#
-# This program is distributed in the hope that it will be useful,
-# but WITHOUT ANY WARRANTY; without even the implied warranty of
-# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
-# GNU General Public License for more details.
-#
-# You should have received a copy of the GNU General Public License
-# along with this program; if not, write to the Free Software
-# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
-#
-from osgeo import gdal
-from osgeo import osr
-import osgeo.gdalconst as gdalc
-from optparse import OptionParser
-import binascii
-import sys, glob
-
-###############################################################################
-# CONSTANTS - DO NOT CHANGE
-
-g_rt_endiannes = 1 # ALWAYS output little-endian
-g_rt_version = 0 # Version of WKTRaster protocol
-g_rt_column = 'rast' # Default name of column, overriden with -n option
-
-###############################################################################
-# UTILITIES
-
-def parse_command_line():
- # Parse command line options
- optprs = OptionParser(version="%prog $Revision: 3752 $")
- # Mandatory parameters
- optprs.add_option("-r", "--raster", dest="raster", action="append", default=None,
- help="append raster to list of input files, at least one raster file required. You may use wildcards (?,*) for specifying multiple files.")
- optprs.add_option("-t", "--table", dest="table", action="store", default=None,
- help="raster destination in form of [<schema>.]<table>, required")
- # Optional parameters
- optprs.add_option("-s", "--srid", dest="srid", action="store", type="int", default=-1,
- help="assign output raster with specified SRID")
- optprs.add_option("-b", "--band", dest="band", action="store", type="int", default=None,
- help="specify number of band to be extracted from raster file")
- optprs.add_option("-c", "--create", dest="create_table", action="store_true", default=False,
- help="create new table and populate it with raster(s), this is the default mode")
- optprs.add_option("-d", "--drop", dest="drop_table", action="store_true", default=False,
- help="drop table, create new one and populate with raster(s)")
- optprs.add_option("-f", "--field", dest="column", action="store", default=g_rt_column,
- help="specify name of raster data column, default is 'rast'")
- optprs.add_option("-I", "--index", dest="index", action="store_true", default=False,
- help="Create a GiST index on the raster column")
- optprs.add_option("-o", "--output", dest="output", action="store", default=sys.stdout,
- help="specify file for conversion output, otherwise send to stdout")
- optprs.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
- help="be excessively verbose and useful for debugging")
-
- (opts, args) = optprs.parse_args()
-
- # Validate options
- if opts.create_table and opts.drop_table:
- optprs.error("options -c and -d are mutually exclusive")
- if not opts.create_table and not opts.drop_table:
- opts.create_table = True
-
- if opts.raster is None:
- optprs.error("use option -r to specify at least one input raster. Wildcards (?,*) are accepted.")
-
- # XXX: Now, if --band=Nth, then only Nth band of all specified rasters is dumped/imported
- # This behavior can be changed to support only single-raster input if --band option used.
- #if opts.band is not None and len(opts.raster) > 1:
- # optprs.error("option -b requires single input raster only, specify -r option only once")
-
- if opts.table is None:
- optprs.error("use option -t to specify raster destination table")
- if len(opts.table.split('.')) > 2:
- optprs.error("invalid format of table name specified with option -t, expected [<schema>.]table")
-
- if opts.output is None:
- optprs.error("failed to initialise output file, try to use option -o explicitly")
-
- return (opts, args)
-
-
-def logit(msg):
- if VERBOSE is True:
- sys.stderr.write(msg)
-
-
-def gdt2pt(gdt):
- logit('MSG: Input GDAL pixel type: %s (%d)\n' % (gdal.GetDataTypeName(gdt), gdt))
- pixtypes = {
- gdalc.GDT_Byte : { 'name': 'PT_8BUI', 'id': 4 },
- gdalc.GDT_Int16 : { 'name': 'PT_16BSI', 'id': 5 },
- gdalc.GDT_UInt16 : { 'name': 'PT_16BUI', 'id': 6 },
- gdalc.GDT_Int32 : { 'name': 'PT_32BSI', 'id': 7 },
- gdalc.GDT_UInt32 : { 'name': 'PT_32BUI', 'id': 8 },
- gdalc.GDT_Float32 : { 'name': 'PT_32BF', 'id': 10 },
- gdalc.GDT_Float64 : { 'name': 'PT_64BF', 'id': 11 }
- }
- logit('MSG: Output WKTRaster pixel type: %(name)s (%(id)d)\n' % (pixtypes.get(gdt, 13)))
- return pixtypes.get(gdt, 13)
-
-
-def pt2fmt(pt):
- fmttypes = {
- 4: 'B', # PT_8BUI
- 5: 'h', # PT_16BSI
- 6: 'H', # PT_16BUI
- 7: 'i', # PT_32BSI
- 8: 'I', # PT_32BUI
- 10: 'f', # PT_32BF
- 11: 'd' # PT_64BF
- }
- return fmttypes.get(pt, 'x')
-
-
-def fmt2printfmt(fmt):
- fmttypes = {
- 'B': '%d', # PT_8BUI
- 'h': '%d', # PT_16BSI
- 'H': '%d', # PT_16BUI
- 'i': '%d', # PT_32BSI
- 'I': '%d', # PT_32BUI
- 'f': '%f', # PT_32BF
- 'd': '%f' # PT_64BF
- }
- return fmttypes.get(fmt, 'f')
-
-
-def wkblify(fmt, data):
- import struct
-
- # Binary to HEX
- fmt_little = '<' +fmt
- hexstr = binascii.hexlify(struct.pack(fmt_little, data)).upper()
-
- # String'ify raw value for log
- valfmt = '\'' + fmt2printfmt(fmt) + '\''
- val = valfmt % data
- logit('HEX (\'fmt=%s\', bytes=%d, val=%s):\t\t%s\n' % (fmt, len(hexstr) / 2, str(val), hexstr))
-
- return hexstr
-
-
-def wkblify_raster_header(rt_version, rt_endian, rt_srid, gdalds):
- assert gdalds is not None, "Error: No GDAL dataset"
-
- # Collect GeoReference information
- gt = gdalds.GetGeoTransform()
- rt_scale = ( gt[1], gt[5] ) # X/Y pixel resolution
- rt_ip = ( gt[0], gt[3] ) # X/Y upper left pixel center
- rt_skew = ( gt[2], gt[4] ) # X-/Y-axis rotation
-
- # TODO: What about looking for SRID based on SRS in WKT?
- #srs = osr.SpatialReference()
- #srs.ImportFromWkt(ds.GetProjection())
-
- # Burn input raster as WKTRaster WKB format
- hexwkb = ''
- ### Endiannes
- hexwkb += wkblify('B', rt_endian)
- ### Version
- hexwkb += wkblify('H', rt_version)
- ### Number of bands
- hexwkb += wkblify('H', gdalds.RasterCount)
- assert 5 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
- ### Georeference
- hexwkb += wkblify('d', rt_scale[0])
- hexwkb += wkblify('d', rt_scale[1])
- hexwkb += wkblify('d', rt_ip[0])
- hexwkb += wkblify('d', rt_ip[1])
- hexwkb += wkblify('d', rt_skew[0])
- hexwkb += wkblify('d', rt_skew[1])
- hexwkb += wkblify('i', rt_srid)
- assert 57 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
- ### Number of columns and rows
- hexwkb += wkblify('H', gdalds.RasterXSize)
- hexwkb += wkblify('H', gdalds.RasterYSize)
- assert 61 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
-
- return hexwkb
-
-
-def wkblify_band(band):
- assert band is not None, "Error: No raster band"
-
- hexwkb = ""
- pixtype = gdt2pt(band.DataType)['id']
- hexwkb += wkblify('B', pixtype)
-
- # FIXME: What if input dataset does not specify NODATA value? Zero?
- nodata = 0
- if band.GetNoDataValue() is not None:
- nodata = band.GetNoDataValue()
- hexwkb += wkblify(pt2fmt(pixtype), nodata)
-
- # Rows of pixels of raster
- for row in range(0, band.YSize):
- scanline = band.ReadRaster(0, row, band.XSize, 1, band.XSize, 1)
- hexwkb += binascii.hexlify(scanline)
- logit('.')
- logit('\n')
-
- return hexwkb
-
-def wkblify_raster(infile, band_num, rt_version, rt_endian, rt_srid):
- assert rt_version == 0, "Error: invalid rt_version"
- assert rt_endian == 1, "Error: invalid rt_endian specifier, use little-endian only"
- assert rt_srid >= -1, "Error: do you really want to specify SRID = %d" % rt_srid
-
- # Open source raster file
- ds = gdal.Open(infile, gdalc.GA_ReadOnly);
- if ds is None:
- sys.exit('Error: Cannot open input file: ' + str(infile))
-
- # Burn input raster as WKTRaster WKB format
- hexwkb = ''
- hexwkb = wkblify_raster_header(rt_version, rt_endian, rt_srid, ds)
-
- ### Raster bands
- band_from = 1
- band_to = ds.RasterCount + 1
- if band_num is not None and band_num > 0:
- band_from = band_num
- band_to = band_num + 1
-
- for bn in range(band_from, band_to):
- band = ds.GetRasterBand(bn)
- if band is None:
- logit("Cannot fetch raster band of index %d" % bn)
- continue
-
- logit('MSG: Band #%d\n' % bn)
- hexwkb += wkblify_band(band)
- band = None
- ds = None
-
- return hexwkb
-
-
-def make_sql_full_table_name(schema_table):
- st = schema_table.split('.')
- if len(st) == 1:
- st.insert(0, "public")
- assert len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
- table = "\"%s\".\"%s\"" % (st[0], st[1])
- return table
-
-def make_sql_table_name(schema_table):
- st = schema_table.split('.')
- assert len(st) == 1 or len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
- if len(st) == 2:
- return st[1]
- return st[0]
-
-def make_sql_drop_table(table):
- sql = "DROP TABLE IF EXISTS %s CASCADE;\n" % make_sql_full_table_name(table)
- logit("SQL: %s" % sql)
- return sql
-
-
-def make_sql_create_table(table, rast):
- sql = "CREATE TABLE %s (rid serial PRIMARY KEY, \"%s\" RASTER );\n" % (make_sql_full_table_name(table), rast)
- logit("SQL: %s" % sql)
- return sql
-
-
-def make_sql_insert_raster(table, rast, hexwkb):
- sql = "INSERT INTO %s ( %s ) VALUES ( (\'%s\')::raster );\n" % (make_sql_full_table_name(table), rast, hexwkb)
- # NOTE: Usually, no need for such detailed verbosity
- #logit("SQL: %s" % sql)
- return sql
-
-###############################################################################
-
-def main():
-
- (opts, args) = parse_command_line()
-
- global VERBOSE
- VERBOSE = opts.verbose
-
- fout = None
- saved_out = None
- if type(opts.output) is str:
- fout = open(opts.output, "w")
- else:
- saved_out = sys.stdout
- fout = opts.output
-
- # DROP TABLE
- if opts.drop_table:
- sql = make_sql_drop_table(opts.table)
- fout.write(sql)
-
- # CREATE TABLE
- sql = make_sql_create_table(opts.table, opts.column)
- fout.write(sql)
-
- # Burn all specified input raster files into single WKTRaster table
- i = 0
-
- for infile in opts.raster:
- filelist = glob.glob(infile)
- for file in filelist:
- logit("MSG: Raster #%d: %s\n" % (i + 1, file))
-
- hexwkb = wkblify_raster(file, opts.band, g_rt_version, g_rt_endiannes, opts.srid)
- assert len(hexwkb) > 0, "Error: No HEX WKB generated"
- assert (len(hexwkb) % 2 == 0), "Error: Invalid size of generated HEX WKB"
-
- # INSERT INTO ... (rast) VALUES (...)
- sql = make_sql_insert_raster(opts.table, opts.column, hexwkb)
- fout.write(sql)
-
- i += 1
-
- if opts.index:
- gist_table = make_sql_table_name(opts.table)
- target_table = make_sql_full_table_name(opts.table)
- logit("MSG: Using GiST spatial index on %s\n" % gist_table)
- sql = "CREATE INDEX \"%s_%s_gist\" ON %s USING GIST (rt_raster_envelope(%s));\n" % \
- (gist_table, opts.column, target_table, opts.column)
- fout.write(sql)
-
- # Cleanup
- sys.stdout = saved_out
-
- logit("MSG: Number of processed raster files: %d\n" % i)
-
-###############################################################################
-
-if __name__ == "__main__":
- main()
Added: spike/wktraster/scripts/gdal2wktraster.py
===================================================================
--- spike/wktraster/scripts/gdal2wktraster.py 2009-03-11 05:31:53 UTC (rev 3854)
+++ spike/wktraster/scripts/gdal2wktraster.py 2009-03-11 17:35:21 UTC (rev 3855)
@@ -0,0 +1,349 @@
+#! /usr/bin/env python
+#
+# $Id: gdal2wktraster.py 3834 2009-03-10 18:01:23Z mloskot $
+#
+# This is a simple utility used to dump GDAL dataset into HEX WKB stream.
+# It's considered as a prototype of raster2pgsql tool planned to develop in future.
+# For more details about raster2pgsql tool, see Specification page:
+# http://www.cef-cfr.ca/index.php?n=Membres.PierreRacineWKTRasterSpecifications
+#
+# For this script to work you'll have to install the GDAL Python library
+# from http://pypi.python.org/pypi/GDAL/
+#
+# Copyright (C) 2009 Mateusz Loskot <mateusz at loskot.net>
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program; if not, write to the Free Software
+# Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+#
+from osgeo import gdal
+from osgeo import osr
+import osgeo.gdalconst as gdalc
+from optparse import OptionParser
+import binascii
+import sys, glob
+
+###############################################################################
+# CONSTANTS - DO NOT CHANGE
+
+g_rt_endiannes = 1 # ALWAYS output little-endian
+g_rt_version = 0 # Version of WKTRaster protocol
+g_rt_column = 'rast' # Default name of column, overriden with -n option
+
+###############################################################################
+# UTILITIES
+
+def parse_command_line():
+ # Parse command line options
+ optprs = OptionParser(version="%prog $Revision: 3752 $")
+ # Mandatory parameters
+ optprs.add_option("-r", "--raster", dest="raster", action="append", default=None,
+ help="append raster to list of input files, at least one raster file required. You may use wildcards (?,*) for specifying multiple files.")
+ optprs.add_option("-t", "--table", dest="table", action="store", default=None,
+ help="raster destination in form of [<schema>.]<table>, required")
+ # Optional parameters
+ optprs.add_option("-s", "--srid", dest="srid", action="store", type="int", default=-1,
+ help="assign output raster with specified SRID")
+ optprs.add_option("-b", "--band", dest="band", action="store", type="int", default=None,
+ help="specify number of band to be extracted from raster file")
+ optprs.add_option("-c", "--create", dest="create_table", action="store_true", default=False,
+ help="create new table and populate it with raster(s), this is the default mode")
+ optprs.add_option("-d", "--drop", dest="drop_table", action="store_true", default=False,
+ help="drop table, create new one and populate with raster(s)")
+ optprs.add_option("-f", "--field", dest="column", action="store", default=g_rt_column,
+ help="specify name of raster data column, default is 'rast'")
+ optprs.add_option("-I", "--index", dest="index", action="store_true", default=False,
+ help="Create a GiST index on the raster column")
+ optprs.add_option("-o", "--output", dest="output", action="store", default=sys.stdout,
+ help="specify file for conversion output, otherwise send to stdout")
+ optprs.add_option("-v", "--verbose", dest="verbose", action="store_true", default=False,
+ help="be excessively verbose and useful for debugging")
+
+ (opts, args) = optprs.parse_args()
+
+ # Validate options
+ if opts.create_table and opts.drop_table:
+ optprs.error("options -c and -d are mutually exclusive")
+ if not opts.create_table and not opts.drop_table:
+ opts.create_table = True
+
+ if opts.raster is None:
+ optprs.error("use option -r to specify at least one input raster. Wildcards (?,*) are accepted.")
+
+ # XXX: Now, if --band=Nth, then only Nth band of all specified rasters is dumped/imported
+ # This behavior can be changed to support only single-raster input if --band option used.
+ #if opts.band is not None and len(opts.raster) > 1:
+ # optprs.error("option -b requires single input raster only, specify -r option only once")
+
+ if opts.table is None:
+ optprs.error("use option -t to specify raster destination table")
+ if len(opts.table.split('.')) > 2:
+ optprs.error("invalid format of table name specified with option -t, expected [<schema>.]table")
+
+ if opts.output is None:
+ optprs.error("failed to initialise output file, try to use option -o explicitly")
+
+ return (opts, args)
+
+
+def logit(msg):
+ if VERBOSE is True:
+ sys.stderr.write(msg)
+
+
+def gdt2pt(gdt):
+ logit('MSG: Input GDAL pixel type: %s (%d)\n' % (gdal.GetDataTypeName(gdt), gdt))
+ pixtypes = {
+ gdalc.GDT_Byte : { 'name': 'PT_8BUI', 'id': 4 },
+ gdalc.GDT_Int16 : { 'name': 'PT_16BSI', 'id': 5 },
+ gdalc.GDT_UInt16 : { 'name': 'PT_16BUI', 'id': 6 },
+ gdalc.GDT_Int32 : { 'name': 'PT_32BSI', 'id': 7 },
+ gdalc.GDT_UInt32 : { 'name': 'PT_32BUI', 'id': 8 },
+ gdalc.GDT_Float32 : { 'name': 'PT_32BF', 'id': 10 },
+ gdalc.GDT_Float64 : { 'name': 'PT_64BF', 'id': 11 }
+ }
+ logit('MSG: Output WKTRaster pixel type: %(name)s (%(id)d)\n' % (pixtypes.get(gdt, 13)))
+ return pixtypes.get(gdt, 13)
+
+
+def pt2fmt(pt):
+ fmttypes = {
+ 4: 'B', # PT_8BUI
+ 5: 'h', # PT_16BSI
+ 6: 'H', # PT_16BUI
+ 7: 'i', # PT_32BSI
+ 8: 'I', # PT_32BUI
+ 10: 'f', # PT_32BF
+ 11: 'd' # PT_64BF
+ }
+ return fmttypes.get(pt, 'x')
+
+
+def fmt2printfmt(fmt):
+ fmttypes = {
+ 'B': '%d', # PT_8BUI
+ 'h': '%d', # PT_16BSI
+ 'H': '%d', # PT_16BUI
+ 'i': '%d', # PT_32BSI
+ 'I': '%d', # PT_32BUI
+ 'f': '%f', # PT_32BF
+ 'd': '%f' # PT_64BF
+ }
+ return fmttypes.get(fmt, 'f')
+
+
+def wkblify(fmt, data):
+ import struct
+
+ # Binary to HEX
+ fmt_little = '<' +fmt
+ hexstr = binascii.hexlify(struct.pack(fmt_little, data)).upper()
+
+ # String'ify raw value for log
+ valfmt = '\'' + fmt2printfmt(fmt) + '\''
+ val = valfmt % data
+ logit('HEX (\'fmt=%s\', bytes=%d, val=%s):\t\t%s\n' % (fmt, len(hexstr) / 2, str(val), hexstr))
+
+ return hexstr
+
+
+def wkblify_raster_header(rt_version, rt_endian, rt_srid, gdalds):
+ assert gdalds is not None, "Error: No GDAL dataset"
+
+ # Collect GeoReference information
+ gt = gdalds.GetGeoTransform()
+ rt_scale = ( gt[1], gt[5] ) # X/Y pixel resolution
+ rt_ip = ( gt[0], gt[3] ) # X/Y upper left pixel center
+ rt_skew = ( gt[2], gt[4] ) # X-/Y-axis rotation
+
+ # TODO: What about looking for SRID based on SRS in WKT?
+ #srs = osr.SpatialReference()
+ #srs.ImportFromWkt(ds.GetProjection())
+
+ # Burn input raster as WKTRaster WKB format
+ hexwkb = ''
+ ### Endiannes
+ hexwkb += wkblify('B', rt_endian)
+ ### Version
+ hexwkb += wkblify('H', rt_version)
+ ### Number of bands
+ hexwkb += wkblify('H', gdalds.RasterCount)
+ assert 5 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
+ ### Georeference
+ hexwkb += wkblify('d', rt_scale[0])
+ hexwkb += wkblify('d', rt_scale[1])
+ hexwkb += wkblify('d', rt_ip[0])
+ hexwkb += wkblify('d', rt_ip[1])
+ hexwkb += wkblify('d', rt_skew[0])
+ hexwkb += wkblify('d', rt_skew[1])
+ hexwkb += wkblify('i', rt_srid)
+ assert 57 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
+ ### Number of columns and rows
+ hexwkb += wkblify('H', gdalds.RasterXSize)
+ hexwkb += wkblify('H', gdalds.RasterYSize)
+ assert 61 == (len(hexwkb) / 2), "Error: Invalid length of HEX WKB"
+
+ return hexwkb
+
+
+def wkblify_band(band):
+ assert band is not None, "Error: No raster band"
+
+ hexwkb = ""
+ pixtype = gdt2pt(band.DataType)['id']
+ hexwkb += wkblify('B', pixtype)
+
+ # FIXME: What if input dataset does not specify NODATA value? Zero?
+ nodata = 0
+ if band.GetNoDataValue() is not None:
+ nodata = band.GetNoDataValue()
+ hexwkb += wkblify(pt2fmt(pixtype), nodata)
+
+ # Rows of pixels of raster
+ for row in range(0, band.YSize):
+ scanline = band.ReadRaster(0, row, band.XSize, 1, band.XSize, 1)
+ hexwkb += binascii.hexlify(scanline)
+ logit('.')
+ logit('\n')
+
+ return hexwkb
+
+def wkblify_raster(infile, band_num, rt_version, rt_endian, rt_srid):
+ assert rt_version == 0, "Error: invalid rt_version"
+ assert rt_endian == 1, "Error: invalid rt_endian specifier, use little-endian only"
+ assert rt_srid >= -1, "Error: do you really want to specify SRID = %d" % rt_srid
+
+ # Open source raster file
+ ds = gdal.Open(infile, gdalc.GA_ReadOnly);
+ if ds is None:
+ sys.exit('Error: Cannot open input file: ' + str(infile))
+
+ # Burn input raster as WKTRaster WKB format
+ hexwkb = ''
+ hexwkb = wkblify_raster_header(rt_version, rt_endian, rt_srid, ds)
+
+ ### Raster bands
+ band_from = 1
+ band_to = ds.RasterCount + 1
+ if band_num is not None and band_num > 0:
+ band_from = band_num
+ band_to = band_num + 1
+
+ for bn in range(band_from, band_to):
+ band = ds.GetRasterBand(bn)
+ if band is None:
+ logit("Cannot fetch raster band of index %d" % bn)
+ continue
+
+ logit('MSG: Band #%d\n' % bn)
+ hexwkb += wkblify_band(band)
+ band = None
+ ds = None
+
+ return hexwkb
+
+
+def make_sql_full_table_name(schema_table):
+ st = schema_table.split('.')
+ if len(st) == 1:
+ st.insert(0, "public")
+ assert len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
+ table = "\"%s\".\"%s\"" % (st[0], st[1])
+ return table
+
+def make_sql_table_name(schema_table):
+ st = schema_table.split('.')
+ assert len(st) == 1 or len(st) == 2, "Invalid format of table name, expected [<schema>.]table"
+ if len(st) == 2:
+ return st[1]
+ return st[0]
+
+def make_sql_drop_table(table):
+ sql = "DROP TABLE IF EXISTS %s CASCADE;\n" % make_sql_full_table_name(table)
+ logit("SQL: %s" % sql)
+ return sql
+
+
+def make_sql_create_table(table, rast):
+ sql = "CREATE TABLE %s (rid serial PRIMARY KEY, \"%s\" RASTER );\n" % (make_sql_full_table_name(table), rast)
+ logit("SQL: %s" % sql)
+ return sql
+
+
+def make_sql_insert_raster(table, rast, hexwkb):
+ sql = "INSERT INTO %s ( %s ) VALUES ( (\'%s\')::raster );\n" % (make_sql_full_table_name(table), rast, hexwkb)
+ # NOTE: Usually, no need for such detailed verbosity
+ #logit("SQL: %s" % sql)
+ return sql
+
+###############################################################################
+
+def main():
+
+ (opts, args) = parse_command_line()
+
+ global VERBOSE
+ VERBOSE = opts.verbose
+
+ fout = None
+ saved_out = None
+ if type(opts.output) is str:
+ fout = open(opts.output, "w")
+ else:
+ saved_out = sys.stdout
+ fout = opts.output
+
+ # DROP TABLE
+ if opts.drop_table:
+ sql = make_sql_drop_table(opts.table)
+ fout.write(sql)
+
+ # CREATE TABLE
+ sql = make_sql_create_table(opts.table, opts.column)
+ fout.write(sql)
+
+ # Burn all specified input raster files into single WKTRaster table
+ i = 0
+
+ for infile in opts.raster:
+ filelist = glob.glob(infile)
+ for file in filelist:
+ logit("MSG: Raster #%d: %s\n" % (i + 1, file))
+
+ hexwkb = wkblify_raster(file, opts.band, g_rt_version, g_rt_endiannes, opts.srid)
+ assert len(hexwkb) > 0, "Error: No HEX WKB generated"
+ assert (len(hexwkb) % 2 == 0), "Error: Invalid size of generated HEX WKB"
+
+ # INSERT INTO ... (rast) VALUES (...)
+ sql = make_sql_insert_raster(opts.table, opts.column, hexwkb)
+ fout.write(sql)
+
+ i += 1
+
+ if opts.index:
+ gist_table = make_sql_table_name(opts.table)
+ target_table = make_sql_full_table_name(opts.table)
+ logit("MSG: Using GiST spatial index on %s\n" % gist_table)
+ sql = "CREATE INDEX \"%s_%s_gist\" ON %s USING GIST (rt_raster_envelope(%s));\n" % \
+ (gist_table, opts.column, target_table, opts.column)
+ fout.write(sql)
+
+ # Cleanup
+ sys.stdout = saved_out
+
+ logit("MSG: Number of processed raster files: %d\n" % i)
+
+###############################################################################
+
+if __name__ == "__main__":
+ main()
Property changes on: spike/wktraster/scripts/gdal2wktraster.py
___________________________________________________________________
Name: svn:executable
+ *
More information about the postgis-commits
mailing list