Ticket #227: gdal2wktraster_outdb_support.patch

*** gdal2wktraster.py   2009-08-16 19:49:33.000000000 +0200
--- gdal2wktraster_outdb.py     2009-08-16 19:59:02.000000000 +0200
***************
*** 2,7 ****
--- 2,11 ----
  #
  # $Id: gdal2wktraster.py 4326 2009-07-22 15:14:41Z mloskot $
  #
+ # 2009-07-26: Support for out-db rasters (experimental). Not under svn control.
+ # This file should be revised by community-
+ # Author: Jorge Arevalo jorgearevalo@gis4free.org
+ # 
  # 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.
***************
*** 63,68 ****
--- 67,80 ----
  g_rt_catalog = ''
  g_rt_schema = 'public'
  
+ # Suffix added to the table name to generate a name for out-db raster
+ g_rt_outdb_sufix = '_outdb.tif'
+ 
+ # Default path for outdb raster
+ g_rt_outdb_path = '/tmp'
+ 
+ # Path separator, changes in each os
+ path_separator = '/'
  ################################################################################
  # UTILITIES
  VERBOSE = False
***************
*** 93,99 ****
      grp_r.add_option("-k", "--block-size", dest="block_size", action="store", default=None,
                       help="uses regular blocking mode with size of block used to tile input raster, specified as WIDTHxHEIGHT")
      grp_r.add_option("-R", "--register", dest="register", action="store_true", default=False, 
!                      help="register the raster as a filesystem (out-db) raster")
      grp_r.add_option("-l", "--overview-level", dest="overview_level", action="store", type="int", default=1,
                       help='Create overview tables named as ov_<RASTER TABLE>_<LEVEL> and '
                       'populate with GDAL-provided overviews (regular blocking only)')
--- 105,113 ----
      grp_r.add_option("-k", "--block-size", dest="block_size", action="store", default=None,
                       help="uses regular blocking mode with size of block used to tile input raster, specified as WIDTHxHEIGHT")
      grp_r.add_option("-R", "--register", dest="register", action="store_true", default=False, 
!                      help="register the raster as a filesystem (out-db) raster")                     
!     grp_r.add_option("-p", "--path", dest="outdb_path", action="store", default=None, 
!                      help="path to register the raster as a filesystem (out-db) raster")                                          
      grp_r.add_option("-l", "--overview-level", dest="overview_level", action="store", type="int", default=1,
                       help='Create overview tables named as ov_<RASTER TABLE>_<LEVEL> and '
                       'populate with GDAL-provided overviews (regular blocking only)')
***************
*** 158,163 ****
--- 172,181 ----
  
      if opts.create_raster_overviews_table and opts.overview_level <= 1:
          prs.error('create table for RASTER_OVERVIEWS available only if overviews import requested')
+         
+     if opts.outdb_path is not None and opts.register == False:
+         prs.error('you need to specify -R option if you want to register the raster as outdb raster')
+         
  
      # 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.
***************
*** 676,682 ****
      rt_skew = ( gt[2], gt[4] )
      rt_scale = ( gt[1] * level, gt[5] * level )
      
-     
      # TODO: Any way to lookup for SRID based on SRS in WKT?
      #srs = osr.SpatialReference()
      #srs.ImportFromWkt(ds.GetProjection())
--- 694,699 ----
***************
*** 738,823 ****
      check_hex(hexwkb)
      return hexwkb
  
! def wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size):
!     """Writes band of given GDAL dataset into HEX-encoded WKB for WKT Raster output."""
      assert band is not None, "Error: Missing GDAL raster band"
  
      hexwkb = ''
      
!     if options.register:
!         # Off-db raster
!         # TODO: Do we want to handle options.overview_level? --mloskot
!         # TODO: Where bandidx and ds come from? --mloskot
!         hexwkb += wkblify('B', bandidx - 1)
!         filepath = os.path.abspath((ds.GetFileList()[0]).replace('\\', '\\\\'))
!         hexwkb += wkblify(str(len(filepath)) + 's', filepath)
!         hexwkb += wkblify('B', 0)
      else:
!         # In-db raster
  
!         # Real NODATA value or Zero'ed byte(s)
!         hexnodata = wkblify_band_nodata(band)
!         bytes_per_pixel = len(hexnodata) / 2 # bytes per pixel used to validate hex string
! 
!         # Right most column and bottom most row of blocks have
!         # portions that extend beyond the raster
!         read_padding_size = calculate_block_pad_size(band, xoff, yoff, read_block_size)
!         valid_read_block_size = ( read_block_size[0] - read_padding_size[0],
!                                   read_block_size[1] - read_padding_size[1] )
! 
! 
!         if read_padding_size[0] > 0 or read_padding_size[1] > 0:
!             target_block_size = (valid_read_block_size[0] / level, valid_read_block_size[1] / level)
!             target_padding_size = (read_padding_size[0] / level, read_padding_size[1] / level)
!         else:
!             target_block_size = block_size
!             target_padding_size = ( 0, 0 )
  
!         logit('MSG: Normalize read_block=%s for level=%d to valid_read_block=%s with padding=%s\n' % \
!               (read_block_size, level, valid_read_block_size, read_padding_size))
!         logit('MSG: Normalize target_block=%s for level=%d to valid_target_block=%s with padding=%s\n' % \
!               (block_size, level, target_block_size, target_padding_size))
!         logit('MSG: ReadAsArray( %d, %d, %s, %s)\n' % \
!               (xoff, yoff, str(valid_read_block_size), str(target_block_size)))
  
!         assert valid_read_block_size[0] > 0 and valid_read_block_size[1] > 0
!         assert target_block_size[0] > 0 and target_block_size[1] > 0
  
!         pixels = band.ReadAsArray(xoff, yoff, valid_read_block_size[0], valid_read_block_size[1],
!                                   target_block_size[0], target_block_size[1])
  
!         # XXX: Use for debugging only
!         #dump_block_numpy(pixels)
  
!         out_pixels = numpy.zeros((block_size[1], block_size[0]), pt2numpy(band.DataType))
  
!         logit('MSG: Read valid source:\t%d x %d\n' % (len(pixels[0]), len(pixels)))
!         logit('MSG: Write into block:\t%d x %d\n' % (len(out_pixels[0]), len(out_pixels)))
          
!         if target_padding_size[0] > 0 or target_padding_size[1] > 0:
  
!             ysize_read_pixels = len(pixels)
!             nodata_value = fetch_band_nodata(band)
  
!             # Apply columns padding
!             pad_cols = numpy.array([nodata_value] * target_padding_size[0])
!             for row in range (0, ysize_read_pixels):
!                 out_line = numpy.append(pixels[row], pad_cols)
!                 out_pixels[row] = out_line
! 
!             # Fill rows padding with nodata value
!             for row in range(ysize_read_pixels, ysize_read_pixels + target_padding_size[1]):
!                 out_pixels[row].fill(nodata_value)
!         else:
!             out_pixels = pixels
  
!         # XXX: Use for debugging only
!         #dump_block_numpy(out_pixels)
  
!         hexwkb = binascii.hexlify(out_pixels)
  
      check_hex(hexwkb)
      return hexwkb
  
  def wkblify_raster_level(options, ds, level, band_range):
      assert ds is not None
--- 755,920 ----
      check_hex(hexwkb)
      return hexwkb
  
! 
! def wkblify_indb_band(band, level, xoff, yoff, read_block_size, block_size):
!     """Writes indb band of given GDAL dataset into HEX-encoded WKB for WKT Raster output."""
      assert band is not None, "Error: Missing GDAL raster band"
  
      hexwkb = ''
      
!     # Real NODATA value or Zero'ed byte(s)
!     hexnodata = wkblify_band_nodata(band)
!     bytes_per_pixel = len(hexnodata) / 2 # bytes per pixel used to validate hex string
! 
!     # Right most column and bottom most row of blocks have
!     # portions that extend beyond the raster
!     read_padding_size = calculate_block_pad_size(band, xoff, yoff, read_block_size)
!     valid_read_block_size = ( read_block_size[0] - read_padding_size[0],
!                                 read_block_size[1] - read_padding_size[1] )
! 
! 
!     if read_padding_size[0] > 0 or read_padding_size[1] > 0:
!         target_block_size = (valid_read_block_size[0] / level, valid_read_block_size[1] / level)
!         target_padding_size = (read_padding_size[0] / level, read_padding_size[1] / level)
      else:
!         target_block_size = block_size
!         target_padding_size = ( 0, 0 )
  
!     logit('MSG: Normalize read_block=%s for level=%d to valid_read_block=%s with padding=%s\n' % \
!             (read_block_size, level, valid_read_block_size, read_padding_size))
!     logit('MSG: Normalize target_block=%s for level=%d to valid_target_block=%s with padding=%s\n' % \
!             (block_size, level, target_block_size, target_padding_size))
!     logit('MSG: ReadAsArray( %d, %d, %s, %s)\n' % \
!             (xoff, yoff, str(valid_read_block_size), str(target_block_size)))
! 
!     assert valid_read_block_size[0] > 0 and valid_read_block_size[1] > 0
!     assert target_block_size[0] > 0 and target_block_size[1] > 0
  
!     pixels = band.ReadAsArray(xoff, yoff, valid_read_block_size[0], valid_read_block_size[1],
!                                 target_block_size[0], target_block_size[1])
  
!     # XXX: Use for debugging only
!     #dump_block_numpy(pixels)
  
!     out_pixels = numpy.zeros((block_size[1], block_size[0]), pt2numpy(band.DataType))
  
!     logit('MSG: Read valid source:\t%d x %d\n' % (len(pixels[0]), len(pixels)))
!     logit('MSG: Write into block:\t%d x %d\n' % (len(out_pixels[0]), len(out_pixels)))
!     
!     if target_padding_size[0] > 0 or target_padding_size[1] > 0:
! 
!         ysize_read_pixels = len(pixels)
!         nodata_value = fetch_band_nodata(band)
  
!         # Apply columns padding
!         pad_cols = numpy.array([nodata_value] * target_padding_size[0])
!         for row in range (0, ysize_read_pixels):
!             out_line = numpy.append(pixels[row], pad_cols)
!             out_pixels[row] = out_line
! 
!         # Fill rows padding with nodata value
!         for row in range(ysize_read_pixels, ysize_read_pixels + target_padding_size[1]):
!             out_pixels[row].fill(nodata_value)
!     else:
!         out_pixels = pixels
  
!     # XXX: Use for debugging only
!     #dump_block_numpy(out_pixels)
! 
!     hexwkb = binascii.hexlify(out_pixels)
! 
!     return hexwkb
!     
! def wkblify_outdb_band(options, level, bandidx, band, out_ds, write_to_file):
!     """Writes outdb band of given GDAL dataset into HEX-encoded WKB for WKT Raster output."""
!     assert band is not None, "Error: Missing GDAL raster band"
!     
!     # If path to outdb file wasn't specified, g_rt_outdb_path is used
!     if options.outdb_path is None:        
!         options.outdb_path = g_rt_outdb_path
          
!     filepath = options.outdb_path + path_separator + options.table + g_rt_outdb_sufix
  
!     # Write the band to the file
!     if write_to_file:
!         # Get band pixels
!         pixels = band.ReadAsArray(0, 0, band.XSize, band.YSize,
!                                   band.XSize, band.YSize)
!                          
!         # Write pixels in file                
!         out_ds.GetRasterBand(bandidx).WriteArray(pixels)
  
!     hexwkb = ''
!     
!     # 0-based band number
!     hexwkb += wkblify('B', bandidx - 1)
!     
!     # Path to the file, as 0-end string
!     hexwkb += wkblify(str(len(filepath)) + 's', filepath)
!     hexwkb += wkblify('B', 0)
  
!     return hexwkb
  
! 
! def wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size, out_ds, bandidx):
!     """Writes band of given GDAL dataset into HEX-encoded WKB for WKT Raster output."""
!     assert band is not None, "Error: Missing GDAL raster band"
! 
!     hexwkb = ''
!     
!     # Off-db raster
!     if options.register:
!         if xoff == 0 and yoff == 0:
!             # Only write the band in file the first time called with each band
!             hexwkb += wkblify_outdb_band(options, level, bandidx, band, out_ds, True)
!         else:
!             hexwkb += wkblify_outdb_band(options, level, bandidx, band, out_ds, False)
!     
!     # In-db raster
!     else:       
!         hexwkb += wkblify_indb_band(band, level, xoff, yoff, read_block_size, block_size)
  
      check_hex(hexwkb)
      return hexwkb
+     
+ def create_empty_outdb_dataset(options, ds, band_range):
+     """Creates a blank tiff raster based on the original one."""
+     assert ds is not None
+     assert len(band_range) == 2
+ 
+     band_from = band_range[0]
+     band_to = band_range[1]
+     
+     # Create new dataset to write TIFF file
+     outdb_file_name = options.table + g_rt_outdb_sufix
+     tiff_width = ds.RasterXSize
+     tiff_height = ds.RasterYSize
+                     
+     tiff_nbands = band_to - band_from
+     raster = gdal.GetDriverByName('GTiff')
+     out_ds = raster.Create(outdb_file_name, tiff_width, tiff_height, tiff_nbands)
+     assert out_ds is not None, "Couldn't create out-db Dataset"
+     
+     #logit("Create empty raster of (%f X %f) with %d bands\n" % (tiff_width, tiff_height, tiff_nbands))
+     
+     # Set data type for all bands
+     for b in range(band_from, band_to):
+         out_band = out_ds.GetRasterBand(b)
+         band = ds.GetRasterBand(b)
+         assert band is not None, "Missing GDAL raster band %d" % b
+         logit("MSG: Band %d\n" % b)
+         out_band.DataType = band.DataType
+         
+     # Create raster GeoTransform
+     gt = ds.GetGeoTransform()
+     out_ds.SetGeoTransform(gt)
+      
+     # Create projection
+     proj = ds.GetProjection()
+     out_ds.SetProjection(proj)
+     
+     return out_ds
+ 
  
  def wkblify_raster_level(options, ds, level, band_range):
      assert ds is not None
***************
*** 874,879 ****
--- 971,981 ----
      # Write (original) raster to hex binary output
      tile_count = 0
      hexwkb = ''
+     
+     if options.register:
+         out_ds = create_empty_outdb_dataset(options, ds, band_range)
+     else:
+         out_ds = None
  
      for ycell in range(0, grid_size[1]):
          for xcell in range(0, grid_size[0]):
***************
*** 896,903 ****
                  assert band is not None, "Missing GDAL raster band %d" % b
                  logit("MSG: Band %d\n" % b)
  
!                 hexwkb += wkblify_band_header(options, band)
!                 hexwkb += wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size)
  
              # INSERT INTO
              check_hex(hexwkb) # TODO: Remove to not to decrease performance
--- 998,1005 ----
                  assert band is not None, "Missing GDAL raster band %d" % b
                  logit("MSG: Band %d\n" % b)
  
!                 hexwkb += wkblify_band_header(options, band)                
!                 hexwkb += wkblify_band(options, band, level, xoff, yoff, read_block_size, block_size, out_ds, b)
  
              # INSERT INTO
              check_hex(hexwkb) # TODO: Remove to not to decrease performance
***************
*** 930,936 ****
          band_range = ( 1, ds.RasterCount + 1 )
  
      # Generate requested overview level (base raster if level = 1)
!     summary = wkblify_raster_level(options, ds, options.overview_level, band_range)
      SUMMARY.append( summary )
      
      # Cleanup
--- 1032,1038 ----
          band_range = ( 1, ds.RasterCount + 1 )
  
      # Generate requested overview level (base raster if level = 1)
!     summary = wkblify_raster_level(options,ds, options.overview_level, band_range)
      SUMMARY.append( summary )
      
      # Cleanup