Changeset 4163 for spike/wktraster/scripts/gdal2wktraster.py

Timestamp:

06/10/09 07:24:31 (3 years ago)

Author:

mloskot

Message:

[wktraster] Fixed stupid bug that was exporting pixel data of non-overview band as it was overview. Simply, band.GetOverview?() call was missing. mloskot's writing 1000 times: don't listen to reggae musing while coding.

Files:

• spike/wktraster/scripts/gdal2wktraster.py (modified) (6 diffs)

spike/wktraster/scripts/gdal2wktraster.py

@@ -470 +470 @@
     assert overview < band.GetOverviewCount()
     
-    band_ov = band.GetOverview(overview)
-    assert band_ov is not None
-    
-    ovf = int(0.5 + ds.RasterXSize / float(band_ov.XSize))
+    ov_band = band.GetOverview(overview)
+    assert ov_band is not None
+    
+    ovf = int(0.5 + ds.RasterXSize / float(ov_band.XSize))
     logit('MSG: Overview factor = %d\n' % ovf)
 
@@ -567 +567 @@
     
     return (xgeo, ygeo)
+
+def calculate_geoxy_overview(gt, xy, level):
+
+    # Update pixel resolution according to overview level
+    newgt = ( gt[0], gt[1] * float(level), gt[2], gt[3], gt[4], gt[5] * float(level) )
+
+    return calculate_geoxy(newgt, xy)
 
 def calculate_bounding_box(ds):
@@ -861 +868 @@
 
         for nov in range(0, noverviews):
-            msg = 'MSG: ---------- OVERVIEW # %d ------------------------------\n' % nov
+            msg = 'MSG: ---------- OVERVIEW # %d ------------------------------\n' % (nov + 1)
             logit(msg)
             options.output.write('\n-- ' + msg)
             
-            band_ov = ds.GetRasterBand(band_from).GetOverview(nov)
-            assert band_ov is not None
+            ov_band = ds.GetRasterBand(band_from).GetOverview(nov)
+            assert ov_band is not None
             ov_factor = calculate_overview_factor(ds, nov)
 
@@ -877 +884 @@
             options.output.write(sql)
 
-            raster_size = ( band_ov.XSize, band_ov.YSize )
+            raster_size = ( ov_band.XSize, ov_band.YSize )
             grid_size = calculate_grid_size(raster_size, block_size)
             
             logit("MSG: Processing overview = %s using block = %s as grid = %s\n" % \
                   (str(raster_size), str(block_size), str(grid_size)))
-            
+
+            ov_band = None
             tile_count = 0
 
@@ -890 +898 @@
                     xoff = xcell * block_size[0]
                     yoff = ycell * block_size[1]
-                    ul = calculate_geoxy(gt, (xoff, yoff))
+                    ul = calculate_geoxy_overview(gt, (xoff, yoff), ov_factor)
 
                     logit("MSG: ---------- CELL %d x %d: (%d x %d) (%d x %d); UL: %s  --------------------\n" % \
@@ -902 +910 @@
                         logit("MSG: Band %d\n" % nband)
 
-                        hexwkb += wkblify_band_header(options, band)
-                        hexwkb += wkblify_band(options, band, xoff, yoff, block_size)
+                        ov_band = band.GetOverview(nov)
+                        assert ov_band is not None, "Missing GDAL overview %d for band %d" % (nov, nband)
+                        
+                        hexwkb += wkblify_band_header(options, ov_band)
+                        hexwkb += wkblify_band(options, ov_band, xoff, yoff, block_size)
 
                     # INSERT INTO